1 static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n"; 2 3 #include <petscmat.h> 4 5 extern PetscErrorCode testPTAPRectangular(void); 6 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 PetscScalar none=-1.; 13 PetscErrorCode ierr; 14 PetscReal fill=4; 15 PetscReal norm; 16 PetscMPIInt size,rank; 17 PetscBool test_hypre=PETSC_FALSE; 18 19 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 20 #if defined(PETSC_HAVE_HYPRE) 21 ierr = PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL);CHKERRQ(ierr); 22 #endif 23 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 24 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 25 26 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 27 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);CHKERRQ(ierr); 28 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 29 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 30 ierr = MatSetUp(A);CHKERRQ(ierr); 31 ierr = MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 32 33 if (!rank) { 34 ierr = MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);CHKERRQ(ierr); 35 } 36 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 37 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 38 ierr = MatSetOptionsPrefix(A,"A_");CHKERRQ(ierr); 39 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 40 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); 41 42 /* Test MatMatMult() */ 43 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); /* B = A^T */ 44 ierr = MatSetOptionsPrefix(B,"B_");CHKERRQ(ierr); 45 ierr = MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); /* C = B*A */ 46 ierr = MatMatMultNumeric(B,A,C);CHKERRQ(ierr); /* recompute C=B*A */ 47 ierr = MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); /* recompute C=B*A */ 48 ierr = MatSetOptionsPrefix(C,"C_");CHKERRQ(ierr); 49 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 50 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);} 51 52 ierr = MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);CHKERRQ(ierr); 53 ierr = MatMatMultNumeric(C,A,D);CHKERRQ(ierr); /* D = C*A = (A^T*A)*A */ 54 ierr = MatSetOptionsPrefix(D,"D_");CHKERRQ(ierr); 55 ierr = MatView(D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 56 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);} 57 58 /* Repeat the numeric product to test reuse of the previous symbolic product */ 59 ierr = MatMatMultNumeric(C,A,D);CHKERRQ(ierr); 60 ierr = MatView(D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 61 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);} 62 63 ierr = MatDestroy(&B);CHKERRQ(ierr); 64 ierr = MatDestroy(&C);CHKERRQ(ierr); 65 66 /* Test PtAP routine. */ 67 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); /* B = A */ 68 ierr = MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); /* C = B^T*A*B */ 69 if (!test_hypre) { 70 ierr = MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 71 ierr = MatNorm(D,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 72 if (norm > 1.e-15 && !rank) { 73 ierr = PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);CHKERRQ(ierr); 74 } 75 } 76 ierr = MatDestroy(&C);CHKERRQ(ierr); 77 ierr = MatDestroy(&D);CHKERRQ(ierr); 78 79 /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */ 80 ierr = MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 81 ierr = MatSetOptionsPrefix(C,"C=BtAB_");CHKERRQ(ierr); 82 if (!test_hypre) { 83 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 84 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); 85 } 86 87 if (!test_hypre) { 88 ierr = MatPtAPSymbolic(A,B,fill,&D);CHKERRQ(ierr); 89 ierr = MatPtAPNumeric(A,B,D);CHKERRQ(ierr); 90 ierr = MatSetOptionsPrefix(D,"D=BtAB_");CHKERRQ(ierr); 91 ierr = MatView(D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 92 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); 93 94 /* Repeat numeric product to test reuse of the previous symbolic product */ 95 ierr = MatPtAPNumeric(A,B,D);CHKERRQ(ierr); 96 ierr = MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 97 ierr = MatNorm(D,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 98 if (norm > 1.e-15) { 99 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);CHKERRQ(ierr); 100 } 101 ierr = MatDestroy(&D);CHKERRQ(ierr); 102 } 103 ierr = MatDestroy(&C);CHKERRQ(ierr); 104 ierr = MatDestroy(&B);CHKERRQ(ierr); 105 106 if (size == 1) { 107 /* A test contributed by Tobias Neckel <neckel@in.tum.de> */ 108 ierr = testPTAPRectangular();CHKERRQ(ierr); 109 110 /* test MatMatTransposeMult(): A*B^T */ 111 ierr = MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D);CHKERRQ(ierr); /* D = A*A^T */ 112 ierr = MatSetOptionsPrefix(D,"D=A*A^T_");CHKERRQ(ierr); 113 ierr = MatView(D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 114 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); 115 116 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); /* B = A^T */ 117 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); /* C=A*B */ 118 ierr = MatSetOptionsPrefix(C,"D=A*B=A*A^T_");CHKERRQ(ierr); 119 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 120 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); 121 } 122 123 ierr = MatDestroy(&A);CHKERRQ(ierr); 124 ierr = MatDestroy(&B);CHKERRQ(ierr); 125 ierr = MatDestroy(&C);CHKERRQ(ierr); 126 ierr = MatDestroy(&D);CHKERRQ(ierr); 127 ierr = PetscFinalize(); 128 return ierr; 129 } 130 131 /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */ 132 PetscErrorCode testPTAPRectangular(void) 133 { 134 135 const int rows = 3; 136 const int cols = 5; 137 PetscErrorCode ierr; 138 int i; 139 Mat A; 140 Mat P; 141 Mat C; 142 143 PetscFunctionBegin; 144 /* set up A */ 145 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);CHKERRQ(ierr); 146 for (i=0; i<rows; i++) { 147 ierr = MatSetValue(A, i, i, 1.0, INSERT_VALUES);CHKERRQ(ierr); 148 } 149 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 152 /* set up P */ 153 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);CHKERRQ(ierr); 154 ierr = MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);CHKERRQ(ierr); 155 ierr = MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);CHKERRQ(ierr); 156 ierr = MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);CHKERRQ(ierr); 157 158 ierr = MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);CHKERRQ(ierr); 159 160 ierr = MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);CHKERRQ(ierr); 161 ierr = MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);CHKERRQ(ierr); 162 ierr = MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);CHKERRQ(ierr); 163 164 ierr = MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);CHKERRQ(ierr); 165 ierr = MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);CHKERRQ(ierr); 166 ierr = MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);CHKERRQ(ierr); 167 168 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 169 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 170 171 /* compute C */ 172 ierr = MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);CHKERRQ(ierr); 173 174 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 175 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 176 177 /* compare results */ 178 /* 179 printf("C:\n"); 180 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 181 182 blitz::Array<double,2> actualC(cols, cols); 183 actualC = 0.0; 184 for (int i=0; i<cols; i++) { 185 for (int j=0; j<cols; j++) { 186 ierr = MatGetValues(C, 1, &i, 1, &j, &actualC(i,j)); 187 CHKERRQ(ierr); ; 188 } 189 } 190 blitz::Array<double,2> expectedC(cols, cols); 191 expectedC = 0.0; 192 193 expectedC(0,0) = 10.0; 194 expectedC(0,1) = 2.0; 195 expectedC(0,2) = -9.0; 196 expectedC(0,3) = -1.0; 197 expectedC(1,0) = 2.0; 198 expectedC(1,1) = 5.0; 199 expectedC(1,2) = -1.0; 200 expectedC(1,3) = -2.0; 201 expectedC(2,0) = -9.0; 202 expectedC(2,1) = -1.0; 203 expectedC(2,2) = 10.0; 204 expectedC(2,3) = 0.0; 205 expectedC(3,0) = -1.0; 206 expectedC(3,1) = -2.0; 207 expectedC(3,2) = 0.0; 208 expectedC(3,3) = 1.0; 209 210 int check = areBlitzArrays2NumericallyEqual(actualC,expectedC); 211 validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check)); 212 */ 213 214 ierr = MatDestroy(&A);CHKERRQ(ierr); 215 ierr = MatDestroy(&P);CHKERRQ(ierr); 216 ierr = MatDestroy(&C);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 /*TEST 221 222 test: 223 224 test: 225 suffix: 2 226 nsize: 2 227 args: -B_matmatmult_via nonscalable 228 229 test: 230 suffix: 3 231 nsize: 2 232 output_file: output/ex93_2.out 233 234 test: 235 suffix: 4 236 nsize: 2 237 args: -A_matptap_via scalable 238 output_file: output/ex93_2.out 239 240 test: 241 suffix: btheap 242 args: -B_matmatmult_via btheap 243 output_file: output/ex93_1.out 244 245 test: 246 suffix: heap 247 args: -B_matmatmult_via heap 248 output_file: output/ex93_1.out 249 250 #HYPRE PtAP is broken for complex numbers 251 test: 252 suffix: hypre 253 nsize: 3 254 requires: hypre !complex 255 args: -B_matmatmult_via hypre -A_matptap_via hypre -test_hypre 256 257 test: 258 suffix: llcondensed 259 args: -B_matmatmult_via llcondensed 260 output_file: output/ex93_1.out 261 262 test: 263 suffix: scalable 264 args: -B_matmatmult_via scalable 265 output_file: output/ex93_1.out 266 267 test: 268 suffix: scalable_fast 269 args: -B_matmatmult_via scalable_fast 270 output_file: output/ex93_1.out 271 272 TEST*/ 273