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