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) { 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 if (!isequal) SETERRQ(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 if (!isequal) SETERRQ(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 if (!isequal) SETERRQ(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 if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP(reuse): C != B^T*A*B"); 67 68 ierr = MatDestroy(&C);CHKERRQ(ierr); 69 ierr = MatDestroy(&B);CHKERRQ(ierr); 70 71 if (size == 1) { 72 /* A test contributed by Tobias Neckel <neckel@in.tum.de> */ 73 ierr = testPTAPRectangular();CHKERRQ(ierr); 74 75 /* test MatMatTransposeMult(): A*B^T */ 76 ierr = MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D);CHKERRQ(ierr); /* D = A*A^T */ 77 ierr = MatScale(A,2.0);CHKERRQ(ierr); 78 ierr = MatMatTransposeMult(A,A,MAT_REUSE_MATRIX,fill,&D);CHKERRQ(ierr); 79 ierr = MatMatTransposeMultEqual(A,A,D,10,&isequal);CHKERRQ(ierr); 80 if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTranspose: D != A*A^T"); 81 } 82 83 ierr = MatDestroy(&A);CHKERRQ(ierr); 84 ierr = MatDestroy(&B);CHKERRQ(ierr); 85 ierr = MatDestroy(&C);CHKERRQ(ierr); 86 ierr = MatDestroy(&D);CHKERRQ(ierr); 87 ierr = PetscFinalize(); 88 return ierr; 89 } 90 91 /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */ 92 PetscErrorCode testPTAPRectangular(void) 93 { 94 const int rows = 3,cols = 5; 95 PetscErrorCode ierr; 96 int i; 97 Mat A,P,C; 98 99 PetscFunctionBegin; 100 /* set up A */ 101 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);CHKERRQ(ierr); 102 for (i=0; i<rows; i++) { 103 ierr = MatSetValue(A, i, i, 1.0, INSERT_VALUES);CHKERRQ(ierr); 104 } 105 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 106 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 108 /* set up P */ 109 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);CHKERRQ(ierr); 110 ierr = MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);CHKERRQ(ierr); 111 ierr = MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);CHKERRQ(ierr); 112 ierr = MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);CHKERRQ(ierr); 113 114 ierr = MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);CHKERRQ(ierr); 115 116 ierr = MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);CHKERRQ(ierr); 117 ierr = MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);CHKERRQ(ierr); 118 ierr = MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);CHKERRQ(ierr); 119 120 ierr = MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);CHKERRQ(ierr); 121 ierr = MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);CHKERRQ(ierr); 122 ierr = MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);CHKERRQ(ierr); 123 124 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 125 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 126 127 /* compute C */ 128 ierr = MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);CHKERRQ(ierr); 129 130 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 131 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 132 133 /* compare results */ 134 /* 135 printf("C:\n"); 136 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 137 138 blitz::Array<double,2> actualC(cols, cols); 139 actualC = 0.0; 140 for (int i=0; i<cols; i++) { 141 for (int j=0; j<cols; j++) { 142 ierr = MatGetValues(C, 1, &i, 1, &j, &actualC(i,j)); 143 CHKERRQ(ierr); ; 144 } 145 } 146 blitz::Array<double,2> expectedC(cols, cols); 147 expectedC = 0.0; 148 149 expectedC(0,0) = 10.0; 150 expectedC(0,1) = 2.0; 151 expectedC(0,2) = -9.0; 152 expectedC(0,3) = -1.0; 153 expectedC(1,0) = 2.0; 154 expectedC(1,1) = 5.0; 155 expectedC(1,2) = -1.0; 156 expectedC(1,3) = -2.0; 157 expectedC(2,0) = -9.0; 158 expectedC(2,1) = -1.0; 159 expectedC(2,2) = 10.0; 160 expectedC(2,3) = 0.0; 161 expectedC(3,0) = -1.0; 162 expectedC(3,1) = -2.0; 163 expectedC(3,2) = 0.0; 164 expectedC(3,3) = 1.0; 165 166 int check = areBlitzArrays2NumericallyEqual(actualC,expectedC); 167 validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check)); 168 */ 169 170 ierr = MatDestroy(&A);CHKERRQ(ierr); 171 ierr = MatDestroy(&P);CHKERRQ(ierr); 172 ierr = MatDestroy(&C);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 /*TEST 177 178 test: 179 180 test: 181 suffix: 2 182 nsize: 2 183 args: -matmatmult_via nonscalable 184 output_file: output/ex93_1.out 185 186 test: 187 suffix: 3 188 nsize: 2 189 output_file: output/ex93_1.out 190 191 test: 192 suffix: 4 193 nsize: 2 194 args: -matptap_via scalable 195 output_file: output/ex93_1.out 196 197 test: 198 suffix: btheap 199 args: -matmatmult_via btheap -matmattransmult_via color 200 output_file: output/ex93_1.out 201 202 test: 203 suffix: heap 204 args: -matmatmult_via heap 205 output_file: output/ex93_1.out 206 207 #HYPRE PtAP is broken for complex numbers 208 test: 209 suffix: hypre 210 nsize: 3 211 requires: hypre !complex 212 args: -matmatmult_via hypre -matptap_via hypre -test_hypre 213 output_file: output/ex93_hypre.out 214 215 test: 216 suffix: llcondensed 217 args: -matmatmult_via llcondensed 218 output_file: output/ex93_1.out 219 220 test: 221 suffix: scalable 222 args: -matmatmult_via scalable 223 output_file: output/ex93_1.out 224 225 test: 226 suffix: scalable_fast 227 args: -matmatmult_via scalable_fast 228 output_file: output/ex93_1.out 229 230 TEST*/ 231