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