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