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, NULL, 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 /* compare results */ 136 /* 137 printf("C:\n"); 138 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 139 140 blitz::Array<double,2> actualC(cols, cols); 141 actualC = 0.0; 142 for (int i=0; i<cols; i++) { 143 for (int j=0; j<cols; j++) PetscCall(MatGetValues(C, 1, &i, 1, &j, &actualC(i,j))); 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 PetscCall(MatDestroy(&A)); 170 PetscCall(MatDestroy(&P)); 171 PetscCall(MatDestroy(&C)); 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174 175 /*TEST 176 177 test: 178 output_file: output/empty.out 179 180 test: 181 suffix: 2 182 nsize: 2 183 args: -matmatmult_via nonscalable 184 output_file: output/empty.out 185 186 test: 187 suffix: 3 188 nsize: 2 189 output_file: output/empty.out 190 191 test: 192 suffix: 4 193 nsize: 2 194 args: -matptap_via scalable 195 output_file: output/empty.out 196 197 test: 198 suffix: btheap 199 args: -matmatmult_via btheap -matmattransmult_via color 200 output_file: output/empty.out 201 202 test: 203 suffix: heap 204 args: -matmatmult_via heap 205 output_file: output/empty.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/empty.out 214 215 test: 216 suffix: llcondensed 217 args: -matmatmult_via llcondensed 218 output_file: output/empty.out 219 220 test: 221 suffix: scalable 222 args: -matmatmult_via scalable 223 output_file: output/empty.out 224 225 test: 226 suffix: scalable_fast 227 args: -matmatmult_via scalable_fast 228 output_file: output/empty.out 229 230 TEST*/ 231