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