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,sizeof(solvertype),&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*/