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 { 15 Mat C,A; 16 PetscInt i,j,m = 5,n = 5,Ii,J,lf = 0; 17 PetscBool LU=PETSC_FALSE,CHOLESKY,TRIANGULAR=PETSC_FALSE,MATDSPL=PETSC_FALSE,flg,matordering; 18 PetscScalar v; 19 IS row,col; 20 PetscViewer viewer1,viewer2; 21 MatFactorInfo info; 22 Vec x,y,b,ytmp; 23 PetscReal norm2,norm2_inplace, tol = 100.*PETSC_MACHINE_EPSILON; 24 PetscRandom rdm; 25 PetscMPIInt size; 26 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; Ii = j + n*i; 46 J = Ii - n; if (J>=0) PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 47 J = Ii + n; if (J<m*n) PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 48 J = Ii - 1; if (J>=0) PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 49 J = Ii + 1; if (J<m*n) PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 50 v = 4.0; PetscCall(MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 51 } 52 } 53 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 54 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 55 56 PetscCall(MatIsSymmetric(C,0.0,&flg)); 57 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"C is non-symmetric"); 58 59 /* Create vectors for error checking */ 60 PetscCall(MatCreateVecs(C,&x,&b)); 61 PetscCall(VecDuplicate(x,&y)); 62 PetscCall(VecDuplicate(x,&ytmp)); 63 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 64 PetscCall(PetscRandomSetFromOptions(rdm)); 65 PetscCall(VecSetRandom(x,rdm)); 66 PetscCall(MatMult(C,x,b)); 67 68 PetscCall(PetscOptionsHasName(NULL,NULL,"-mat_ordering",&matordering)); 69 if (matordering) { 70 PetscCall(MatGetOrdering(C,MATORDERINGRCM,&row,&col)); 71 } else { 72 PetscCall(MatGetOrdering(C,MATORDERINGNATURAL,&row,&col)); 73 } 74 75 PetscCall(PetscOptionsHasName(NULL,NULL,"-display_matrices",&MATDSPL)); 76 if (MATDSPL) { 77 printf("original matrix:\n"); 78 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO)); 79 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF)); 80 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF)); 81 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF)); 82 PetscCall(MatView(C,viewer1)); 83 } 84 85 /* Compute LU or ILU factor A */ 86 PetscCall(MatFactorInfoInitialize(&info)); 87 88 info.fill = 1.0; 89 info.diagonal_fill = 0; 90 info.zeropivot = 0.0; 91 92 PetscCall(PetscOptionsHasName(NULL,NULL,"-lu",&LU)); 93 if (LU) { 94 printf("Test LU...\n"); 95 PetscCall(MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A)); 96 PetscCall(MatLUFactorSymbolic(A,C,row,col,&info)); 97 } else { 98 printf("Test ILU...\n"); 99 info.levels = lf; 100 101 PetscCall(MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ILU,&A)); 102 PetscCall(MatILUFactorSymbolic(A,C,row,col,&info)); 103 } 104 PetscCall(MatLUFactorNumeric(A,C,&info)); 105 106 /* Solve A*y = b, then check the error */ 107 PetscCall(MatSolve(A,b,y)); 108 PetscCall(VecAXPY(y,-1.0,x)); 109 PetscCall(VecNorm(y,NORM_2,&norm2)); 110 PetscCall(MatDestroy(&A)); 111 112 /* Test in-place ILU(0) and compare it with the out-place ILU(0) */ 113 if (!LU && lf==0) { 114 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&A)); 115 PetscCall(MatILUFactor(A,row,col,&info)); 116 /* 117 printf("In-place factored matrix:\n"); 118 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF)); 119 */ 120 PetscCall(MatSolve(A,b,y)); 121 PetscCall(VecAXPY(y,-1.0,x)); 122 PetscCall(VecNorm(y,NORM_2,&norm2_inplace)); 123 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); 124 PetscCall(MatDestroy(&A)); 125 } 126 127 /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */ 128 CHOLESKY = LU; 129 if (CHOLESKY) { 130 printf("Test Cholesky...\n"); 131 lf = -1; 132 PetscCall(MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A)); 133 PetscCall(MatCholeskyFactorSymbolic(A,C,row,&info)); 134 } else { 135 printf("Test ICC...\n"); 136 info.levels = lf; 137 info.fill = 1.0; 138 info.diagonal_fill = 0; 139 info.zeropivot = 0.0; 140 141 PetscCall(MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A)); 142 PetscCall(MatICCFactorSymbolic(A,C,row,&info)); 143 } 144 PetscCall(MatCholeskyFactorNumeric(A,C,&info)); 145 146 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */ 147 if (lf == -1) { 148 PetscCall(PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR)); 149 if (TRIANGULAR) { 150 printf("Test MatForwardSolve...\n"); 151 PetscCall(MatForwardSolve(A,b,ytmp)); 152 printf("Test MatBackwardSolve...\n"); 153 PetscCall(MatBackwardSolve(A,ytmp,y)); 154 PetscCall(VecAXPY(y,-1.0,x)); 155 PetscCall(VecNorm(y,NORM_2,&norm2)); 156 if (norm2 > tol) { 157 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2)); 158 } 159 } 160 } 161 162 PetscCall(MatSolve(A,b,y)); 163 PetscCall(MatDestroy(&A)); 164 PetscCall(VecAXPY(y,-1.0,x)); 165 PetscCall(VecNorm(y,NORM_2,&norm2)); 166 if (lf == -1 && norm2 > tol) { 167 PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n",lf,(double)norm2)); 168 } 169 170 /* Test in-place ICC(0) and compare it with the out-place ICC(0) */ 171 if (!CHOLESKY && lf==0 && !matordering) { 172 PetscCall(MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A)); 173 PetscCall(MatICCFactor(A,row,&info)); 174 /* 175 printf("In-place factored matrix:\n"); 176 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 177 */ 178 PetscCall(MatSolve(A,b,y)); 179 PetscCall(VecAXPY(y,-1.0,x)); 180 PetscCall(VecNorm(y,NORM_2,&norm2_inplace)); 181 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); 182 PetscCall(MatDestroy(&A)); 183 } 184 185 /* Free data structures */ 186 PetscCall(ISDestroy(&row)); 187 PetscCall(ISDestroy(&col)); 188 PetscCall(MatDestroy(&C)); 189 PetscCall(PetscViewerDestroy(&viewer1)); 190 PetscCall(PetscViewerDestroy(&viewer2)); 191 PetscCall(PetscRandomDestroy(&rdm)); 192 PetscCall(VecDestroy(&x)); 193 PetscCall(VecDestroy(&y)); 194 PetscCall(VecDestroy(&ytmp)); 195 PetscCall(VecDestroy(&b)); 196 PetscCall(PetscFinalize()); 197 return 0; 198 } 199 200 /*TEST 201 202 test: 203 args: -mat_ordering -display_matrices -nox 204 filter: grep -v " MPI process" 205 206 test: 207 suffix: 2 208 args: -mat_ordering -display_matrices -nox -lu 209 210 test: 211 suffix: 3 212 args: -mat_ordering -lu -triangular_solve 213 214 test: 215 suffix: 4 216 217 test: 218 suffix: 5 219 args: -lu 220 221 test: 222 suffix: 6 223 args: -lu -triangular_solve 224 output_file: output/ex30_3.out 225 226 TEST*/ 227