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

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

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

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

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

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

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

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

 67:   MatDestroy(&C);

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

 75:     MatPtAPMultEqual(Adense,B,C,10,&isequal);
 77:     MatDestroy(&Adense);
 78:   }

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

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

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

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

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

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

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

123:   MatSetValue(P, 1, 0,  0.0, INSERT_VALUES);
124:   MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
125:   MatSetValue(P, 1, 2,  1.0, INSERT_VALUES);

127:   MatSetValue(P, 2, 0,  3.0, INSERT_VALUES);
128:   MatSetValue(P, 2, 1,  0.0, INSERT_VALUES);
129:   MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);

131:   MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
132:   MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

134:   /* compute C */
135:   MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);

137:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

140:   /* compare results */
141:   /*
142:   printf("C:\n");
143:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);

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

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

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

176:   MatDestroy(&A);
177:   MatDestroy(&P);
178:   MatDestroy(&C);
179:   return 0;
180: }

182: /*TEST

184:    test:

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

192:    test:
193:       suffix: 3
194:       nsize: 2
195:       output_file: output/ex93_1.out

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

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

208:    test:
209:       suffix: heap
210:       args: -matmatmult_via heap
211:       output_file: output/ex93_1.out

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

221:    test:
222:       suffix: llcondensed
223:       args: -matmatmult_via llcondensed
224:       output_file: output/ex93_1.out

226:    test:
227:       suffix: scalable
228:       args: -matmatmult_via scalable
229:       output_file: output/ex93_1.out

231:    test:
232:       suffix: scalable_fast
233:       args: -matmatmult_via scalable_fast
234:       output_file: output/ex93_1.out

236: TEST*/