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