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

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

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

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

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

 46:   MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D); /* D = C*A = (A^T*A)*A */
 47:   MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D);
 48:   MatMatMultEqual(C, A, D, 10, &isequal);

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

 55:   /* Test MatPtAP */
 56:   MatDuplicate(A, MAT_COPY_VALUES, &B);        /* B = A */
 57:   MatPtAP(A, B, MAT_INITIAL_MATRIX, fill, &C); /* C = B^T*A*B */
 58:   MatPtAPMultEqual(A, B, C, 10, &isequal);

 61:   /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
 62:   MatPtAP(A, B, MAT_REUSE_MATRIX, fill, &C);
 63:   MatPtAPMultEqual(A, B, C, 10, &isequal);

 66:   MatDestroy(&C);

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

 74:     MatPtAPMultEqual(Adense, B, C, 10, &isequal);
 76:     MatDestroy(&Adense);
 77:   }

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

 83:     /* test MatMatTransposeMult(): A*B^T */
 84:     MatMatTransposeMult(A, A, MAT_INITIAL_MATRIX, fill, &D); /* D = A*A^T */
 85:     MatScale(A, 2.0);
 86:     MatMatTransposeMult(A, A, MAT_REUSE_MATRIX, fill, &D);
 87:     MatMatTransposeMultEqual(A, A, D, 10, &isequal);
 89:   }

 91:   MatDestroy(&A);
 92:   MatDestroy(&B);
 93:   MatDestroy(&C);
 94:   MatDestroy(&D);
 95:   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:   /* set up A  */
107:   MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows, 1, NULL, &A);
108:   for (i = 0; i < rows; i++) MatSetValue(A, i, i, 1.0, INSERT_VALUES);
109:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
110:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

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

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

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

128:   MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
129:   MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);

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

134:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
135:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

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

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

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

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

173:   MatDestroy(&A);
174:   MatDestroy(&P);
175:   MatDestroy(&C);
176:   return 0;
177: }

179: /*TEST

181:    test:

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

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

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

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

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

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

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

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

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

233: TEST*/