Actual source code: ex88.c
2: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";
4: #include <petscmat.h>
6: typedef struct _n_User *User;
7: struct _n_User {
8: Mat B;
9: };
13: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
14: {
15: User user;
19: MatShellGetContext(A,&user);
20: MatMult(user->B,X,Y);
21: return(0);
22: }
26: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
27: {
28: User user;
32: MatShellGetContext(A,&user);
33: MatMultTranspose(user->B,X,Y);
34: return(0);
35: }
39: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
40: {
41: User user;
45: MatShellGetContext(A,&user);
46: MatGetDiagonal(user->B,X);
47: return(0);
48: }
52: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
53: {
55: Vec W1,W2;
56: Mat E;
57: const char *mattypename;
58: PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
61: VecDuplicate(X,&W1);
62: VecDuplicate(X,&W2);
63: MatScale(A,31);
64: MatShift(A,37);
65: MatDiagonalScale(A,X,Y);
66: MatScale(A,41);
67: MatDiagonalScale(A,Y,Z);
68: MatComputeExplicitOperator(A,&E);
70: PetscObjectGetType((PetscObject)A,&mattypename);
71: PetscViewerASCIIPrintf(viewer,"Matrix of type: %s\n",mattypename);
72: MatView(E,viewer);
73: MatMult(A,Z,W1);
74: MatMultTranspose(A,W1,W2);
75: VecView(W2,viewer);
76: MatGetDiagonal(A,W2);
77: VecView(W2,viewer);
78: MatDestroy(&E);
79: VecDestroy(&W1);
80: VecDestroy(&W2);
81: return(0);
82: }
86: int main(int argc,char **args)
87: {
88: const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
89: const PetscInt inds[] = {0,1};
90: PetscScalar avals[] = {2,3,5,7};
91: Mat A,S,D[4],N;
92: Vec X,Y,Z;
93: User user;
94: PetscInt i;
97: PetscInitialize(&argc,&args,(char *)0,help);
98: MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,PETSC_NULL,&A);
99: MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
100: VecCreateSeq(PETSC_COMM_WORLD,2,&X);
101: VecDuplicate(X,&Y);
102: VecDuplicate(X,&Z);
103: VecSetValues(X,2,inds,xvals,INSERT_VALUES);
104: VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
105: VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
106: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108: VecAssemblyBegin(X);
109: VecAssemblyBegin(Y);
110: VecAssemblyBegin(Z);
111: VecAssemblyEnd(X);
112: VecAssemblyEnd(Y);
113: VecAssemblyEnd(Z);
115: PetscNew(struct _n_User,&user);
116: user->B = A;
118: MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
119: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_User);
120: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_User);
121: MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_User);
123: for (i=0; i<4; i++) {
124: MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
125: }
126: MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,D,&N);
128: TestMatrix(S,X,Y,Z);
129: TestMatrix(A,X,Y,Z);
130: TestMatrix(N,X,Y,Z);
132: for (i=0; i<4; i++) {MatDestroy(&D[i]);}
133: MatDestroy(&A);
134: MatDestroy(&S);
135: MatDestroy(&N);
136: VecDestroy(&X);
137: VecDestroy(&Y);
138: VecDestroy(&Z);
139: PetscFree(user);
140: PetscFinalize();
141: return 0;
142: }