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