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