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 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 18 19 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 20 #if defined(PETSC_HAVE_HYPRE) 21 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL)); 22 #endif 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 25 26 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 27 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3)); 28 PetscCall(MatSetType(A,MATAIJ)); 29 PetscCall(MatSetFromOptions(A)); 30 PetscCall(MatSetUp(A)); 31 PetscCall(MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 32 33 if (rank == 0) { 34 PetscCall(MatSetValues(A,3,ij,3,ij,a,ADD_VALUES)); 35 } 36 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 37 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 38 39 /* Test MatMatMult() */ 40 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 41 PetscCall(MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C)); /* C = B*A */ 42 PetscCall(MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C)); /* recompute C=B*A */ 43 PetscCall(MatSetOptionsPrefix(C,"C_")); 44 PetscCall(MatMatMultEqual(B,A,C,10,&isequal)); 45 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: C != B*A"); 46 47 PetscCall(MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = C*A = (A^T*A)*A */ 48 PetscCall(MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D)); 49 PetscCall(MatMatMultEqual(C,A,D,10,&isequal)); 50 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: D != C*A"); 51 52 PetscCall(MatDestroy(&B)); 53 PetscCall(MatDestroy(&C)); 54 PetscCall(MatDestroy(&D)); 55 56 /* Test MatPtAP */ 57 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); /* B = A */ 58 PetscCall(MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C)); /* C = B^T*A*B */ 59 PetscCall(MatPtAPMultEqual(A,B,C,10,&isequal)); 60 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP: C != B^T*A*B"); 61 62 /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */ 63 PetscCall(MatPtAP(A,B,MAT_REUSE_MATRIX,fill,&C)); 64 PetscCall(MatPtAPMultEqual(A,B,C,10,&isequal)); 65 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP(reuse): C != B^T*A*B"); 66 67 PetscCall(MatDestroy(&C)); 68 69 /* Test MatPtAP with A as a dense matrix */ 70 { 71 Mat Adense; 72 PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense)); 73 PetscCall(MatPtAP(Adense,B,MAT_INITIAL_MATRIX,fill,&C)); 74 75 PetscCall(MatPtAPMultEqual(Adense,B,C,10,&isequal)); 76 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP(reuse): C != B^T*Adense*B"); 77 PetscCall(MatDestroy(&Adense)); 78 } 79 80 if (size == 1) { 81 /* A test contributed by Tobias Neckel <neckel@in.tum.de> */ 82 PetscCall(testPTAPRectangular()); 83 84 /* test MatMatTransposeMult(): A*B^T */ 85 PetscCall(MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = A*A^T */ 86 PetscCall(MatScale(A,2.0)); 87 PetscCall(MatMatTransposeMult(A,A,MAT_REUSE_MATRIX,fill,&D)); 88 PetscCall(MatMatTransposeMultEqual(A,A,D,10,&isequal)); 89 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTranspose: D != A*A^T"); 90 } 91 92 PetscCall(MatDestroy(&A)); 93 PetscCall(MatDestroy(&B)); 94 PetscCall(MatDestroy(&C)); 95 PetscCall(MatDestroy(&D)); 96 PetscCall(PetscFinalize()); 97 return 0; 98 } 99 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; 106 107 PetscFunctionBegin; 108 /* set up A */ 109 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A)); 110 for (i=0; i<rows; i++) { 111 PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES)); 112 } 113 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 114 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 115 116 /* set up P */ 117 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P)); 118 PetscCall(MatSetValue(P, 0, 0, 1.0, INSERT_VALUES)); 119 PetscCall(MatSetValue(P, 0, 1, 2.0, INSERT_VALUES)); 120 PetscCall(MatSetValue(P, 0, 2, 0.0, INSERT_VALUES)); 121 122 PetscCall(MatSetValue(P, 0, 3, -1.0, INSERT_VALUES)); 123 124 PetscCall(MatSetValue(P, 1, 0, 0.0, INSERT_VALUES)); 125 PetscCall(MatSetValue(P, 1, 1, -1.0, INSERT_VALUES)); 126 PetscCall(MatSetValue(P, 1, 2, 1.0, INSERT_VALUES)); 127 128 PetscCall(MatSetValue(P, 2, 0, 3.0, INSERT_VALUES)); 129 PetscCall(MatSetValue(P, 2, 1, 0.0, INSERT_VALUES)); 130 PetscCall(MatSetValue(P, 2, 2, -3.0, INSERT_VALUES)); 131 132 PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 133 PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 134 135 /* compute C */ 136 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C)); 137 138 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 139 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 140 141 /* compare results */ 142 /* 143 printf("C:\n"); 144 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 145 146 blitz::Array<double,2> actualC(cols, cols); 147 actualC = 0.0; 148 for (int i=0; i<cols; i++) { 149 for (int j=0; j<cols; j++) { 150 PetscCall(MatGetValues(C, 1, &i, 1, &j, &actualC(i,j))); 151 } 152 } 153 blitz::Array<double,2> expectedC(cols, cols); 154 expectedC = 0.0; 155 156 expectedC(0,0) = 10.0; 157 expectedC(0,1) = 2.0; 158 expectedC(0,2) = -9.0; 159 expectedC(0,3) = -1.0; 160 expectedC(1,0) = 2.0; 161 expectedC(1,1) = 5.0; 162 expectedC(1,2) = -1.0; 163 expectedC(1,3) = -2.0; 164 expectedC(2,0) = -9.0; 165 expectedC(2,1) = -1.0; 166 expectedC(2,2) = 10.0; 167 expectedC(2,3) = 0.0; 168 expectedC(3,0) = -1.0; 169 expectedC(3,1) = -2.0; 170 expectedC(3,2) = 0.0; 171 expectedC(3,3) = 1.0; 172 173 int check = areBlitzArrays2NumericallyEqual(actualC,expectedC); 174 validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check)); 175 */ 176 177 PetscCall(MatDestroy(&A)); 178 PetscCall(MatDestroy(&P)); 179 PetscCall(MatDestroy(&C)); 180 PetscFunctionReturn(0); 181 } 182 183 /*TEST 184 185 test: 186 187 test: 188 suffix: 2 189 nsize: 2 190 args: -matmatmult_via nonscalable 191 output_file: output/ex93_1.out 192 193 test: 194 suffix: 3 195 nsize: 2 196 output_file: output/ex93_1.out 197 198 test: 199 suffix: 4 200 nsize: 2 201 args: -matptap_via scalable 202 output_file: output/ex93_1.out 203 204 test: 205 suffix: btheap 206 args: -matmatmult_via btheap -matmattransmult_via color 207 output_file: output/ex93_1.out 208 209 test: 210 suffix: heap 211 args: -matmatmult_via heap 212 output_file: output/ex93_1.out 213 214 #HYPRE PtAP is broken for complex numbers 215 test: 216 suffix: hypre 217 nsize: 3 218 requires: hypre !complex 219 args: -matmatmult_via hypre -matptap_via hypre -test_hypre 220 output_file: output/ex93_hypre.out 221 222 test: 223 suffix: llcondensed 224 args: -matmatmult_via llcondensed 225 output_file: output/ex93_1.out 226 227 test: 228 suffix: scalable 229 args: -matmatmult_via scalable 230 output_file: output/ex93_1.out 231 232 test: 233 suffix: scalable_fast 234 args: -matmatmult_via scalable_fast 235 output_file: output/ex93_1.out 236 237 TEST*/ 238