Actual source code: ex88.c

  1: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";

  3: #include <petscmat.h>

  5: typedef struct _n_User *User;
  6: struct _n_User {
  7:   Mat B;
  8: };

 10: static PetscErrorCode MatView_User(Mat A, PetscViewer viewer)
 11: {
 12:   User user;

 14:   PetscFunctionBegin;
 15:   PetscCall(MatShellGetContext(A, &user));
 16:   PetscCall(MatView(user->B, viewer));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
 21: {
 22:   User user;

 24:   PetscFunctionBegin;
 25:   PetscCall(MatShellGetContext(A, &user));
 26:   PetscCall(MatMult(user->B, X, Y));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
 31: {
 32:   User user;

 34:   PetscFunctionBegin;
 35:   PetscCall(MatShellGetContext(A, &user));
 36:   PetscCall(MatMultTranspose(user->B, X, Y));
 37:   PetscFunctionReturn(PETSC_SUCCESS);
 38: }

 40: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
 41: {
 42:   User user;

 44:   PetscFunctionBegin;
 45:   PetscCall(MatShellGetContext(A, &user));
 46:   PetscCall(MatGetDiagonal(user->B, X));
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z)
 51: {
 52:   Vec         W1, W2, W3, diff;
 53:   Mat         E;
 54:   const char *mattypename;
 55:   PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
 56:   PetscReal   nrm;
 57: #if defined(PETSC_USE_COMPLEX)
 58:   const PetscScalar diag[2]     = {PetscCMPLX(-6.2902938000000000e+07, 4.5741953400000000e+08), PetscCMPLX(1.0828994620000000e+09, 1.2955916360000000e+09)};
 59:   const PetscScalar multadd[2]  = {PetscCMPLX(1.4926230300000000e+08, -1.2811063360000000e+09), PetscCMPLX(-1.2985220710000000e+09, -2.1029893020000000e+09)};
 60:   const PetscScalar multtadd[2] = {PetscCMPLX(-1.5271967100000000e+08, -1.2648172000000000e+09), PetscCMPLX(-9.9654009700000000e+08, -2.1192784380000000e+09)};
 61: #else
 62:   const PetscScalar diag[2]     = {2.9678190300000000e+08, 1.4173141580000000e+09};
 63:   const PetscScalar multadd[2]  = {-6.8966198500000000e+08, -2.0310609940000000e+09};
 64:   const PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09};
 65: #endif

 67:   PetscFunctionBegin;
 68:   PetscCall(PetscObjectGetType((PetscObject)A, &mattypename));
 69:   PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename));
 70:   PetscCall(VecDuplicate(X, &W1));
 71:   PetscCall(VecDuplicate(X, &W2));
 72:   PetscCall(VecDuplicate(X, &W3));
 73:   PetscCall(MatScale(A, 31));
 74:   PetscCall(MatShift(A, 37));
 75:   PetscCall(MatDiagonalScale(A, X, Y));
 76:   PetscCall(MatScale(A, 41));
 77:   PetscCall(MatDiagonalScale(A, Y, Z));
 78:   PetscCall(MatComputeOperator(A, MATDENSE, &E));

 80:   PetscCall(MatView(E, viewer));
 81:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n"));
 82:   PetscCall(MatMult(A, Z, W1));
 83:   PetscCall(MatMultTranspose(A, W1, W2));
 84:   PetscCall(VecView(W2, viewer));
 85:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTranspose\n"));
 86:   PetscCall(VecConjugate(W1));
 87:   PetscCall(MatMultHermitianTranspose(A, W1, W2));
 88:   PetscCall(VecConjugate(W2));
 89:   PetscCall(VecView(W2, viewer));

 91:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n"));
 92:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff));
 93:   PetscCall(VecSet(W1, -1.0));
 94:   PetscCall(MatMultAdd(A, W1, W1, W2));
 95:   PetscCall(VecView(W2, viewer));
 96:   PetscCall(VecAXPY(W2, -1.0, diff));
 97:   PetscCall(VecNorm(W2, NORM_2, &nrm));
 98: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 99:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result");
100: #endif

102:   PetscCall(VecSet(W2, -1.0));
103:   PetscCall(MatMultAdd(A, W1, W2, W2));
104:   PetscCall(VecView(W2, viewer));
105:   PetscCall(VecAXPY(W2, -1.0, diff));
106:   PetscCall(VecNorm(W2, NORM_2, &nrm));
107: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
108:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,y,y) produces incorrect result");
109: #endif
110:   PetscCall(VecDestroy(&diff));

112:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n"));
113:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));

115:   PetscCall(VecSet(W1, -1.0));
116:   PetscCall(MatMultTransposeAdd(A, W1, W1, W2));
117:   PetscCall(VecView(W2, viewer));
118:   PetscCall(VecAXPY(W2, -1.0, diff));
119:   PetscCall(VecNorm(W2, NORM_2, &nrm));
120: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
121:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,x,y) produces incorrect result");
122: #endif

