1 static char help[] = "Tests sequential and parallel MatMatMatMult() and MatPtAP(). Modified from ex96.c \n\ 2 -Mx <xg>, where <xg> = number of coarse grid points in the x-direction\n\ 3 -My <yg>, where <yg> = number of coarse grid points in the y-direction\n\ 4 -Mz <zg>, where <zg> = number of coarse grid points in the z-direction\n\ 5 -Npx <npx>, where <npx> = number of processors in the x-direction\n\ 6 -Npy <npy>, where <npy> = number of processors in the y-direction\n\ 7 -Npz <npz>, where <npz> = number of processors in the z-direction\n\n"; 8 9 /* 10 Example of usage: mpiexec -n 3 ./ex41 -Mx 10 -My 10 -Mz 10 11 */ 12 13 #include <petscdm.h> 14 #include <petscdmda.h> 15 16 /* User-defined application contexts */ 17 typedef struct { 18 PetscInt mx, my, mz; /* number grid points in x, y and z direction */ 19 Vec localX, localF; /* local vectors with ghost region */ 20 DM da; 21 Vec x, b, r; /* global vectors */ 22 Mat J; /* Jacobian on grid */ 23 } GridCtx; 24 typedef struct { 25 GridCtx fine; 26 GridCtx coarse; 27 PetscInt ratio; 28 Mat Ii; /* interpolation from coarse to fine */ 29 } AppCtx; 30 31 #define COARSE_LEVEL 0 32 #define FINE_LEVEL 1 33 34 /* 35 Mm_ratio - ration of grid lines between fine and coarse grids. 36 */ 37 int main(int argc, char **argv) 38 { 39 AppCtx user; 40 PetscMPIInt size, rank; 41 PetscInt m, n, M, N, i, nrows; 42 PetscScalar one = 1.0; 43 PetscReal fill = 2.0; 44 Mat A, P, R, C, PtAP, D; 45 PetscScalar *array; 46 PetscRandom rdm; 47 PetscBool Test_3D = PETSC_FALSE, flg; 48 const PetscInt *ia, *ja; 49 50 PetscFunctionBeginUser; 51 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 52 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 53 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 54 55 /* Get size of fine grids and coarse grids */ 56 user.ratio = 2; 57 user.coarse.mx = 4; 58 user.coarse.my = 4; 59 user.coarse.mz = 4; 60 61 PetscCall(PetscOptionsGetInt(NULL, NULL, "-Mx", &user.coarse.mx, NULL)); 62 PetscCall(PetscOptionsGetInt(NULL, NULL, "-My", &user.coarse.my, NULL)); 63 PetscCall(PetscOptionsGetInt(NULL, NULL, "-Mz", &user.coarse.mz, NULL)); 64 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ratio", &user.ratio, NULL)); 65 if (user.coarse.mz) Test_3D = PETSC_TRUE; 66 67 user.fine.mx = user.ratio * (user.coarse.mx - 1) + 1; 68 user.fine.my = user.ratio * (user.coarse.my - 1) + 1; 69 user.fine.mz = user.ratio * (user.coarse.mz - 1) + 1; 70 71 if (rank == 0) { 72 if (!Test_3D) { 73 PetscCall(PetscPrintf(PETSC_COMM_SELF, "coarse grids: %" PetscInt_FMT " %" PetscInt_FMT "; fine grids: %" PetscInt_FMT " %" PetscInt_FMT "\n", user.coarse.mx, user.coarse.my, user.fine.mx, user.fine.my)); 74 } else { 75 PetscCall(PetscPrintf(PETSC_COMM_SELF, "coarse grids: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "; fine grids: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", user.coarse.mx, user.coarse.my, user.coarse.mz, user.fine.mx, 76 user.fine.my, user.fine.mz)); 77 } 78 } 79 80 /* Set up distributed array for fine grid */ 81 if (!Test_3D) { 82 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.fine.mx, user.fine.my, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.fine.da)); 83 } else { 84 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.fine.mx, user.fine.my, user.fine.mz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &user.fine.da)); 85 } 86 PetscCall(DMSetFromOptions(user.fine.da)); 87 PetscCall(DMSetUp(user.fine.da)); 88 89 /* Create and set A at fine grids */ 90 PetscCall(DMSetMatType(user.fine.da, MATAIJ)); 91 PetscCall(DMCreateMatrix(user.fine.da, &A)); 92 PetscCall(MatGetLocalSize(A, &m, &n)); 93 PetscCall(MatGetSize(A, &M, &N)); 94 95 /* set val=one to A (replace with random values!) */ 96 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 97 PetscCall(PetscRandomSetFromOptions(rdm)); 98 if (size == 1) { 99 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 100 if (flg) { 101 PetscCall(MatSeqAIJGetArray(A, &array)); 102 for (i = 0; i < ia[nrows]; i++) array[i] = one; 103 PetscCall(MatSeqAIJRestoreArray(A, &array)); 104 } 105 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 106 } else { 107 Mat AA, AB; 108 PetscCall(MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL)); 109 PetscCall(MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 110 if (flg) { 111 PetscCall(MatSeqAIJGetArray(AA, &array)); 112 for (i = 0; i < ia[nrows]; i++) array[i] = one; 113 PetscCall(MatSeqAIJRestoreArray(AA, &array)); 114 } 115 PetscCall(MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 116 PetscCall(MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 117 if (flg) { 118 PetscCall(MatSeqAIJGetArray(AB, &array)); 119 for (i = 0; i < ia[nrows]; i++) array[i] = one; 120 PetscCall(MatSeqAIJRestoreArray(AB, &array)); 121 } 122 PetscCall(MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg)); 123 } 124 /* Set up distributed array for coarse grid */ 125 if (!Test_3D) { 126 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.coarse.mx, user.coarse.my, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.coarse.da)); 127 } else { 128 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.coarse.mx, user.coarse.my, user.coarse.mz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &user.coarse.da)); 129 } 130 PetscCall(DMSetFromOptions(user.coarse.da)); 131 PetscCall(DMSetUp(user.coarse.da)); 132 133 /* Create interpolation between the fine and coarse grids */ 134 PetscCall(DMCreateInterpolation(user.coarse.da, user.fine.da, &P, NULL)); 135 136 /* Get R = P^T */ 137 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R)); 138 139 /* C = R*A*P */ 140 /* Developer's API */ 141 PetscCall(MatProductCreate(R, A, P, &D)); 142 PetscCall(MatProductSetType(D, MATPRODUCT_ABC)); 143 PetscCall(MatProductSetFromOptions(D)); 144 PetscCall(MatProductSymbolic(D)); 145 PetscCall(MatProductNumeric(D)); 146 PetscCall(MatProductNumeric(D)); /* Test reuse symbolic D */ 147 148 /* User's API */ 149 { /* Test MatMatMatMult_Basic() */ 150 Mat Adense, Cdense; 151 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense)); 152 PetscCall(MatMatMatMult(R, Adense, P, MAT_INITIAL_MATRIX, fill, &Cdense)); 153 PetscCall(MatMatMatMult(R, Adense, P, MAT_REUSE_MATRIX, fill, &Cdense)); 154 155 PetscCall(MatMultEqual(D, Cdense, 10, &flg)); 156 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D*v != Cdense*v"); 157 PetscCall(MatDestroy(&Adense)); 158 PetscCall(MatDestroy(&Cdense)); 159 } 160 161 PetscCall(MatMatMatMult(R, A, P, MAT_INITIAL_MATRIX, fill, &C)); 162 PetscCall(MatMatMatMult(R, A, P, MAT_REUSE_MATRIX, fill, &C)); 163 PetscCall(MatProductClear(C)); 164 165 /* Test D == C */ 166 PetscCall(MatEqual(D, C, &flg)); 167 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "D != C"); 168 169 /* Test C == PtAP */ 170 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &PtAP)); 171 PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &PtAP)); 172 PetscCall(MatEqual(C, PtAP, &flg)); 173 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "C != PtAP"); 174 PetscCall(MatDestroy(&PtAP)); 175 176 /* Clean up */ 177 PetscCall(MatDestroy(&A)); 178 PetscCall(PetscRandomDestroy(&rdm)); 179 PetscCall(DMDestroy(&user.fine.da)); 180 PetscCall(DMDestroy(&user.coarse.da)); 181 PetscCall(MatDestroy(&P)); 182 PetscCall(MatDestroy(&R)); 183 PetscCall(MatDestroy(&C)); 184 PetscCall(MatDestroy(&D)); 185 PetscCall(PetscFinalize()); 186 return 0; 187 } 188 189 /*TEST 190 191 test: 192 193 test: 194 suffix: 2 195 nsize: 2 196 args: -matmatmatmult_via scalable 197 198 test: 199 suffix: 3 200 nsize: 2 201 args: -matmatmatmult_via nonscalable 202 output_file: output/ex111_1.out 203 204 TEST*/ 205