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