Actual source code: ex93.c

  1: static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";

  3: #include <petscmat.h>

  5: extern PetscErrorCode testPTAPRectangular(void);

  7: int main(int argc, char **argv)
  8: {
  9:   Mat         A, B, C, D;
 10:   PetscScalar a[]  = {1., 1., 0., 0., 1., 1., 0., 0., 1.};
 11:   PetscInt    ij[] = {0, 1, 2};
 12:   PetscReal   fill = 4.0;
 13:   PetscMPIInt size, rank;
 14:   PetscBool   isequal;
 15: #if defined(PETSC_HAVE_HYPRE)
 16:   PetscBool test_hypre = PETSC_FALSE;
 17: #endif

 19:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 21: #if defined(PETSC_HAVE_HYPRE)
 22:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_hypre", &test_hypre, NULL));
 23: #endif
 24:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 25:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 27:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 28:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3));
 29:   PetscCall(MatSetType(A, MATAIJ));
 30:   PetscCall(MatSetFromOptions(A));
 31:   PetscCall(MatSetUp(A));
 32:   PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));

 34:   if (rank == 0) PetscCall(MatSetValues(A, 3, ij, 3, ij, a, ADD_VALUES));
 35:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 36:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 38:   /* Test MatMatMult() */
 39:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));        /* B = A^T */
 40:   PetscCall(MatMatMult(B, A, MAT_INITIAL_MATRIX, fill, &C)); /* C = B*A */
 41:   PetscCall(MatMatMult(B, A, MAT_REUSE_MATRIX, fill, &C));   /* recompute C=B*A */
 42:   PetscCall(MatSetOptionsPrefix(C, "C_"));
 43:   PetscCall(MatMatMultEqual(B, A, C, 10, &isequal));
 44:   PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatMult: C != B*A");

 46:   PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D)); /* D = C*A = (A^T*A)*A */
 47:   PetscCall(MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D));
 48:   PetscCall(MatMatMultEqual(C, A, D, 10, &isequal));
 49:   PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatMult: D != C*A");

 51:   PetscCall(MatDestroy(&B));
 52:   PetscCall(MatDestroy(&C));
 53:   PetscCall(MatDestroy(&D));

 55:   /* Test MatPtAP */
 56:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));        /* B = A */
 57:   PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, fill, &C)); /* C = B^T*A*B */
 58:   PetscCall(MatPtAPMultEqual(A, B, C, 10, &isequal));
 59:   PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP: C != B^T*A*B");

 61:   /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
 62:   PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, fill, &C));
 63:   PetscCall(MatPtAPMultEqual(A, B, C, 10, &isequal));
 64:   PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP(reuse): C != B^T*A*B");

 66:   PetscCall(MatDestroy(&C));

 68:   /* Test MatPtAP with A as a dense matrix */
 69:   {
 70:     Mat Adense;
 71:     PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense));
 72:     PetscCall(MatPtAP(Adense, B, MAT_INITIAL_MATRIX, fill, &C));

 74:     PetscCall(MatPtAPMultEqual(Adense, B, C, 10, &isequal));
 75:     PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatPtAP(reuse): C != B^T*Adense*B");
 76:     PetscCall(MatDestroy(&Adense));
 77:   }

 79:   if (size == 1) {
 80:     /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
 81:     PetscCall(testPTAPRectangular());

 83:     /* test MatMatTransposeMult(): A*B^T */
 84:     PetscCall(MatMatTransposeMult(A, A, MAT_INITIAL_MATRIX, fill, &D)); /* D = A*A^T */
 85:     PetscCall(MatScale(A, 2.0));
 86:     PetscCall(MatMatTransposeMult(A, A, MAT_REUSE_MATRIX, fill, &D));
 87:     PetscCall(MatMatTransposeMultEqual(A, A, D, 10, &isequal));
 88:     PetscCheck(isequal, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTranspose: D != A*A^T");
 89:   }

 91:   PetscCall(MatDestroy(&A));
 92:   PetscCall(MatDestroy(&B));
 93:   PetscCall(MatDestroy(&C));
 94:   PetscCall(MatDestroy(&D));
 95:   PetscCall(PetscFinalize());
 96:   return 0;
 97: }

 99: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
