1 2 /* 3 Laplacian in 3D. Use for testing MatSolve routines. 4 Modeled by the partial differential equation 5 6 - Laplacian u = 1,0 < x,y,z < 1, 7 8 with boundary conditions 9 u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. 10 */ 11 12 static char help[] = "This example is for testing different MatSolve routines :MatSolve(), MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), and MatMatSolve().\n\ 13 Example usage: ./ex129 -mat_type aij -dof 2\n\n"; 14 15 #include <petscdm.h> 16 #include <petscdmda.h> 17 18 extern PetscErrorCode ComputeMatrix(DM, Mat); 19 extern PetscErrorCode ComputeRHS(DM, Vec); 20 extern PetscErrorCode ComputeRHSMatrix(PetscInt, PetscInt, Mat *); 21 22 int main(int argc, char **args) { 23 PetscMPIInt size; 24 Vec x, b, y, b1; 25 DM da; 26 Mat A, F, RHS, X, C1; 27 MatFactorInfo info; 28 IS perm, iperm; 29 PetscInt dof = 1, M = 8, m, n, nrhs; 30 PetscScalar one = 1.0; 31 PetscReal norm, tol = 1000 * PETSC_MACHINE_EPSILON; 32 PetscBool InplaceLU = PETSC_FALSE; 33 34 PetscFunctionBeginUser; 35 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 36 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 37 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only"); 38 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL)); 39 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 40 41 PetscCall(DMDACreate(PETSC_COMM_WORLD, &da)); 42 PetscCall(DMSetDimension(da, 3)); 43 PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE)); 44 PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR)); 45 PetscCall(DMDASetSizes(da, M, M, M)); 46 PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE)); 47 PetscCall(DMDASetDof(da, dof)); 48 PetscCall(DMDASetStencilWidth(da, 1)); 49 PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL)); 50 PetscCall(DMSetMatType(da, MATBAIJ)); 51 PetscCall(DMSetFromOptions(da)); 52 PetscCall(DMSetUp(da)); 53 54 PetscCall(DMCreateGlobalVector(da, &x)); 55 PetscCall(DMCreateGlobalVector(da, &b)); 56 PetscCall(VecDuplicate(b, &y)); 57 PetscCall(ComputeRHS(da, b)); 58 PetscCall(VecSet(y, one)); 59 PetscCall(DMCreateMatrix(da, &A)); 60 PetscCall(ComputeMatrix(da, A)); 61 PetscCall(MatGetSize(A, &m, &n)); 62 nrhs = 2; 63 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 64 PetscCall(ComputeRHSMatrix(m, nrhs, &RHS)); 65 PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &X)); 66 67 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); 68 69 PetscCall(PetscOptionsGetBool(NULL, NULL, "-inplacelu", &InplaceLU, NULL)); 70 PetscCall(MatFactorInfoInitialize(&info)); 71 if (!InplaceLU) { 72 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F)); 73 info.fill = 5.0; 74 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info)); 75 PetscCall(MatLUFactorNumeric(F, A, &info)); 76 } else { /* Test inplace factorization */ 77 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F)); 78 PetscCall(MatLUFactor(F, perm, iperm, &info)); 79 } 80 81 PetscCall(VecDuplicate(y, &b1)); 82 83 /* MatSolve */ 84 PetscCall(MatSolve(F, b, x)); 85 PetscCall(MatMult(A, x, b1)); 86 PetscCall(VecAXPY(b1, -1.0, b)); 87 PetscCall(VecNorm(b1, NORM_2, &norm)); 88 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve : Error of norm %g\n", (double)norm)); } 89 90 /* MatSolveTranspose */ 91 PetscCall(MatSolveTranspose(F, b, x)); 92 PetscCall(MatMultTranspose(A, x, b1)); 93 PetscCall(VecAXPY(b1, -1.0, b)); 94 PetscCall(VecNorm(b1, NORM_2, &norm)); 95 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose : Error of norm %g\n", (double)norm)); } 96 97 /* MatSolveAdd */ 98 PetscCall(MatSolveAdd(F, b, y, x)); 99 PetscCall(MatMult(A, y, b1)); 100 PetscCall(VecScale(b1, -1.0)); 101 PetscCall(MatMultAdd(A, x, b1, b1)); 102 PetscCall(VecAXPY(b1, -1.0, b)); 103 PetscCall(VecNorm(b1, NORM_2, &norm)); 104 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveAdd : Error of norm %g\n", (double)norm)); } 105 106 /* MatSolveTransposeAdd */ 107 PetscCall(MatSolveTransposeAdd(F, b, y, x)); 108 PetscCall(MatMultTranspose(A, y, b1)); 109 PetscCall(VecScale(b1, -1.0)); 110 PetscCall(MatMultTransposeAdd(A, x, b1, b1)); 111 PetscCall(VecAXPY(b1, -1.0, b)); 112 PetscCall(VecNorm(b1, NORM_2, &norm)); 113 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTransposeAdd : Error of norm %g\n", (double)norm)); } 114 115 /* MatMatSolve */ 116 PetscCall(MatMatSolve(F, RHS, X)); 117 PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, 2.0, &C1)); 118 PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN)); 119 PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm)); 120 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve : Error of norm %g\n", (double)norm)); } 121 122 PetscCall(VecDestroy(&x)); 123 PetscCall(VecDestroy(&b)); 124 PetscCall(VecDestroy(&b1)); 125 PetscCall(VecDestroy(&y)); 126 PetscCall(MatDestroy(&A)); 127 PetscCall(MatDestroy(&F)); 128 PetscCall(MatDestroy(&RHS)); 129 PetscCall(MatDestroy(&C1)); 130 PetscCall(MatDestroy(&X)); 131 PetscCall(ISDestroy(&perm)); 132 PetscCall(ISDestroy(&iperm)); 133 PetscCall(DMDestroy(&da)); 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 PetscErrorCode ComputeRHS(DM da, Vec b) { 139 PetscInt mx, my, mz; 140 PetscScalar h; 141 142 PetscFunctionBegin; 143 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 144 h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1)); 145 PetscCall(VecSet(b, h)); 146 PetscFunctionReturn(0); 147 } 148 149 PetscErrorCode ComputeRHSMatrix(PetscInt m, PetscInt nrhs, Mat *C) { 150 PetscRandom rand; 151 Mat RHS; 152 PetscScalar *array, rval; 153 PetscInt i, k; 154 155 PetscFunctionBegin; 156 PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS)); 157 PetscCall(MatSetSizes(RHS, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 158 PetscCall(MatSetType(RHS, MATSEQDENSE)); 159 PetscCall(MatSetUp(RHS)); 160 161 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 162 PetscCall(PetscRandomSetFromOptions(rand)); 163 PetscCall(MatDenseGetArray(RHS, &array)); 164 for (i = 0; i < m; i++) { 165 PetscCall(PetscRandomGetValue(rand, &rval)); 166 array[i] = rval; 167 } 168 if (nrhs > 1) { 169 for (k = 1; k < nrhs; k++) { 170 for (i = 0; i < m; i++) { array[m * k + i] = array[i]; } 171 } 172 } 173 PetscCall(MatDenseRestoreArray(RHS, &array)); 174 PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY)); 175 PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY)); 176 *C = RHS; 177 PetscCall(PetscRandomDestroy(&rand)); 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode ComputeMatrix(DM da, Mat B) { 182 PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3; 183 PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy, r1, r2; 184 MatStencil row, col; 185 PetscRandom rand; 186 187 PetscFunctionBegin; 188 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 189 PetscCall(PetscRandomSetSeed(rand, 1)); 190 PetscCall(PetscRandomSetInterval(rand, -.001, .001)); 191 PetscCall(PetscRandomSetFromOptions(rand)); 192 193 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 194 /* For simplicity, this example only works on mx=my=mz */ 195 PetscCheck(mx == my && mx == mz, PETSC_COMM_SELF, PETSC_ERR_SUP, "This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT, mx, my, mz); 196 197 Hx = 1.0 / (PetscReal)(mx - 1); 198 Hy = 1.0 / (PetscReal)(my - 1); 199 Hz = 1.0 / (PetscReal)(mz - 1); 200 HxHydHz = Hx * Hy / Hz; 201 HxHzdHy = Hx * Hz / Hy; 202 HyHzdHx = Hy * Hz / Hx; 203 204 PetscCall(PetscMalloc1(2 * dof * dof + 1, &v)); 205 v_neighbor = v + dof * dof; 206 PetscCall(PetscArrayzero(v, 2 * dof * dof + 1)); 207 k3 = 0; 208 for (k1 = 0; k1 < dof; k1++) { 209 for (k2 = 0; k2 < dof; k2++) { 210 if (k1 == k2) { 211 v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx); 212 v_neighbor[k3] = -HxHydHz; 213 } else { 214 PetscCall(PetscRandomGetValue(rand, &r1)); 215 PetscCall(PetscRandomGetValue(rand, &r2)); 216 217 v[k3] = r1; 218 v_neighbor[k3] = r2; 219 } 220 k3++; 221 } 222 } 223 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm)); 224 225 for (k = zs; k < zs + zm; k++) { 226 for (j = ys; j < ys + ym; j++) { 227 for (i = xs; i < xs + xm; i++) { 228 row.i = i; 229 row.j = j; 230 row.k = k; 231 if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */ 232 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES)); 233 } else { /* interior points */ 234 /* center */ 235 col.i = i; 236 col.j = j; 237 col.k = k; 238 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES)); 239 240 /* x neighbors */ 241 col.i = i - 1; 242 col.j = j; 243 col.k = k; 244 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 245 col.i = i + 1; 246 col.j = j; 247 col.k = k; 248 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 249 250 /* y neighbors */ 251 col.i = i; 252 col.j = j - 1; 253 col.k = k; 254 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 255 col.i = i; 256 col.j = j + 1; 257 col.k = k; 258 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 259 260 /* z neighbors */ 261 col.i = i; 262 col.j = j; 263 col.k = k - 1; 264 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 265 col.i = i; 266 col.j = j; 267 col.k = k + 1; 268 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 269 } 270 } 271 } 272 } 273 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 274 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 275 PetscCall(PetscFree(v)); 276 PetscCall(PetscRandomDestroy(&rand)); 277 PetscFunctionReturn(0); 278 } 279 280 /*TEST 281 282 test: 283 args: -dm_mat_type aij -dof 1 284 output_file: output/ex129.out 285 286 test: 287 suffix: 2 288 args: -dm_mat_type aij -dof 1 -inplacelu 289 output_file: output/ex129.out 290 291 TEST*/ 292