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;

 16:   PetscInitialize(&argc, &argv, (char *)0, help);
 17:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 18:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 19:   PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL);

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

 24:   /* Create a aij matrix A */
 25:   M = N = m * n;
 26:   MatCreate(PETSC_COMM_WORLD, &A);
 27:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
 28:   MatSetType(A, MATAIJ);
 29:   MatSetFromOptions(A);
 30:   MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
 31:   MatSeqAIJSetPreallocation(A, 5, NULL);

 33:   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:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 42:     }
 43:     if (i < m - 1) {
 44:       J = Ii + n;
 45:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 46:     }
 47:     if (j > 0) {
 48:       J = Ii - 1;
 49:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 50:     }
 51:     if (j < n - 1) {
 52:       J = Ii + 1;
 53:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 54:     }
 55:     v = 4.0;
 56:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
 57:   }
 58:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 61:   /* Create a dense matrix B */
 62:   MatGetLocalSize(A, &am, &an);
 63:   MatCreate(PETSC_COMM_WORLD, &B);
 64:   MatSetSizes(B, an, PETSC_DECIDE, PETSC_DECIDE, M);
 65:   MatSetType(B, MATDENSE);
 66:   MatSeqDenseSetPreallocation(B, NULL);
 67:   MatMPIDenseSetPreallocation(B, NULL);
 68:   MatSetFromOptions(B);
 69:   MatSetRandom(B, r);
 70:   PetscRandomDestroy(&r);

 72:   /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */
 73:   MatCreate(PETSC_COMM_WORLD, &C);
 74:   MatSetType(C, MATDENSE);
 75:   MatSetSizes(C, am, PETSC_DECIDE, PETSC_DECIDE, N);
 76:   MatSetFromOptions(C);
 77:   MatSetUp(C);
 78:   MatZeroEntries(C);
 79:   MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C);
 80:   MatNorm(C, NORM_INFINITY, &norm);
 81:   MatDestroy(&C);

 83:   /* Test C = A*B (aij*dense) */
 84:   MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C);
 85:   MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C);

 87:   /* Test developer API */
 88:   MatProductCreate(A, B, NULL, &D);
 89:   MatProductSetType(D, MATPRODUCT_AB);
 90:   MatProductSetAlgorithm(D, "default");
 91:   MatProductSetFill(D, fill);
 92:   MatProductSetFromOptions(D);
 93:   MatProductSymbolic(D);
 94:   for (i = 0; i < 2; i++) MatProductNumeric(D);
 95:   MatEqual(C, D, &equal);
 97:   MatDestroy(&D);

 99:   /* Test D = AT*B (transpose(aij)*dense) */
100:   MatCreateTranspose(A, &AT);
101:   MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &D);
102:   MatMatMultEqual(AT, B, D, 10, &equal);
104:   MatDestroy(&D);
105:   MatDestroy(&AT);

107:   /* Test D = C*A (dense*aij) */
108:   MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D);
109:   MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D);
110:   MatMatMultEqual(C, A, D, 10, &equal);
112:   MatDestroy(&D);

114:   /* Test D = A*C (aij*dense) */
115:   MatMatMult(A, C, MAT_INITIAL_MATRIX, fill, &D);
116:   MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &D);
117:   MatMatMultEqual(A, C, D, 10, &equal);
119:   MatDestroy(&D);

121:   /* Test D = B*C (dense*dense) */
122:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
123:   if (size == 1) {
124:     MatMatMult(B, C, MAT_INITIAL_MATRIX, fill, &D);
125:     MatMatMult(B, C, MAT_REUSE_MATRIX, fill, &D);
126:     MatMatMultEqual(B, C, D, 10, &equal);
128:     MatDestroy(&D);
129:   }

131:   /* Test D = B*C^T (dense*dense) */
132:   MatMatTransposeMult(B, C, MAT_INITIAL_MATRIX, fill, &D);
133:   MatMatTransposeMult(B, C, MAT_REUSE_MATRIX, fill, &D);
134:   MatMatTransposeMultEqual(B, C, D, 10, &equal);
136:   MatDestroy(&D);

138:   /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */
139:   flg = PETSC_FALSE;
140:   PetscOptionsHasName(NULL, NULL, "-test_userAPI", &flg);
141:   if (flg) {
142:     /* user driver */
143:     MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &B);
144:   } else {
145:     /* clear internal data structures related with previous products to avoid circular references */
146:     MatProductClear(A);
147:     MatProductClear(B);
148:     MatProductClear(C);
149:     MatProductCreateWithMat(A, C, NULL, B);
150:     MatProductSetType(B, MATPRODUCT_AB);
151:     MatProductSetFromOptions(B);
152:     MatProductSymbolic(B);
153:     MatProductNumeric(B);
154:   }

156:   MatDestroy(&C);
157:   MatDestroy(&B);
158:   MatDestroy(&A);
159:   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*/