1 /* 2 Laplacian in 3D. Use for testing MatSolve routines. 3 Modeled by the partial differential equation 4 5 - Laplacian u = 1,0 < x,y,z < 1, 6 7 with boundary conditions 8 u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. 9 */ 10 11 static char help[] = "This example is for testing different MatSolve routines :MatSolve(), MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), and MatMatSolve().\n\ 12 Example usage: ./ex129 -mat_type aij -dof 2\n\n"; 13 14 #include <petscdm.h> 15 #include <petscdmda.h> 16 17 extern PetscErrorCode ComputeMatrix(DM, Mat); 18 extern PetscErrorCode ComputeRHS(DM, Vec); 19 extern PetscErrorCode ComputeRHSMatrix(PetscInt, PetscInt, Mat *); 20 21 int main(int argc, char **args) 22 { 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 { 140 PetscInt mx, my, mz; 141 PetscScalar h; 142 143 PetscFunctionBegin; 144 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 145 h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1)); 146 PetscCall(VecSet(b, h)); 147 PetscFunctionReturn(PETSC_SUCCESS); 148 } 149 150 PetscErrorCode ComputeRHSMatrix(PetscInt m, PetscInt nrhs, Mat *C) 151 { 152 PetscRandom rand; 153 Mat RHS; 154 PetscScalar *array, rval; 155 PetscInt i, k; 156 157 PetscFunctionBegin; 158 PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS)); 159 PetscCall(MatSetSizes(RHS, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 160 PetscCall(MatSetType(RHS, MATSEQDENSE)); 161 PetscCall(MatSetUp(RHS)); 162 163 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 164 PetscCall(PetscRandomSetFromOptions(rand)); 165 PetscCall(MatDenseGetArray(RHS, &array)); 166 for (i = 0; i < m; i++) { 167 PetscCall(PetscRandomGetValue(rand, &rval)); 168 array[i] = rval; 169 } 170 if (nrhs > 1) { 171 for (k = 1; k < nrhs; k++) { 172 for (i = 0; i < m; i++) array[m * k + i] = array[i]; 173 } 174 } 175 PetscCall(MatDenseRestoreArray(RHS, &array)); 176 PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY)); 177 PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY)); 178 *C = RHS; 179 PetscCall(PetscRandomDestroy(&rand)); 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 PetscErrorCode ComputeMatrix(DM da, Mat B) 184 { 185 PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3; 186 PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy, r1, r2; 187 MatStencil row, col; 188 PetscRandom rand; 189 190 PetscFunctionBegin; 191 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 192 PetscCall(PetscRandomSetSeed(rand, 1)); 193 PetscCall(PetscRandomSetInterval(rand, -.001, .001)); 194 PetscCall(PetscRandomSetFromOptions(rand)); 195 196 PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 197 /* For simplicity, this example only works on mx=my=mz */ 198 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); 199 200 Hx = 1.0 / (PetscReal)(mx - 1); 201 Hy = 1.0 / (PetscReal)(my - 1); 202 Hz = 1.0 / (PetscReal)(mz - 1); 203 HxHydHz = Hx * Hy / Hz; 204 HxHzdHy = Hx * Hz / Hy; 205 HyHzdHx = Hy * Hz / Hx; 206 207 PetscCall(PetscMalloc1(2 * dof * dof + 1, &v)); 208 v_neighbor = v + dof * dof; 209 PetscCall(PetscArrayzero(v, 2 * dof * dof + 1)); 210 k3 = 0; 211 for (k1 = 0; k1 < dof; k1++) { 212 for (k2 = 0; k2 < dof; k2++) { 213 if (k1 == k2) { 214 v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx); 215 v_neighbor[k3] = -HxHydHz; 216 } else { 217 PetscCall(PetscRandomGetValue(rand, &r1)); 218 PetscCall(PetscRandomGetValue(rand, &r2)); 219 220 v[k3] = r1; 221 v_neighbor[k3] = r2; 222 } 223 k3++; 224 } 225 } 226 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm)); 227 228 for (k = zs; k < zs + zm; k++) { 229 for (j = ys; j < ys + ym; j++) { 230 for (i = xs; i < xs + xm; i++) { 231 row.i = i; 232 row.j = j; 233 row.k = k; 234 if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */ 235 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES)); 236 } else { /* interior points */ 237 /* center */ 238 col.i = i; 239 col.j = j; 240 col.k = k; 241 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES)); 242 243 /* x neighbors */ 244 col.i = i - 1; 245 col.j = j; 246 col.k = k; 247 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 248 col.i = i + 1; 249 col.j = j; 250 col.k = k; 251 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 252 253 /* y neighbors */ 254 col.i = i; 255 col.j = j - 1; 256 col.k = k; 257 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 258 col.i = i; 259 col.j = j + 1; 260 col.k = k; 261 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 262 263 /* z neighbors */ 264 col.i = i; 265 col.j = j; 266 col.k = k - 1; 267 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 268 col.i = i; 269 col.j = j; 270 col.k = k + 1; 271 PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES)); 272 } 273 } 274 } 275 } 276 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 277 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 278 PetscCall(PetscFree(v)); 279 PetscCall(PetscRandomDestroy(&rand)); 280 PetscFunctionReturn(PETSC_SUCCESS); 281 } 282 283 /*TEST 284 285 test: 286 args: -dm_mat_type aij -dof 1 287 output_file: output/ex129.out 288 289 test: 290 suffix: 2 291 args: -dm_mat_type aij -dof 1 -inplacelu 292 output_file: output/ex129.out 293 294 TEST*/ 295