100: PetscErrorCode testPTAPRectangular(void)
101: {
102:   const int rows = 3, cols = 5;
103:   int       i;
104:   Mat       A, P, C;

106:   PetscFunctionBegin;
107:   /* set up A  */
108:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows, 1, NULL, &A));
109:   for (i = 0; i < rows; i++) PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
110:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
111:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

113:   /* set up P */
114:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols, 5, NULL, &P));
115:   PetscCall(MatSetValue(P, 0, 0, 1.0, INSERT_VALUES));
116:   PetscCall(MatSetValue(P, 0, 1, 2.0, INSERT_VALUES));
117:   PetscCall(MatSetValue(P, 0, 2, 0.0, INSERT_VALUES));

119:   PetscCall(MatSetValue(P, 0, 3, -1.0, INSERT_VALUES));

121:   PetscCall(MatSetValue(P, 1, 0, 0.0, INSERT_VALUES));
122:   PetscCall(MatSetValue(P, 1, 1, -1.0, INSERT_VALUES));
123:   PetscCall(MatSetValue(P, 1, 2, 1.0, INSERT_VALUES));

125:   PetscCall(MatSetValue(P, 2, 0, 3.0, INSERT_VALUES));
126:   PetscCall(MatSetValue(P, 2, 1, 0.0, INSERT_VALUES));
127:   PetscCall(MatSetValue(P, 2, 2, -3.0, INSERT_VALUES));

129:   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
130:   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));

132:   /* compute C */
133:   PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C));

135:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
136:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

138:   /* compare results */
139:   /*
140:   printf("C:\n");
141:   PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD));

143:   blitz::Array<double,2> actualC(cols, cols);
144:   actualC = 0.0;
145:   for (int i=0; i<cols; i++) {
146:     for (int j=0; j<cols; j++) {
147:       PetscCall(MatGetValues(C, 1, &i, 1, &j, &actualC(i,j)));
148:     }
149:   }
150:   blitz::Array<double,2> expectedC(cols, cols);
151:   expectedC = 0.0;

153:   expectedC(0,0) = 10.0;
154:   expectedC(0,1) =  2.0;
155:   expectedC(0,2) = -9.0;
156:   expectedC(0,3) = -1.0;
157:   expectedC(1,0) =  2.0;
158:   expectedC(1,1) =  5.0;
159:   expectedC(1,2) = -1.0;
160:   expectedC(1,3) = -2.0;
161:   expectedC(2,0) = -9.0;
162:   expectedC(2,1) = -1.0;
163:   expectedC(2,2) = 10.0;
164:   expectedC(2,3) =  0.0;
165:   expectedC(3,0) = -1.0;
166:   expectedC(3,1) = -2.0;
167:   expectedC(3,2) =  0.0;
168:   expectedC(3,3) =  1.0;

170:   int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
171:   validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
172:   */

174:   PetscCall(MatDestroy(&A));
175:   PetscCall(MatDestroy(&P));
176:   PetscCall(MatDestroy(&C));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*TEST

182:    test:

184:    test:
185:       suffix: 2
186:       nsize: 2
187:       args: -matmatmult_via nonscalable
188:       output_file: output/ex93_1.out

190:    test:
191:       suffix: 3
192:       nsize: 2
193:       output_file: output/ex93_1.out

195:    test:
196:       suffix: 4
197:       nsize: 2
198:       args: -matptap_via scalable
199:       output_file: output/ex93_1.out

201:    test:
202:       suffix: btheap
203:       args: -matmatmult_via btheap -matmattransmult_via color
204:       output_file: output/ex93_1.out

206:    test:
207:       suffix: heap
208:       args: -matmatmult_via heap
209:       output_file: output/ex93_1.out

211:    #HYPRE PtAP is broken for complex numbers
212:    test:
213:       suffix: hypre
214:       nsize: 3
215:       requires: hypre !complex
216:       args: -matmatmult_via hypre -matptap_via hypre -test_hypre
217:       output_file: output/ex93_hypre.out

219:    test:
220:       suffix: llcondensed
221:       args: -matmatmult_via llcondensed
222:       output_file: output/ex93_1.out

224:    test:
225:       suffix: scalable
226:       args: -matmatmult_via scalable
227:       output_file: output/ex93_1.out

229:    test:
230:       suffix: scalable_fast
231:       args: -matmatmult_via scalable_fast
232:       output_file: output/ex93_1.out

234: TEST*/