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 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));CHKERRQ(ierr); 143 } 144 } 145 blitz::Array<double,2> expectedC(cols, cols); 146 expectedC = 0.0; 147 148 expectedC(0,0) = 10.0; 149 expectedC(0,1) = 2.0; 150 expectedC(0,2) = -9.0; 151 expectedC(0,3) = -1.0; 152 expectedC(1,0) = 2.0; 153 expectedC(1,1) = 5.0; 154 expectedC(1,2) = -1.0; 155 expectedC(1,3) = -2.0; 156 expectedC(2,0) = -9.0; 157 expectedC(2,1) = -1.0; 158 expectedC(2,2) = 10.0; 159 expectedC(2,3) = 0.0; 160 expectedC(3,0) = -1.0; 161 expectedC(3,1) = -2.0; 162 expectedC(3,2) = 0.0; 163 expectedC(3,3) = 1.0; 164 165 int check = areBlitzArrays2NumericallyEqual(actualC,expectedC); 166 validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check)); 167 */ 168 169 ierr = MatDestroy(&A);CHKERRQ(ierr); 170 ierr = MatDestroy(&P);CHKERRQ(ierr); 171 ierr = MatDestroy(&C);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 /*TEST 176 177 test: 178 179 test: 180 suffix: 2 181 nsize: 2 182 args: -matmatmult_via nonscalable 183 output_file: output/ex93_1.out 184 185 test: 186 suffix: 3 187 nsize: 2 188 output_file: output/ex93_1.out 189 190 test: 191 suffix: 4 192 nsize: 2 193 args: -matptap_via scalable 194 output_file: output/ex93_1.out 195 196 test: 197 suffix: btheap 198 args: -matmatmult_via btheap -matmattransmult_via color 199 output_file: output/ex93_1.out 200 201 test: 202 suffix: heap 203 args: -matmatmult_via heap 204 output_file: output/ex93_1.out 205 206 #HYPRE PtAP is broken for complex numbers 207 test: 208 suffix: hypre 209 nsize: 3 210 requires: hypre !complex 211 args: -matmatmult_via hypre -matptap_via hypre -test_hypre 212 output_file: output/ex93_hypre.out 213 214 test: 215 suffix: llcondensed 216 args: -matmatmult_via llcondensed 217 output_file: output/ex93_1.out 218 219 test: 220 suffix: scalable 221 args: -matmatmult_via scalable 222 output_file: output/ex93_1.out 223 224 test: 225 suffix: scalable_fast 226 args: -matmatmult_via scalable_fast 227 output_file: output/ex93_1.out 228 229 TEST*/ 230