Actual source code: ex89.c

  1: static char help[] = "Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";

  3: #include <petscdmda.h>

  5: int main(int argc, char **argv)
  6: {
  7:   DM              coarsedm, finedm;
  8:   PetscMPIInt     size, rank;
  9:   PetscInt        M, N, Z, i, nrows;
 10:   PetscScalar     one  = 1.0;
 11:   PetscReal       fill = 2.0;
 12:   Mat             A, P, C;
 13:   PetscScalar    *array, alpha;
 14:   PetscBool       Test_3D = PETSC_FALSE, flg;
 15:   const PetscInt *ia, *ja;
 16:   PetscInt        dof;
 17:   MPI_Comm        comm;

 20:   PetscInitialize(&argc, &argv, NULL, help);
 21:   comm = PETSC_COMM_WORLD;
 22:   MPI_Comm_rank(comm, &rank);
 23:   MPI_Comm_size(comm, &size);
 24:   M   = 10;
 25:   N   = 10;
 26:   Z   = 10;
 27:   dof = 10;

 29:   PetscOptionsGetBool(NULL, NULL, "-test_3D", &Test_3D, NULL);
 30:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 31:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
 32:   PetscOptionsGetInt(NULL, NULL, "-Z", &Z, NULL);
 33:   /* Set up distributed array for fine grid */
 34:   if (!Test_3D) {
 35:     DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &coarsedm);
 36:   } else {
 37:     DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, Z, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, NULL, &coarsedm);
 38:   }
 39:   DMSetFromOptions(coarsedm);
 40:   DMSetUp(coarsedm);

 42:   /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
 43:   DMRefine(coarsedm, PetscObjectComm((PetscObject)coarsedm), &finedm);

 45:   /*------------------------------------------------------------*/
 46:   DMSetMatType(finedm, MATAIJ);
 47:   DMCreateMatrix(finedm, &A);

 49:   /* set val=one to A */
 50:   if (size == 1) {
 51:     MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 52:     if (flg) {
 53:       MatSeqAIJGetArray(A, &array);
 54:       for (i = 0; i < ia[nrows]; i++) array[i] = one;
 55:       MatSeqAIJRestoreArray(A, &array);
 56:     }
 57:     MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 58:   } else {
 59:     Mat AA, AB;
 60:     MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL);
 61:     MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 62:     if (flg) {
 63:       MatSeqAIJGetArray(AA, &array);
 64:       for (i = 0; i < ia[nrows]; i++) array[i] = one;
 65:       MatSeqAIJRestoreArray(AA, &array);
 66:     }
 67:     MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 68:     MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 69:     if (flg) {
 70:       MatSeqAIJGetArray(AB, &array);
 71:       for (i = 0; i < ia[nrows]; i++) array[i] = one;
 72:       MatSeqAIJRestoreArray(AB, &array);
 73:     }
 74:     MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
 75:   }
 76:   /* Create interpolation between the fine and coarse grids */
 77:   DMCreateInterpolation(coarsedm, finedm, &P, NULL);

 79:   /* Test P^T * A * P - MatPtAP() */
 80:   /*------------------------------*/
 81:   /* (1) Developer API */
 82:   MatProductCreate(A, P, NULL, &C);
 83:   MatProductSetType(C, MATPRODUCT_PtAP);
 84:   MatProductSetAlgorithm(C, "allatonce");
 85:   MatProductSetFill(C, PETSC_DEFAULT);
 86:   MatProductSetFromOptions(C);
 87:   MatProductSymbolic(C);
 88:   MatProductNumeric(C);
 89:   MatProductNumeric(C); /* Test reuse of symbolic C */

 91:   { /* Test MatProductView() */
 92:     PetscViewer viewer;
 93:     PetscViewerASCIIOpen(comm, NULL, &viewer);
 94:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
 95:     MatProductView(C, viewer);
 96:     PetscViewerPopFormat(viewer);
 97:     PetscViewerDestroy(&viewer);
 98:   }

100:   MatPtAPMultEqual(A, P, C, 10, &flg);
102:   MatDestroy(&C);

104:   /* (2) User API */
105:   MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C);
106:   /* Test MAT_REUSE_MATRIX - reuse symbolic C */
107:   alpha = 1.0;
108:   for (i = 0; i < 1; i++) {
109:     alpha -= 0.1;
110:     MatScale(A, alpha);
111:     MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C);
112:   }

114:   /* Free intermediate data structures created for reuse of C=Pt*A*P */
115:   MatProductClear(C);

117:   MatPtAPMultEqual(A, P, C, 10, &flg);

120:   MatDestroy(&C);
121:   MatDestroy(&A);
122:   MatDestroy(&P);
123:   DMDestroy(&finedm);
124:   DMDestroy(&coarsedm);
125:   PetscFinalize();
126:   return 0;
127: }

129: /*TEST

131:    test:
132:       args: -M 10 -N 10 -Z 10
133:       output_file: output/ex89_1.out

135:    test:
136:       suffix: allatonce
137:       nsize: 4
138:       args: -M 10 -N 10 -Z 10
139:       output_file: output/ex89_2.out

141:    test:
142:       suffix: allatonce_merged
143:       nsize: 4
144:       args: -M 10 -M 5 -M 10 -mat_product_algorithm allatonce_merged
145:       output_file: output/ex89_3.out

147:    test:
148:       suffix: nonscalable_3D
149:       nsize: 4
150:       args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm nonscalable
151:       output_file: output/ex89_4.out

153:    test:
154:       suffix: allatonce_merged_3D
155:       nsize: 4
156:       args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm allatonce_merged
157:       output_file: output/ex89_3.out

159:    test:
160:       suffix: nonscalable
161:       nsize: 4
162:       args: -M 10 -N 10 -Z 10 -mat_product_algorithm nonscalable
163:       output_file: output/ex89_5.out

165: TEST*/