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;

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

 29:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_3D", &Test_3D, NULL));
 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 31:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Z", &Z, NULL));
 33:   /* Set up distributed array for fine grid */
 34:   if (!Test_3D) {
 35:     PetscCall(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:     PetscCall(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:   PetscCall(DMSetFromOptions(coarsedm));
 40:   PetscCall(DMSetUp(coarsedm));

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

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

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

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

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

100:   PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
101:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
102:   PetscCall(MatDestroy(&C));

104:   /* (2) User API */
105:   PetscCall(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:     PetscCall(MatScale(A, alpha));
111:     PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C));
112:   }

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

117:   PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
118:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP");

120:   PetscCall(MatDestroy(&C));
121:   PetscCall(MatDestroy(&A));
122:   PetscCall(MatDestroy(&P));
123:   PetscCall(DMDestroy(&finedm));
124:   PetscCall(DMDestroy(&coarsedm));
125:   PetscCall(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*/