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 PetscErrorCode ierr; 13 PetscReal fill=4.0; 14 PetscMPIInt size,rank; 15 PetscBool isequal; 16 #if defined(PETSC_HAVE_HYPRE) 17 PetscBool test_hypre=PETSC_FALSE; 18 #endif 19 20 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 21 #if defined(PETSC_HAVE_HYPRE) 22 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL)); 23 #endif 24 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 25 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 26 27 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 28 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3)); 29 CHKERRQ(MatSetType(A,MATAIJ)); 30 CHKERRQ(MatSetFromOptions(A)); 31 CHKERRQ(MatSetUp(A)); 32 CHKERRQ(MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 33 34 if (rank == 0) { 35 CHKERRQ(MatSetValues(A,3,ij,3,ij,a,ADD_VALUES)); 36 } 37 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 38 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 39 40 /* Test MatMatMult() */ 41 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 42 CHKERRQ(MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C)); /* C = B*A */ 43 CHKERRQ(MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C)); /* recompute C=B*A */ 44 CHKERRQ(MatSetOptionsPrefix(C,"C_")); 45 CHKERRQ(MatMatMultEqual(B,A,C,10,&isequal)); 46 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: C != B*A"); 47 48 CHKERRQ(MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = C*A = (A^T*A)*A */ 49 CHKERRQ(MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D)); 50 CHKERRQ(MatMatMultEqual(C,A,D,10,&isequal)); 51 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: D != C*A"); 52 53 CHKERRQ(MatDestroy(&B)); 54 CHKERRQ(MatDestroy(&C)); 55 CHKERRQ(MatDestroy(&D)); 56 57 /* Test MatPtAP */ 58 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&B)); /* B = A */ 59 CHKERRQ(MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C)); /* C = B^T*A*B */ 60 CHKERRQ(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 CHKERRQ(MatPtAP(A,B,MAT_REUSE_MATRIX,fill,&C)); 65 CHKERRQ(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 CHKERRQ(MatDestroy(&C)); 69 70 /* Test MatPtAP with A as a dense matrix */ 71 { 72 Mat Adense; 73 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense)); 74 CHKERRQ(MatPtAP(Adense,B,MAT_INITIAL_MATRIX,fill,&C)); 75 76 CHKERRQ(MatPtAPMultEqual(Adense,B,C,10,&isequal)); 77 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP(reuse): C != B^T*Adense*B"); 78 CHKERRQ(MatDestroy(&Adense)); 79 } 80 81 if (size == 1) { 82 /* A test contributed by Tobias Neckel <neckel@in.tum.de> */ 83 CHKERRQ(testPTAPRectangular()); 84 85 /* test MatMatTransposeMult(): A*B^T */ 86 CHKERRQ(MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = A*A^T */ 87 CHKERRQ(MatScale(A,2.0)); 88 CHKERRQ(MatMatTransposeMult(A,A,MAT_REUSE_MATRIX,fill,&D)); 89 CHKERRQ(MatMatTransposeMultEqual(A,A,D,10,&isequal)); 90 PetscCheck(isequal,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTranspose: D != A*A^T"); 91 } 92 93 CHKERRQ(MatDestroy(&A)); 94 CHKERRQ(MatDestroy(&B)); 95 CHKERRQ(MatDestroy(&C)); 96 CHKERRQ(MatDestroy(&D)); 97 ierr = PetscFinalize(); 98 return ierr; 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 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A)); 111 for (i=0; i<rows; i++) { 112 CHKERRQ(MatSetValue(A, i, i, 1.0, INSERT_VALUES)); 113 } 114 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 115 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 116 117 /* set up P */ 118 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P)); 119 CHKERRQ(MatSetValue(P, 0, 0, 1.0, INSERT_VALUES)); 120 CHKERRQ(MatSetValue(P, 0, 1, 2.0, INSERT_VALUES)); 121 CHKERRQ(MatSetValue(P, 0, 2, 0.0, INSERT_VALUES)); 122 123 CHKERRQ(MatSetValue(P, 0, 3, -1.0, INSERT_VALUES)); 124 125 CHKERRQ(MatSetValue(P, 1, 0, 0.0, INSERT_VALUES)); 126 CHKERRQ(MatSetValue(P, 1, 1, -1.0, INSERT_VALUES)); 127 CHKERRQ(MatSetValue(P, 1, 2, 1.0, INSERT_VALUES)); 128 129 CHKERRQ(MatSetValue(P, 2, 0, 3.0, INSERT_VALUES)); 130 CHKERRQ(MatSetValue(P, 2, 1, 0.0, INSERT_VALUES)); 131 CHKERRQ(MatSetValue(P, 2, 2, -3.0, INSERT_VALUES)); 132 133 CHKERRQ(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 134 CHKERRQ(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 135 136 /* compute C */ 137 CHKERRQ(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C)); 138 139 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 140 CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 141 142 /* compare results */ 143 /* 144 printf("C:\n"); 145 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatDestroy(&A)); 179 CHKERRQ(MatDestroy(&P)); 180 CHKERRQ(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