1 static char help[] = "Tests MatForwardSolve and MatBackwardSolve for LU and Cholesky decompositions using C/Pardiso.\n\n"; 2 3 #include <petscdmda.h> 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 DM da; 9 DMDALocalInfo info; 10 Mat A, F; 11 Vec x, y, b, ytmp; 12 IS rowis, colis; 13 PetscInt i, j, k, n = 5; 14 PetscBool CHOL = PETSC_FALSE; 15 MatStencil row, cols[5]; 16 PetscScalar vals[5]; 17 PetscReal norm2, tol = 100. * PETSC_MACHINE_EPSILON; 18 PetscRandom rdm; 19 MatFactorInfo finfo; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 24 25 PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, n, n, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da)); 26 PetscCall(DMSetUp(da)); 27 28 PetscCall(DMCreateMatrix(da, &A)); 29 PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOL)); 30 if (CHOL) PetscCall(MatSetType(A, MATSBAIJ)); 31 else PetscCall(MatSetType(A, MATBAIJ)); 32 PetscCall(MatSetUp(A)); 33 34 PetscCall(DMDAGetLocalInfo(da, &info)); 35 for (j = info.ys; j < info.ys + info.ym; j++) { 36 for (i = info.xs; i < info.xs + info.xm; i++) { 37 row.j = j; 38 row.i = i; 39 40 k = 0; 41 if (j != 0) { 42 cols[k].j = j - 1; 43 cols[k].i = i; 44 vals[k] = -1; 45 ++k; 46 } 47 if (i != 0) { 48 cols[k].j = j; 49 cols[k].i = i - 1; 50 vals[k] = -1; 51 ++k; 52 } 53 cols[k].j = j; 54 cols[k].i = i; 55 vals[k] = 4; 56 ++k; 57 if (j != info.my - 1) { 58 cols[k].j = j + 1; 59 cols[k].i = i; 60 vals[k] = -1; 61 ++k; 62 } 63 if (i != info.mx - 1) { 64 cols[k].j = j; 65 cols[k].i = i + 1; 66 vals[k] = -1; 67 ++k; 68 } 69 PetscCall(MatSetValuesStencil(A, 1, &row, k, cols, vals, INSERT_VALUES)); 70 } 71 } 72 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 73 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 74 75 if (CHOL) PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE)); 76 77 /* Create vectors for error checking */ 78 PetscCall(MatCreateVecs(A, &x, &b)); 79 PetscCall(VecDuplicate(x, &y)); 80 PetscCall(VecDuplicate(x, &ytmp)); 81 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 82 PetscCall(PetscRandomSetFromOptions(rdm)); 83 PetscCall(VecSetRandom(x, rdm)); 84 PetscCall(MatMult(A, x, b)); 85 86 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &rowis, &colis)); 87 if (CHOL) { 88 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test Cholesky...\n")); 89 PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_CHOLESKY, &F)); 90 PetscCall(MatCholeskyFactorSymbolic(F, A, rowis, &finfo)); 91 PetscCall(MatCholeskyFactorNumeric(F, A, &finfo)); 92 } else { 93 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test LU...\n")); 94 PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_LU, &F)); 95 PetscCall(MatLUFactorSymbolic(F, A, rowis, colis, &finfo)); 96 PetscCall(MatLUFactorNumeric(F, A, &finfo)); 97 } 98 99 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatForwardSolve...\n")); 100 PetscCall(MatForwardSolve(F, b, ytmp)); 101 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatBackwardSolve...\n")); 102 PetscCall(MatBackwardSolve(F, ytmp, y)); 103 PetscCall(VecAXPY(y, -1.0, x)); 104 PetscCall(VecNorm(y, NORM_2, &norm2)); 105 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2)); 106 107 /* Free data structures */ 108 PetscCall(ISDestroy(&rowis)); 109 PetscCall(ISDestroy(&colis)); 110 PetscCall(MatDestroy(&F)); 111 PetscCall(MatDestroy(&A)); 112 PetscCall(DMDestroy(&da)); 113 PetscCall(PetscRandomDestroy(&rdm)); 114 PetscCall(VecDestroy(&x)); 115 PetscCall(VecDestroy(&y)); 116 PetscCall(VecDestroy(&ytmp)); 117 PetscCall(VecDestroy(&b)); 118 PetscCall(PetscFinalize()); 119 return 0; 120 } 121 122 /*TEST 123 124 test: 125 requires: mkl_cpardiso 126 nsize: 4 127 128 test: 129 suffix: 2 130 requires: !complex mkl_cpardiso 131 nsize: 4 132 args: -chol 133 134 TEST*/ 135