Actual source code: ex128.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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";

     11:  #include <petscmat.h>

     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;
     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;

     27:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     28:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     29:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
     30:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
     31:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
     32:   PetscOptionsGetInt(NULL,NULL,"-lf",&lf,NULL);

     34:   MatCreate(PETSC_COMM_SELF,&C);
     35:   MatSetSizes(C,m*n,m*n,m*n,m*n);
     36:   MatSetFromOptions(C);
     37:   MatSetUp(C);

     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)  {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     44:       J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     45:       J = Ii - 1; if (J>=0)  {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     46:       J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     47:       v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
     48:     }
     49:   }
     50:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
     51:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

     53:   MatIsSymmetric(C,0.0,&flg);
     54:   if (!flg) SETERRQ(PETSC_COMM_SELF,1,"C is non-symmetric");
     55:   MatConvert(C,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sC);

     57:   /* Create vectors for error checking */
     58:   MatCreateVecs(C,&x,&b);
     59:   VecDuplicate(x,&y);
     60:   VecDuplicate(x,&ytmp);
     61:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
     62:   PetscRandomSetFromOptions(rdm);
     63:   VecSetRandom(x,rdm);
     64:   MatMult(C,x,b);

     66:   MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);

     68:   /* Compute CHOLESKY or ICC factor sA */
     69:   MatFactorInfoInitialize(&info);

     71:   info.fill          = 1.0;
     72:   info.diagonal_fill = 0;
     73:   info.zeropivot     = 0.0;

     75:   PetscOptionsHasName(NULL,NULL,"-cholesky",&CHOLESKY);
     76:   if (CHOLESKY) {
     77:     PetscPrintf(PETSC_COMM_SELF,"Test CHOLESKY...\n");
     78:     MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sA);
     79:     MatCholeskyFactorSymbolic(sA,sC,row,&info);
     80:   } else {
     81:     PetscPrintf(PETSC_COMM_SELF,"Test ICC...\n");
     82:     info.levels = lf;

     84:     MatGetFactor(sC,MATSOLVERPETSC,MAT_FACTOR_ICC,&sA);
     85:     MatICCFactorSymbolic(sA,sC,row,&info);
     86:   }
     87:   MatCholeskyFactorNumeric(sA,sC,&info);

     89:   /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
     90:   if (CHOLESKY) {
     91:     PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR);
     92:     if (TRIANGULAR) {
     93:       PetscPrintf(PETSC_COMM_SELF,"Test MatForwardSolve...\n");
     94:       MatForwardSolve(sA,b,ytmp);
     95:       PetscPrintf(PETSC_COMM_SELF,"Test MatBackwardSolve...\n");
     96:       MatBackwardSolve(sA,ytmp,y);
     97:       VecAXPY(y,-1.0,x);
     98:       VecNorm(y,NORM_2,&norm2);
     99:       if (norm2 > tol) {
    100:         PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2);
    101:       }
    102:     }
    103:   }

    105:   MatSolve(sA,b,y);
    106:   MatDestroy(&sC);
    107:   MatDestroy(&sA);
    108:   VecAXPY(y,-1.0,x);
    109:   VecNorm(y,NORM_2,&norm2);
    110:   if (lf == -1 && norm2 > tol) {
    111:     PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ:   Cholesky/ICC levels %D, residual %g\n",lf,(double)norm2);
    112:   }

    114:   /* Free data structures */
    115:   MatDestroy(&C);
    116:   ISDestroy(&row);
    117:   ISDestroy(&col);
    118:   PetscRandomDestroy(&rdm);
    119:   VecDestroy(&x);
    120:   VecDestroy(&y);
    121:   VecDestroy(&ytmp);
    122:   VecDestroy(&b);
    123:   PetscFinalize();
    124:   return ierr;
    125: }

    128: /*TEST

    130:    test:
    131:       output_file: output/ex128.out

    133:    test:
    134:       suffix: 2
    135:       args: -cholesky -triangular_solve

    137: TEST*/