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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 29 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 30 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsGetInt(NULL,NULL,"-lf",&lf,NULL);CHKERRQ(ierr); 33 34 ierr = MatCreate(PETSC_COMM_SELF,&C);CHKERRQ(ierr); 35 ierr = MatSetSizes(C,m*n,m*n,m*n,m*n);CHKERRQ(ierr); 36 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 37 ierr = MatSetUp(C);CHKERRQ(ierr); 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) {ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 44 J = Ii + n; if (J<m*n) {ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 45 J = Ii - 1; if (J>=0) {ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 46 J = Ii + 1; if (J<m*n) {ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} 47 v = 4.0; ierr = MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr); 48 } 49 } 50 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 51 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 52 53 ierr = MatIsSymmetric(C,0.0,&flg);CHKERRQ(ierr); 54 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"C is non-symmetric"); 55 ierr = MatConvert(C,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sC);CHKERRQ(ierr); 56 57 /* Create vectors for error checking */ 58 ierr = MatCreateVecs(C,&x,&b);CHKERRQ(ierr); 59 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 60 ierr = VecDuplicate(x,&ytmp);CHKERRQ(ierr); 61 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rdm);CHKERRQ(ierr); 62 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 63 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 64 ierr = MatMult(C,x,b);CHKERRQ(ierr); 65 66 ierr = MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);CHKERRQ(ierr); 67 68 /* Compute CHOLESKY or ICC factor sA */ 69 ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr); 70 71 info.fill = 1.0; 72 info.diagonal_fill = 0; 73 info.zeropivot = 0.0; 74 75 ierr = PetscOptionsHasName(NULL,NULL,"-cholesky",&CHOLESKY);CHKERRQ(ierr); 76 if (CHOLESKY) { 77 ierr = PetscPrintf(PETSC_COMM_SELF,"Test CHOLESKY...\n");CHKERRQ(ierr); 78 ierr = MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sA);CHKERRQ(ierr); 79 ierr = MatCholeskyFactorSymbolic(sA,sC,row,&info);CHKERRQ(ierr); 80 } else { 81 ierr = PetscPrintf(PETSC_COMM_SELF,"Test ICC...\n");CHKERRQ(ierr); 82 info.levels = lf; 83 84 ierr = MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_ICC,&sA);CHKERRQ(ierr); 85 ierr = MatICCFactorSymbolic(sA,sC,row,&info);CHKERRQ(ierr); 86 } 87 ierr = MatCholeskyFactorNumeric(sA,sC,&info);CHKERRQ(ierr); 88 89 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */ 90 if (CHOLESKY) { 91 ierr = PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR);CHKERRQ(ierr); 92 if (TRIANGULAR) { 93 ierr = PetscPrintf(PETSC_COMM_SELF,"Test MatForwardSolve...\n");CHKERRQ(ierr); 94 ierr = MatForwardSolve(sA,b,ytmp);CHKERRQ(ierr); 95 ierr = PetscPrintf(PETSC_COMM_SELF,"Test MatBackwardSolve...\n");CHKERRQ(ierr); 96 ierr = MatBackwardSolve(sA,ytmp,y);CHKERRQ(ierr); 97 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 98 ierr = VecNorm(y,NORM_2,&norm2);CHKERRQ(ierr); 99 if (norm2 > tol) { 100 ierr = PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2);CHKERRQ(ierr); 101 } 102 } 103 } 104 105 ierr = MatSolve(sA,b,y);CHKERRQ(ierr); 106 ierr = MatDestroy(&sC);CHKERRQ(ierr); 107 ierr = MatDestroy(&sA);CHKERRQ(ierr); 108 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 109 ierr = VecNorm(y,NORM_2,&norm2);CHKERRQ(ierr); 110 if (lf == -1 && norm2 > tol) { 111 ierr = PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %D, residual %g\n",lf,(double)norm2);CHKERRQ(ierr); 112 } 113 114 /* Free data structures */ 115 ierr = MatDestroy(&C);CHKERRQ(ierr); 116 ierr = ISDestroy(&row);CHKERRQ(ierr); 117 ierr = ISDestroy(&col);CHKERRQ(ierr); 118 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 119 ierr = VecDestroy(&x);CHKERRQ(ierr); 120 ierr = VecDestroy(&y);CHKERRQ(ierr); 121 ierr = VecDestroy(&ytmp);CHKERRQ(ierr); 122 ierr = VecDestroy(&b);CHKERRQ(ierr); 123 ierr = PetscFinalize(); 124 return ierr; 125 } 126 127 128 /*TEST 129 130 test: 131 output_file: output/ex128.out 132 133 test: 134 suffix: 2 135 args: -cholesky -triangular_solve 136 137 TEST*/ 138