124:   PetscCall(VecSet(W2, -1.0));
125:   PetscCall(MatMultTransposeAdd(A, W1, W2, W2));
126:   PetscCall(VecView(W2, viewer));
127:   PetscCall(VecAXPY(W2, -1.0, diff));
128:   PetscCall(VecNorm(W2, NORM_2, &nrm));
129: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
130:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,y,y) produces incorrect result");
131: #endif
132:   PetscCall(VecDestroy(&diff));

134:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTransposeAdd\n"));
135:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));

137:   PetscCall(VecSet(W1, -1.0));
138:   PetscCall(MatMultHermitianTransposeAdd(A, W1, W1, W3));
139:   PetscCall(VecConjugate(W3));
140:   PetscCall(VecView(W3, viewer));
141:   PetscCall(VecAXPY(W3, -1.0, diff));
142:   PetscCall(VecNorm(W3, NORM_2, &nrm));
143: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
144:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,x,y) produces incorrect result");
145: #endif

147:   PetscCall(VecSet(W3, -1.0));
148:   PetscCall(MatMultHermitianTransposeAdd(A, W1, W3, W3));
149:   PetscCall(VecConjugate(W3));
150:   PetscCall(VecView(W3, viewer));
151:   PetscCall(VecAXPY(W3, -1.0, diff));
152:   PetscCall(VecNorm(W3, NORM_2, &nrm));
153: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
154:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,y,y) produces incorrect result");
155: #endif
156:   PetscCall(VecDestroy(&diff));

158:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n"));
159:   PetscCall(MatGetDiagonal(A, W2));
160:   PetscCall(VecView(W2, viewer));
161:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff));
162:   PetscCall(VecAXPY(diff, -1.0, W2));
163:   PetscCall(VecNorm(diff, NORM_2, &nrm));
164: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
165:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result");
166: #endif
167:   PetscCall(VecDestroy(&diff));

169:   /* MATSHELL does not support MatDiagonalSet after MatScale */
170:   if (strncmp(mattypename, "shell", 5)) {
171:     PetscCall(MatDiagonalSet(A, X, INSERT_VALUES));
172:     PetscCall(MatGetDiagonal(A, W1));
173:     PetscCall(VecView(W1, viewer));
174:   } else {
175:     PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n"));
176:   }

178:   PetscCall(MatDestroy(&E));
179:   PetscCall(VecDestroy(&W1));
180:   PetscCall(VecDestroy(&W2));
181:   PetscCall(VecDestroy(&W3));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: int main(int argc, char **args)
186: {
187:   const PetscInt inds[] = {0, 1};
188: #if defined(PETSC_USE_COMPLEX)
189:   const PetscScalar xvals[] = {PetscCMPLX(11, 4), PetscCMPLX(13, 2)}, yvals[] = {PetscCMPLX(17, 3), PetscCMPLX(19, 1)}, zvals[] = {PetscCMPLX(23, 6), PetscCMPLX(29, 2)};
190:   PetscScalar       avals[] = {PetscCMPLX(2, 3), PetscCMPLX(3, 5), PetscCMPLX(5, 4), PetscCMPLX(7, 5)};
191: #else
192:   const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29};
193:   PetscScalar       avals[] = {2, 3, 5, 7};
194: #endif
195:   Mat      A, S, D[4], N;
196:   Vec      X, Y, Z;
197:   User     user;
198:   PetscInt i;

200:   PetscFunctionBeginUser;
201:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
202:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A));
203:   PetscCall(MatSetUp(A));
204:   PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
205:   PetscCall(VecDuplicate(X, &Y));
206:   PetscCall(VecDuplicate(X, &Z));
207:   PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));
208:   PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
209:   PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES));
210:   PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES));
211:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
212:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
213:   PetscCall(VecAssemblyBegin(X));
214:   PetscCall(VecAssemblyBegin(Y));
215:   PetscCall(VecAssemblyBegin(Z));
216:   PetscCall(VecAssemblyEnd(X));
217:   PetscCall(VecAssemblyEnd(Y));
218:   PetscCall(VecAssemblyEnd(Z));

220:   PetscCall(PetscNew(&user));
221:   user->B = A;

223:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
224:   PetscCall(MatSetUp(S));
225:   PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User));
226:   PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User));
227:   PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User));
228:   PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User));

230:   for (i = 0; i < 4; i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i]));
231:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N));
232:   PetscCall(MatSetUp(N));

234:   PetscCall(TestMatrix(S, X, Y, Z));
235:   PetscCall(TestMatrix(A, X, Y, Z));
236:   PetscCall(TestMatrix(N, X, Y, Z));

238:   for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i]));
239:   PetscCall(MatDestroy(&A));
240:   PetscCall(MatDestroy(&S));
241:   PetscCall(MatDestroy(&N));
242:   PetscCall(VecDestroy(&X));
243:   PetscCall(VecDestroy(&Y));
244:   PetscCall(VecDestroy(&Z));
245:   PetscCall(PetscFree(user));
246:   PetscCall(PetscFinalize());
247:   return 0;
248: }

250: /*TEST

252:    testset:
253:      test:
254:        suffix: 1
255:        requires:!complex
256:      test:
257:        suffix: 2
258:        requires: complex

260: TEST*/