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

      1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n";

      3:  #include <petscmat.h>

      5: int main(int argc,char **args)
      6: {
      7:   PetscInt       ipack,i,rstart,rend,N=10,num_numfac=5,col[3],k;
      8:   Mat            A[5],F;
      9:   Vec            u,x,b;
     11:   PetscMPIInt    rank;
     12:   PetscScalar    value[3];
     13:   PetscReal      norm,tol=100*PETSC_MACHINE_EPSILON;
     14:   IS             perm,iperm;
     15:   MatFactorInfo  info;
     16:   char           solvertype[64]="petsc";
     17:   PetscBool      flg,flg_superlu,flg_mumps;

     19:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     20:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

     22:   /* Create and assemble matrices, all have same data structure */
     23:   for (k=0; k<num_numfac; k++) {
     24:     MatCreate(PETSC_COMM_WORLD,&A[k]);
     25:     MatSetSizes(A[k],PETSC_DECIDE,PETSC_DECIDE,N,N);
     26:     MatSetFromOptions(A[k]);
     27:     MatSetUp(A[k]);
     28:     MatGetOwnershipRange(A[k],&rstart,&rend);

     30:     value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
     31:     for (i=rstart; i<rend; i++) {
     32:       col[0] = i-1; col[1] = i; col[2] = i+1;
     33:       if (i == 0) {
     34:         MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);
     35:       } else if (i == N-1) {
     36:         MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);
     37:       } else {
     38:         MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);
     39:       }
     40:     }
     41:     MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);
     42:     MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);
     43:     MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
     44:   }

     46:   /* Create vectors */
     47:   MatCreateVecs(A[0],&x,&b);
     48:   VecDuplicate(x,&u);

     50:   /* Set rhs vector b */
     51:   VecSet(b,1.0);

     53:   /* Get a symbolic factor F from A[0] */
     54:   PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,64,&flg);
     55:   PetscStrcmp(solvertype,"superlu",&flg_superlu);
     56:   PetscStrcmp(solvertype,"mumps",&flg_mumps);

     58:   ipack = 0;
     59:   if (flg_superlu) ipack = 1;
     60:   else {
     61:     if (flg_mumps) ipack = 2;
     62:   }

     64:   switch (ipack) {
     65:   case 1:
     66: #if defined(PETSC_HAVE_SUPERLU)
     67:     PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
     68:     MatGetFactor(A[0],MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
     69:     break;
     70: #else
     71:     PetscPrintf(PETSC_COMM_WORLD," SUPERLU is not installed, using PETSC LU\n");
     72: #endif
     73:   case 2:
     74: #if defined(PETSC_HAVE_MUMPS)
     75:     PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
     76:     MatGetFactor(A[0],MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
     77:     {
     78:       /* test mumps options */
     79:       PetscInt icntl_7 = 5;
     80:       MatMumpsSetIcntl(F,7,icntl_7);
     81:     }
     82:     break;
     83: #else
     84:     PetscPrintf(PETSC_COMM_WORLD," MUMPS is not installed, use PETSC LU\n");
     85: #endif
     86:   default:
     87:     PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
     88:     MatGetFactor(A[0],MATSOLVERPETSC,MAT_FACTOR_LU,&F);
     89:   }

     91:   MatFactorInfoInitialize(&info);
     92:   info.fill = 5.0;
     93:   MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);
     94:   MatLUFactorSymbolic(F,A[0],perm,iperm,&info);

     96:   /* Compute numeric factors using same F, then solve */
     97:   for (k = 0; k < num_numfac; k++) {
     98:     /* Get numeric factor of A[k] */
     99:     MatLUFactorNumeric(F,A[k],&info);

    101:     /* Solve A[k] * x = b */
    102:     MatSolve(F,b,x);

    104:     /* Check the residual */
    105:     MatMult(A[k],x,u);
    106:     VecAXPY(u,-1.0,b);
    107:     VecNorm(u,NORM_INFINITY,&norm);
    108:     if (norm > tol) {
    109:       PetscPrintf(PETSC_COMM_WORLD,"%D-the LU numfact and solve: residual %g\n",k,(double)norm);
    110:     }
    111:   }

    113:   /* Free data structures */
    114:   for (k=0; k<num_numfac; k++) {
    115:     MatDestroy(&A[k]);
    116:   }
    117:   MatDestroy(&F);
    118:   ISDestroy(&perm);
    119:   ISDestroy(&iperm);
    120:   VecDestroy(&x);
    121:   VecDestroy(&b);
    122:   VecDestroy(&u);
    123:   PetscFinalize();
    124:   return ierr;
    125: }

    127: /*TEST

    129:    test:

    131:    test:
    132:       suffix: 2
    133:       args: -mat_solver_type superlu
    134:       requires: superlu

    136:    test:
    137:       suffix: 3
    138:       nsize: 2
    139:       requires: mumps
    140:       args: -mat_solver_type mumps

    142: TEST*/