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: };

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

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

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

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

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

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

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

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

 51: static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z)
 52: {
 53:   Vec         W1, W2, diff;
 54:   Mat         E;
 55:   const char *mattypename;
 56:   PetscViewer viewer      = PETSC_VIEWER_STDOUT_WORLD;
 57:   PetscScalar diag[2]     = {2.9678190300000000e+08, 1.4173141580000000e+09};
 58:   PetscScalar multadd[2]  = {-6.8966198500000000e+08, -2.0310609940000000e+09};
 59:   PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09};
 60:   PetscReal   nrm;

 62:   PetscFunctionBegin;
 63:   PetscCall(PetscObjectGetType((PetscObject)A, &mattypename));
 64:   PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename));
 65:   PetscCall(VecDuplicate(X, &W1));
 66:   PetscCall(VecDuplicate(X, &W2));
 67:   PetscCall(MatScale(A, 31));
 68:   PetscCall(MatShift(A, 37));
 69:   PetscCall(MatDiagonalScale(A, X, Y));
 70:   PetscCall(MatScale(A, 41));
 71:   PetscCall(MatDiagonalScale(A, Y, Z));
 72:   PetscCall(MatComputeOperator(A, MATDENSE, &E));

 74:   PetscCall(MatView(E, viewer));
 75:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n"));
 76:   PetscCall(MatMult(A, Z, W1));
 77:   PetscCall(MatMultTranspose(A, W1, W2));
 78:   PetscCall(VecView(W2, viewer));

 80:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n"));
 81:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff));
 82:   PetscCall(VecSet(W1, -1.0));
 83:   PetscCall(MatMultAdd(A, W1, W1, W2));
 84:   PetscCall(VecView(W2, viewer));
 85:   PetscCall(VecAXPY(W2, -1.0, diff));
 86:   PetscCall(VecNorm(W2, NORM_2, &nrm));
 87: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 88:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result");
 89: #endif

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

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

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

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

123:   PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n"));
124:   PetscCall(MatGetDiagonal(A, W2));
125:   PetscCall(VecView(W2, viewer));
126:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff));
127:   PetscCall(VecAXPY(diff, -1.0, W2));
128:   PetscCall(VecNorm(diff, NORM_2, &nrm));
129:   PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result");
130:   PetscCall(VecDestroy(&diff));

132:   /* MATSHELL does not support MatDiagonalSet after MatScale */
133:   if (strncmp(mattypename, "shell", 5)) {
134:     PetscCall(MatDiagonalSet(A, X, INSERT_VALUES));
135:     PetscCall(MatGetDiagonal(A, W1));
136:     PetscCall(VecView(W1, viewer));
137:   } else {
138:     PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n"));
139:   }

141:   PetscCall(MatDestroy(&E));
142:   PetscCall(VecDestroy(&W1));
143:   PetscCall(VecDestroy(&W2));
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

147: int main(int argc, char **args)
148: {
149:   const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29};
150:   const PetscInt    inds[]  = {0, 1};
151:   PetscScalar       avals[] = {2, 3, 5, 7};
152:   Mat               A, S, D[4], N;
153:   Vec               X, Y, Z;
154:   User              user;
155:   PetscInt          i;

157:   PetscFunctionBeginUser;
158:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
159:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A));
160:   PetscCall(MatSetUp(A));
161:   PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));
162:   PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
163:   PetscCall(VecDuplicate(X, &Y));
164:   PetscCall(VecDuplicate(X, &Z));
165:   PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
166:   PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES));
167:   PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES));
168:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
169:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
170:   PetscCall(VecAssemblyBegin(X));
171:   PetscCall(VecAssemblyBegin(Y));
172:   PetscCall(VecAssemblyBegin(Z));
173:   PetscCall(VecAssemblyEnd(X));
174:   PetscCall(VecAssemblyEnd(Y));
175:   PetscCall(VecAssemblyEnd(Z));

177:   PetscCall(PetscNew(&user));
178:   user->B = A;

180:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
181:   PetscCall(MatSetUp(S));
182:   PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User));
183:   PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User));
184:   PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User));
185:   PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User));

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

191:   PetscCall(TestMatrix(S, X, Y, Z));
192:   PetscCall(TestMatrix(A, X, Y, Z));
193:   PetscCall(TestMatrix(N, X, Y, Z));

195:   for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i]));
196:   PetscCall(MatDestroy(&A));
197:   PetscCall(MatDestroy(&S));
198:   PetscCall(MatDestroy(&N));
199:   PetscCall(VecDestroy(&X));
200:   PetscCall(VecDestroy(&Y));
201:   PetscCall(VecDestroy(&Z));
202:   PetscCall(PetscFree(user));
203:   PetscCall(PetscFinalize());
204:   return 0;
205: }

207: /*TEST

209:    test:

211: TEST*/