Actual source code: ex88.c
petsc-3.14.6 2021-03-30
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: };
11: static PetscErrorCode MatView_User(Mat A,PetscViewer viewer)
12: {
13: User user;
17: MatShellGetContext(A,&user);
18: MatView(user->B,viewer);
19: return(0);
20: }
22: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
23: {
24: User user;
28: MatShellGetContext(A,&user);
29: MatMult(user->B,X,Y);
30: return(0);
31: }
33: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
34: {
35: User user;
39: MatShellGetContext(A,&user);
40: MatMultTranspose(user->B,X,Y);
41: return(0);
42: }
44: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
45: {
46: User user;
50: MatShellGetContext(A,&user);
51: MatGetDiagonal(user->B,X);
52: return(0);
53: }
55: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
56: {
58: Vec W1,W2,diff;
59: Mat E;
60: const char *mattypename;
61: PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
62: PetscScalar diag[2] = { 2.9678190300000000e+08, 1.4173141580000000e+09};
63: PetscScalar multadd[2] = {-6.8966198500000000e+08,-2.0310609940000000e+09};
64: PetscScalar multtadd[2] = {-9.1052873900000000e+08,-1.8101942400000000e+09};
65: PetscReal nrm;
68: PetscObjectGetType((PetscObject)A,&mattypename);
69: PetscViewerASCIIPrintf(viewer,"\nMatrix of type: %s\n",mattypename);
70: VecDuplicate(X,&W1);
71: VecDuplicate(X,&W2);
72: MatScale(A,31);
73: MatShift(A,37);
74: MatDiagonalScale(A,X,Y);
75: MatScale(A,41);
76: MatDiagonalScale(A,Y,Z);
77: MatComputeOperator(A,MATDENSE,&E);
79: MatView(E,viewer);
80: PetscViewerASCIIPrintf(viewer,"Testing MatMult + MatMultTranspose\n");
81: MatMult(A,Z,W1);
82: MatMultTranspose(A,W1,W2);
83: VecView(W2,viewer);
85: PetscViewerASCIIPrintf(viewer,"Testing MatMultAdd\n");
86: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multadd,&diff);
87: VecSet(W1,-1.0);
88: MatMultAdd(A,W1,W1,W2);
89: VecView(W2,viewer);
90: VecAXPY(W2,-1.0,diff);
91: VecNorm(W2,NORM_2,&nrm);
92: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
93: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultAdd(A,x,x,y) produces incorrect result");
94: #endif
96: VecSet(W2,-1.0);
97: MatMultAdd(A,W1,W2,W2);
98: VecView(W2,viewer);
99: VecAXPY(W2,-1.0,diff);
100: VecNorm(W2,NORM_2,&nrm);
101: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
102: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultAdd(A,x,y,y) produces incorrect result");
103: #endif
104: VecDestroy(&diff);
106: PetscViewerASCIIPrintf(viewer,"Testing MatMultTranposeAdd\n");
107: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multtadd,&diff);
109: VecSet(W1,-1.0);
110: MatMultTransposeAdd(A,W1,W1,W2);
111: VecView(W2,viewer);
112: VecAXPY(W2,-1.0,diff);
113: VecNorm(W2,NORM_2,&nrm);
114: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
115: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultTranposeAdd(A,x,x,y) produces incorrect result");
116: #endif
118: VecSet(W2,-1.0);
119: MatMultTransposeAdd(A,W1,W2,W2);
120: VecView(W2,viewer);
121: VecAXPY(W2,-1.0,diff);
122: VecNorm(W2,NORM_2,&nrm);
123: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
124: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultTranposeAdd(A,x,y,y) produces incorrect result");
125: #endif
126: VecDestroy(&diff);
128: PetscViewerASCIIPrintf(viewer,"Testing MatGetDiagonal\n");
129: MatGetDiagonal(A,W2);
130: VecView(W2,viewer);
131: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,diag,&diff);
132: VecAXPY(diff,-1.0,W2);
133: VecNorm(diff,NORM_2,&nrm);
134: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetDiagonal() produces incorrect result");
135: VecDestroy(&diff);
137: /* MATSHELL does not support MatDiagonalSet after MatScale */
138: if (strncmp(mattypename, "shell", 5)) {
139: MatDiagonalSet(A,X,INSERT_VALUES);
140: MatGetDiagonal(A,W1);
141: VecView(W1,viewer);
142: } else {
143: PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n");
144: }
146: MatDestroy(&E);
147: VecDestroy(&W1);
148: VecDestroy(&W2);
149: return(0);
150: }
152: int main(int argc,char **args)
153: {
154: const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
155: const PetscInt inds[] = {0,1};
156: PetscScalar avals[] = {2,3,5,7};
157: Mat A,S,D[4],N;
158: Vec X,Y,Z;
159: User user;
160: PetscInt i;
161: PetscErrorCode ierr;
163: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
164: MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
165: MatSetUp(A);
166: MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
167: VecCreateSeq(PETSC_COMM_WORLD,2,&X);
168: VecDuplicate(X,&Y);
169: VecDuplicate(X,&Z);
170: VecSetValues(X,2,inds,xvals,INSERT_VALUES);
171: VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
172: VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
173: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
175: VecAssemblyBegin(X);
176: VecAssemblyBegin(Y);
177: VecAssemblyBegin(Z);
178: VecAssemblyEnd(X);
179: VecAssemblyEnd(Y);
180: VecAssemblyEnd(Z);
182: PetscNew(&user);
183: user->B = A;
185: MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
186: MatSetUp(S);
187: MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User);
188: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
189: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
190: MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
192: for (i=0; i<4; i++) {
193: MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
194: }
195: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N);
196: MatSetUp(N);
198: TestMatrix(S,X,Y,Z);
199: TestMatrix(A,X,Y,Z);
200: TestMatrix(N,X,Y,Z);
202: for (i=0; i<4; i++) {MatDestroy(&D[i]);}
203: MatDestroy(&A);
204: MatDestroy(&S);
205: MatDestroy(&N);
206: VecDestroy(&X);
207: VecDestroy(&Y);
208: VecDestroy(&Z);
209: PetscFree(user);
210: PetscFinalize();
211: return ierr;
212: }
215: /*TEST
217: test:
219: TEST*/