1 static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **argv) 6 { 7 Mat A, B, C, D, AT; 8 PetscInt i, M, N, Istart, Iend, n = 7, j, J, Ii, m = 8, am, an; 9 PetscScalar v; 10 PetscRandom r; 11 PetscBool equal = PETSC_FALSE, flg; 12 PetscReal fill = 1.0, norm; 13 PetscMPIInt size; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 17 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 19 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL)); 20 21 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r)); 22 PetscCall(PetscRandomSetFromOptions(r)); 23 24 /* Create a aij matrix A */ 25 M = N = m * n; 26 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 27 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 28 PetscCall(MatSetType(A, MATAIJ)); 29 PetscCall(MatSetFromOptions(A)); 30 PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL)); 31 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 32 33 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend)); 34 am = Iend - Istart; 35 for (Ii = Istart; Ii < Iend; Ii++) { 36 v = -1.0; 37 i = Ii / n; 38 j = Ii - i * n; 39 if (i > 0) { 40 J = Ii - n; 41 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 42 } 43 if (i < m - 1) { 44 J = Ii + n; 45 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 46 } 47 if (j > 0) { 48 J = Ii - 1; 49 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 50 } 51 if (j < n - 1) { 52 J = Ii + 1; 53 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 54 } 55 v = 4.0; 56 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 57 } 58 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 59 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 60 61 /* Create a dense matrix B */ 62 PetscCall(MatGetLocalSize(A, &am, &an)); 63 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 64 PetscCall(MatSetSizes(B, an, PETSC_DECIDE, PETSC_DECIDE, M)); 65 PetscCall(MatSetType(B, MATDENSE)); 66 PetscCall(MatSeqDenseSetPreallocation(B, NULL)); 67 PetscCall(MatMPIDenseSetPreallocation(B, NULL)); 68 PetscCall(MatSetFromOptions(B)); 69 PetscCall(MatSetRandom(B, r)); 70 PetscCall(PetscRandomDestroy(&r)); 71 72 /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */ 73 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 74 PetscCall(MatSetType(C, MATDENSE)); 75 PetscCall(MatSetSizes(C, am, PETSC_DECIDE, PETSC_DECIDE, N)); 76 PetscCall(MatSetFromOptions(C)); 77 PetscCall(MatSetUp(C)); 78 PetscCall(MatZeroEntries(C)); 79 PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C)); 80 PetscCall(MatNorm(C, NORM_INFINITY, &norm)); 81 PetscCall(MatDestroy(&C)); 82 83 /* Test C = A*B (aij*dense) */ 84 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C)); 85 PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C)); 86 87 /* Test developer API */ 88 PetscCall(MatProductCreate(A, B, NULL, &D)); 89 PetscCall(MatProductSetType(D, MATPRODUCT_AB)); 90 PetscCall(MatProductSetAlgorithm(D, "default")); 91 PetscCall(MatProductSetFill(D, fill)); 92 PetscCall(MatProductSetFromOptions(D)); 93 PetscCall(MatProductSymbolic(D)); 94 for (i = 0; i < 2; i++) PetscCall(MatProductNumeric(D)); 95 PetscCall(MatEqual(C, D, &equal)); 96 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "C != D"); 97 PetscCall(MatDestroy(&D)); 98 99 /* Test D = AT*B (transpose(aij)*dense) */ 100 PetscCall(MatCreateTranspose(A, &AT)); 101 PetscCall(MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &D)); 102 PetscCall(MatMatMultEqual(AT, B, D, 10, &equal)); 103 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != AT*B (transpose(aij)*dense)"); 104 PetscCall(MatDestroy(&D)); 105 PetscCall(MatDestroy(&AT)); 106 107 /* Test D = C*A (dense*aij) */ 108 PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D)); 109 PetscCall(MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D)); 110 PetscCall(MatMatMultEqual(C, A, D, 10, &equal)); 111 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != C*A (dense*aij)"); 112 PetscCall(MatDestroy(&D)); 113 114 /* Test D = A*C (aij*dense) */ 115 PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, fill, &D)); 116 PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &D)); 117 PetscCall(MatMatMultEqual(A, C, D, 10, &equal)); 118 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != A*C (aij*dense)"); 119 PetscCall(MatDestroy(&D)); 120 121 /* Test D = B*C (dense*dense) */ 122 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 123 if (size == 1) { 124 PetscCall(MatMatMult(B, C, MAT_INITIAL_MATRIX, fill, &D)); 125 PetscCall(MatMatMult(B, C, MAT_REUSE_MATRIX, fill, &D)); 126 PetscCall(MatMatMultEqual(B, C, D, 10, &equal)); 127 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != B*C (dense*dense)"); 128 PetscCall(MatDestroy(&D)); 129 } 130 131 /* Test D = B*C^T (dense*dense) */ 132 PetscCall(MatMatTransposeMult(B, C, MAT_INITIAL_MATRIX, fill, &D)); 133 PetscCall(MatMatTransposeMult(B, C, MAT_REUSE_MATRIX, fill, &D)); 134 PetscCall(MatMatTransposeMultEqual(B, C, D, 10, &equal)); 135 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "D != B*C^T (dense*dense)"); 136 PetscCall(MatDestroy(&D)); 137 138 /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */ 139 flg = PETSC_FALSE; 140 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_userAPI", &flg)); 141 if (flg) { 142 /* user driver */ 143 PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &B)); 144 } else { 145 /* clear internal data structures related with previous products to avoid circular references */ 146 PetscCall(MatProductClear(A)); 147 PetscCall(MatProductClear(B)); 148 PetscCall(MatProductClear(C)); 149 PetscCall(MatProductCreateWithMat(A, C, NULL, B)); 150 PetscCall(MatProductSetType(B, MATPRODUCT_AB)); 151 PetscCall(MatProductSetFromOptions(B)); 152 PetscCall(MatProductSymbolic(B)); 153 PetscCall(MatProductNumeric(B)); 154 } 155 156 PetscCall(MatDestroy(&C)); 157 PetscCall(MatDestroy(&B)); 158 PetscCall(MatDestroy(&A)); 159 PetscCall(PetscFinalize()); 160 return 0; 161 } 162 163 /*TEST 164 165 test: 166 args: -M 10 -N 10 167 output_file: output/ex109.out 168 169 test: 170 suffix: 2 171 nsize: 3 172 output_file: output/ex109.out 173 174 test: 175 suffix: 3 176 nsize: 2 177 args: -matmattransmult_mpidense_mpidense_via cyclic 178 output_file: output/ex109.out 179 180 test: 181 suffix: 4 182 args: -test_userAPI 183 output_file: output/ex109.out 184 185 test: 186 suffix: 5 187 nsize: 3 188 args: -test_userAPI 189 output_file: output/ex109.out 190 191 TEST*/ 192