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