Actual source code: ex109.c

  1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";

  3: #include <petscmat.h>

  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;

 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));

 21:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
 22:   PetscCall(PetscRandomSetFromOptions(r));

 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));

 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));

 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));

 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));

 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));

 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));

 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));

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));

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));

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:   }

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));

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:   }

156:   PetscCall(MatDestroy(&C));
157:   PetscCall(MatDestroy(&B));
158:   PetscCall(MatDestroy(&A));
159:   PetscCall(PetscFinalize());
160:   return 0;
161: }

163: /*TEST

165:    test:
166:       args: -M 10 -N 10
167:       output_file: output/ex109.out

169:    test:
170:       suffix: 2
171:       nsize: 3
172:       output_file: output/ex109.out

174:    test:
175:       suffix: 3
176:       nsize: 2
177:       args: -matmattransmult_mpidense_mpidense_via cyclic
178:       output_file: output/ex109.out

180:    test:
181:       suffix: 4
182:       args: -test_userAPI
183:       output_file: output/ex109.out

185:    test:
186:       suffix: 5
187:       nsize: 3
188:       args: -test_userAPI
189:       output_file: output/ex109.out

191: TEST*/