1 static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqaij format, and illustrates drawing of matrix sparsity structure with MatView().\n\ 2 Input parameters are:\n\ 3 -lf <level> : level of fill for ILU (default is 0)\n\ 4 -lu : use full LU or Cholesky factorization\n\ 5 -m <value>,-n <value> : grid dimensions\n\ 6 Note that most users should employ the KSP interface to the\n\ 7 linear solvers instead of using the factorization routines\n\ 8 directly.\n\n"; 9 10 #include <petscmat.h> 11 12 int main(int argc, char **args) 13 { 14 Mat C, A; 15 PetscInt i, j, m = 5, n = 5, Ii, J, lf = 0; 16 PetscBool LU = PETSC_FALSE, CHOLESKY = PETSC_FALSE, TRIANGULAR = PETSC_FALSE, MATDSPL = PETSC_FALSE, flg, matordering, use_mkl_pardiso = PETSC_FALSE; 17 PetscScalar v; 18 IS row, col; 19 PetscViewer viewer1, viewer2; 20 MatFactorInfo info; 21 Vec x, y, b, ytmp; 22 PetscReal norm2, norm2_inplace, tol = 100. * PETSC_MACHINE_EPSILON; 23 PetscRandom rdm; 24 PetscMPIInt size; 25 char pack[PETSC_MAX_PATH_LEN]; 26 27 PetscFunctionBeginUser; 28 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 29 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 30 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL)); 34 35 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 0, 0, 400, 400, &viewer1)); 36 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 400, 0, 400, 400, &viewer2)); 37 38 PetscCall(MatCreate(PETSC_COMM_SELF, &C)); 39 PetscCall(MatSetSizes(C, m * n, m * n, m * n, m * n)); 40 PetscCall(MatSetFromOptions(C)); 41 PetscCall(MatSetUp(C)); 42 43 /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */ 44 for (i = 0; i < m; i++) { 45 for (j = 0; j < n; j++) { 46 v = -1.0; 47 Ii = j + n * i; 48 J = Ii - n; 49 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 50 J = Ii + n; 51 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 52 J = Ii - 1; 53 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 54 J = Ii + 1; 55 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 56 v = 4.0; 57 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 58 } 59 } 60 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 61 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 62 63 PetscCall(MatIsSymmetric(C, 0.0, &flg)); 64 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric"); 65 66 PetscCall(MatSetOption(C, MAT_SPD, PETSC_TRUE)); 67 68 /* Create vectors for error checking */ 69 PetscCall(MatCreateVecs(C, &x, &b)); 70 PetscCall(VecDuplicate(x, &y)); 71 PetscCall(VecDuplicate(x, &ytmp)); 72 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 73 PetscCall(PetscRandomSetFromOptions(rdm)); 74 PetscCall(VecSetRandom(x, rdm)); 75 PetscCall(MatMult(C, x, b)); 76 77 PetscCall(PetscOptionsHasName(NULL, NULL, "-mat_ordering", &matordering)); 78 if (matordering) { 79 PetscCall(MatGetOrdering(C, MATORDERINGRCM, &row, &col)); 80 } else { 81 PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col)); 82 } 83 84 PetscCall(PetscOptionsHasName(NULL, NULL, "-display_matrices", &MATDSPL)); 85 if (MATDSPL) { 86 printf("original matrix:\n"); 87 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO)); 88 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF)); 89 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF)); 90 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF)); 91 PetscCall(MatView(C, viewer1)); 92 } 93 94 /* Compute LU or ILU factor A */ 95 PetscCall(MatFactorInfoInitialize(&info)); 96 97 info.fill = 1.0; 98 info.diagonal_fill = 0; 99 info.zeropivot = 0.0; 100 101 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL)); 102 #if defined(PETSC_HAVE_MKL_PARDISO) 103 PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &use_mkl_pardiso)); 104 #endif 105 106 PetscCall(PetscOptionsHasName(NULL, NULL, "-lu", &LU)); 107 if (LU) { 108 printf("Test LU...\n"); 109 if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &A)); 110 else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A)); 111 PetscCall(MatLUFactorSymbolic(A, C, row, col, &info)); 112 } else { 113 printf("Test ILU...\n"); 114 info.levels = lf; 115 116 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A)); 117 PetscCall(MatILUFactorSymbolic(A, C, row, col, &info)); 118 } 119 PetscCall(MatLUFactorNumeric(A, C, &info)); 120 121 /* test MatForwardSolve() and MatBackwardSolve() with MKL Pardiso*/ 122 if (LU && use_mkl_pardiso) { 123 PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR)); 124 if (TRIANGULAR) { 125 printf("Test MatForwardSolve...\n"); 126 PetscCall(MatForwardSolve(A, b, ytmp)); 127 printf("Test MatBackwardSolve...\n"); 128 PetscCall(MatBackwardSolve(A, ytmp, y)); 129 PetscCall(VecAXPY(y, -1.0, x)); 130 PetscCall(VecNorm(y, NORM_2, &norm2)); 131 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2)); 132 } 133 } 134 135 /* Solve A*y = b, then check the error */ 136 PetscCall(MatSolve(A, b, y)); 137 PetscCall(VecAXPY(y, -1.0, x)); 138 PetscCall(VecNorm(y, NORM_2, &norm2)); 139 PetscCall(MatDestroy(&A)); 140 141 /* Test in-place ILU(0) and compare it with the out-place ILU(0) */ 142 if (!LU && lf == 0) { 143 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &A)); 144 PetscCall(MatILUFactor(A, row, col, &info)); 145 /* 146 printf("In-place factored matrix:\n"); 147 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF)); 148 */ 149 PetscCall(MatSolve(A, b, y)); 150 PetscCall(VecAXPY(y, -1.0, x)); 151 PetscCall(VecNorm(y, NORM_2, &norm2_inplace)); 152 PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ILU(0) %g and in-place ILU(0) %g give different residuals", (double)norm2, (double)norm2_inplace); 153 PetscCall(MatDestroy(&A)); 154 } 155 156 /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */ 157 PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOLESKY)); 158 if (CHOLESKY) { 159 printf("Test Cholesky...\n"); 160 lf = -1; 161 if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &A)); 162 else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &A)); 163 PetscCall(MatCholeskyFactorSymbolic(A, C, row, &info)); 164 } else { 165 printf("Test ICC...\n"); 166 info.levels = lf; 167 info.fill = 1.0; 168 info.diagonal_fill = 0; 169 info.zeropivot = 0.0; 170 171 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A)); 172 PetscCall(MatICCFactorSymbolic(A, C, row, &info)); 173 } 174 PetscCall(MatCholeskyFactorNumeric(A, C, &info)); 175 176 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */ 177 if (CHOLESKY) { 178 PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR)); 179 if (TRIANGULAR) { 180 printf("Test MatForwardSolve...\n"); 181 PetscCall(MatForwardSolve(A, b, ytmp)); 182 printf("Test MatBackwardSolve...\n"); 183 PetscCall(MatBackwardSolve(A, ytmp, y)); 184 PetscCall(VecAXPY(y, -1.0, x)); 185 PetscCall(VecNorm(y, NORM_2, &norm2)); 186 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2)); 187 } 188 } 189 190 PetscCall(MatSolve(A, b, y)); 191 PetscCall(MatDestroy(&A)); 192 PetscCall(VecAXPY(y, -1.0, x)); 193 PetscCall(VecNorm(y, NORM_2, &norm2)); 194 if (lf == -1 && norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2)); 195 196 /* Test in-place ICC(0) and compare it with the out-place ICC(0) */ 197 if (!CHOLESKY && lf == 0 && !matordering) { 198 PetscCall(MatConvert(C, MATSBAIJ, MAT_INITIAL_MATRIX, &A)); 199 PetscCall(MatICCFactor(A, row, &info)); 200 /* 201 printf("In-place factored matrix:\n"); 202 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 203 */ 204 PetscCall(MatSolve(A, b, y)); 205 PetscCall(VecAXPY(y, -1.0, x)); 206 PetscCall(VecNorm(y, NORM_2, &norm2_inplace)); 207 PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ICC(0) %g and in-place ICC(0) %g give different residuals", (double)norm2, (double)norm2_inplace); 208 PetscCall(MatDestroy(&A)); 209 } 210 211 /* Free data structures */ 212 PetscCall(ISDestroy(&row)); 213 PetscCall(ISDestroy(&col)); 214 PetscCall(MatDestroy(&C)); 215 PetscCall(PetscViewerDestroy(&viewer1)); 216 PetscCall(PetscViewerDestroy(&viewer2)); 217 PetscCall(PetscRandomDestroy(&rdm)); 218 PetscCall(VecDestroy(&x)); 219 PetscCall(VecDestroy(&y)); 220 PetscCall(VecDestroy(&ytmp)); 221 PetscCall(VecDestroy(&b)); 222 PetscCall(PetscFinalize()); 223 return 0; 224 } 225 226 /*TEST 227 228 test: 229 args: -mat_ordering -display_matrices -nox 230 filter: grep -v " MPI process" 231 232 test: 233 suffix: 2 234 args: -mat_ordering -display_matrices -nox -lu -chol 235 236 test: 237 suffix: 3 238 args: -mat_ordering -lu -chol -triangular_solve 239 240 test: 241 suffix: 4 242 243 test: 244 suffix: 5 245 args: -lu -chol 246 247 test: 248 suffix: 6 249 args: -lu -chol -triangular_solve 250 output_file: output/ex30_3.out 251 252 test: 253 suffix: 7 254 requires: mkl_pardiso 255 args: -lu -chol -mat_solver_type mkl_pardiso 256 output_file: output/ex30_5.out 257 258 test: 259 suffix: 8 260 requires: mkl_pardiso 261 args: -lu -mat_solver_type mkl_pardiso -triangular_solve 262 263 TEST*/ 264