1 2 static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqsbaij format. Modified from ex30.c\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,sC,sA; 16 PetscInt i,j,m = 5,n = 5,Ii,J,lf = 0; 17 PetscErrorCode ierr; 18 PetscBool CHOLESKY=PETSC_FALSE,TRIANGULAR=PETSC_FALSE,flg; 19 PetscScalar v; 20 IS row,col; 21 MatFactorInfo info; 22 Vec x,y,b,ytmp; 23 PetscReal norm2,tol = 100*PETSC_MACHINE_EPSILON; 24 PetscRandom rdm; 25 PetscMPIInt size; 26 27 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 29 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 30 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 31 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 32 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-lf",&lf,NULL)); 33 34 CHKERRQ(MatCreate(PETSC_COMM_SELF,&C)); 35 CHKERRQ(MatSetSizes(C,m*n,m*n,m*n,m*n)); 36 CHKERRQ(MatSetFromOptions(C)); 37 CHKERRQ(MatSetUp(C)); 38 39 /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */ 40 for (i=0; i<m; i++) { 41 for (j=0; j<n; j++) { 42 v = -1.0; Ii = j + n*i; 43 J = Ii - n; if (J>=0) CHKERRQ(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 44 J = Ii + n; if (J<m*n) CHKERRQ(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 45 J = Ii - 1; if (J>=0) CHKERRQ(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 46 J = Ii + 1; if (J<m*n) CHKERRQ(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 47 v = 4.0; CHKERRQ(MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 48 } 49 } 50 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 51 CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 52 53 CHKERRQ(MatIsSymmetric(C,0.0,&flg)); 54 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"C is non-symmetric"); 55 CHKERRQ(MatConvert(C,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sC)); 56 57 /* Create vectors for error checking */ 58 CHKERRQ(MatCreateVecs(C,&x,&b)); 59 CHKERRQ(VecDuplicate(x,&y)); 60 CHKERRQ(VecDuplicate(x,&ytmp)); 61 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 62 CHKERRQ(PetscRandomSetFromOptions(rdm)); 63 CHKERRQ(VecSetRandom(x,rdm)); 64 CHKERRQ(MatMult(C,x,b)); 65 66 CHKERRQ(MatGetOrdering(C,MATORDERINGNATURAL,&row,&col)); 67 68 /* Compute CHOLESKY or ICC factor sA */ 69 CHKERRQ(MatFactorInfoInitialize(&info)); 70 71 info.fill = 1.0; 72 info.diagonal_fill = 0; 73 info.zeropivot = 0.0; 74 75 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-cholesky",&CHOLESKY)); 76 if (CHOLESKY) { 77 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Test CHOLESKY...\n")); 78 CHKERRQ(MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sA)); 79 CHKERRQ(MatCholeskyFactorSymbolic(sA,sC,row,&info)); 80 } else { 81 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Test ICC...\n")); 82 info.levels = lf; 83 84 CHKERRQ(MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_ICC,&sA)); 85 CHKERRQ(MatICCFactorSymbolic(sA,sC,row,&info)); 86 } 87 CHKERRQ(MatCholeskyFactorNumeric(sA,sC,&info)); 88 89 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */ 90 if (CHOLESKY) { 91 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR)); 92 if (TRIANGULAR) { 93 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Test MatForwardSolve...\n")); 94 CHKERRQ(MatForwardSolve(sA,b,ytmp)); 95 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Test MatBackwardSolve...\n")); 96 CHKERRQ(MatBackwardSolve(sA,ytmp,y)); 97 CHKERRQ(VecAXPY(y,-1.0,x)); 98 CHKERRQ(VecNorm(y,NORM_2,&norm2)); 99 if (norm2 > tol) { 100 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2)); 101 } 102 } 103 } 104 105 CHKERRQ(MatSolve(sA,b,y)); 106 CHKERRQ(MatDestroy(&sC)); 107 CHKERRQ(MatDestroy(&sA)); 108 CHKERRQ(VecAXPY(y,-1.0,x)); 109 CHKERRQ(VecNorm(y,NORM_2,&norm2)); 110 if (lf == -1 && norm2 > tol) { 111 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n",lf,(double)norm2)); 112 } 113 114 /* Free data structures */ 115 CHKERRQ(MatDestroy(&C)); 116 CHKERRQ(ISDestroy(&row)); 117 CHKERRQ(ISDestroy(&col)); 118 CHKERRQ(PetscRandomDestroy(&rdm)); 119 CHKERRQ(VecDestroy(&x)); 120 CHKERRQ(VecDestroy(&y)); 121 CHKERRQ(VecDestroy(&ytmp)); 122 CHKERRQ(VecDestroy(&b)); 123 ierr = PetscFinalize(); 124 return ierr; 125 } 126 127 /*TEST 128 129 test: 130 output_file: output/ex128.out 131 132 test: 133 suffix: 2 134 args: -cholesky -triangular_solve 135 136 TEST*/ 137