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:   MatShellGetContext(A, &user);
 16:   MatView(user->B, viewer);
 17:   return 0;
 18: }

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

 24:   MatShellGetContext(A, &user);
 25:   MatMult(user->B, X, Y);
 26:   return 0;
 27: }

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

 33:   MatShellGetContext(A, &user);
 34:   MatMultTranspose(user->B, X, Y);
 35:   return 0;
 36: }

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

 42:   MatShellGetContext(A, &user);
 43:   MatGetDiagonal(user->B, X);
 44:   return 0;
 45: }

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

 58:   PetscObjectGetType((PetscObject)A, &mattypename);
 59:   PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename);
 60:   VecDuplicate(X, &W1);
 61:   VecDuplicate(X, &W2);
 62:   MatScale(A, 31);
 63:   MatShift(A, 37);
 64:   MatDiagonalScale(A, X, Y);
 65:   MatScale(A, 41);
 66:   MatDiagonalScale(A, Y, Z);
 67:   MatComputeOperator(A, MATDENSE, &E);

 69:   MatView(E, viewer);
 70:   PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n");
 71:   MatMult(A, Z, W1);
 72:   MatMultTranspose(A, W1, W2);
 73:   VecView(W2, viewer);

 75:   PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n");
 76:   VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff);
 77:   VecSet(W1, -1.0);
 78:   MatMultAdd(A, W1, W1, W2);
 79:   VecView(W2, viewer);
 80:   VecAXPY(W2, -1.0, diff);
 81:   VecNorm(W2, NORM_2, &nrm);
 82: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 84: #endif

 86:   VecSet(W2, -1.0);
 87:   MatMultAdd(A, W1, W2, W2);
 88:   VecView(W2, viewer);
 89:   VecAXPY(W2, -1.0, diff);
 90:   VecNorm(W2, NORM_2, &nrm);
 91: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 93: #endif
 94:   VecDestroy(&diff);

 96:   PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n");
 97:   VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff);

 99:   VecSet(W1, -1.0);
100:   MatMultTransposeAdd(A, W1, W1, W2);
101:   VecView(W2, viewer);
102:   VecAXPY(W2, -1.0, diff);
103:   VecNorm(W2, NORM_2, &nrm);
104: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
106: #endif

108:   VecSet(W2, -1.0);
109:   MatMultTransposeAdd(A, W1, W2, W2);
110:   VecView(W2, viewer);
111:   VecAXPY(W2, -1.0, diff);
112:   VecNorm(W2, NORM_2, &nrm);
113: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
115: #endif
116:   VecDestroy(&diff);

118:   PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n");
119:   MatGetDiagonal(A, W2);
120:   VecView(W2, viewer);
121:   VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff);
122:   VecAXPY(diff, -1.0, W2);
123:   VecNorm(diff, NORM_2, &nrm);
125:   VecDestroy(&diff);

127:   /* MATSHELL does not support MatDiagonalSet after MatScale */
128:   if (strncmp(mattypename, "shell", 5)) {
129:     MatDiagonalSet(A, X, INSERT_VALUES);
130:     MatGetDiagonal(A, W1);
131:     VecView(W1, viewer);
132:   } else {
133:     PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n");
134:   }

136:   MatDestroy(&E);
137:   VecDestroy(&W1);
138:   VecDestroy(&W2);
139:   return 0;
140: }

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

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

172:   PetscNew(&user);
173:   user->B = A;

175:   MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S);
176:   MatSetUp(S);
177:   MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User);
178:   MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User);
179:   MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User);
180:   MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User);

182:   for (i = 0; i < 4; i++) MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i]);
183:   MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N);
184:   MatSetUp(N);

186:   TestMatrix(S, X, Y, Z);
187:   TestMatrix(A, X, Y, Z);
188:   TestMatrix(N, X, Y, Z);

190:   for (i = 0; i < 4; i++) MatDestroy(&D[i]);
191:   MatDestroy(&A);
192:   MatDestroy(&S);
193:   MatDestroy(&N);
194:   VecDestroy(&X);
195:   VecDestroy(&Y);
196:   VecDestroy(&Z);
197:   PetscFree(user);
198:   PetscFinalize();
199:   return 0;
200: }

202: /*TEST

204:    test:

206: TEST*/