Actual source code: ex28.c

  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       i,rstart,rend,N=10,num_numfac=5,col[3],k;
  8:   Mat            A[5],F;
  9:   Vec            u,x,b;
 10:   PetscMPIInt    rank;
 11:   PetscScalar    value[3];
 12:   PetscReal      norm,tol=100*PETSC_MACHINE_EPSILON;
 13:   IS             perm,iperm;
 14:   MatFactorInfo  info;
 15:   MatFactorType  facttype = MAT_FACTOR_LU;
 16:   char           solvertype[64];
 17:   char           factortype[64];

 19:   PetscInitialize(&argc,&args,(char*)0,help);
 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*(k+1);
 31:     value[1] =  2.0*(k+1);
 32:     value[2] = -1.0*(k+1);
 33:     for (i=rstart; i<rend; i++) {
 34:       col[0] = i-1; col[1] = i; col[2] = i+1;
 35:       if (i == 0) {
 36:         MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);
 37:       } else if (i == N-1) {
 38:         MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);
 39:       } else {
 40:         MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);
 41:       }
 42:     }
 43:     MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);
 44:     MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);
 45:     MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
 46:   }

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

 52:   /* Set rhs vector b */
 53:   VecSet(b,1.0);

 55:   /* Get a symbolic factor F from A[0] */
 56:   PetscStrncpy(solvertype,"petsc",sizeof(solvertype));
 57:   PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,sizeof(solvertype),NULL);
 58:   PetscOptionsGetEnum(NULL,NULL,"-mat_factor_type",MatFactorTypes,(PetscEnum*)&facttype,NULL);

 60:   MatGetFactor(A[0],solvertype,facttype,&F);
 61:   /* test mumps options */
 62: #if defined(PETSC_HAVE_MUMPS)
 63:   MatMumpsSetIcntl(F,7,5);
 64: #endif
 65:   PetscStrncpy(factortype,MatFactorTypes[facttype],sizeof(factortype));
 66:   PetscStrtoupper(solvertype);
 67:   PetscStrtoupper(factortype);
 68:   PetscPrintf(PETSC_COMM_WORLD," %s %s:\n",solvertype,factortype);

 70:   MatFactorInfoInitialize(&info);
 71:   info.fill = 5.0;
 72:   MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);
 73:   switch (facttype) {
 74:   case MAT_FACTOR_LU:
 75:     MatLUFactorSymbolic(F,A[0],perm,iperm,&info);
 76:     break;
 77:   case MAT_FACTOR_ILU:
 78:     MatILUFactorSymbolic(F,A[0],perm,iperm,&info);
 79:     break;
 80:   case MAT_FACTOR_ICC:
 81:     MatICCFactorSymbolic(F,A[0],perm,&info);
 82:     break;
 83:   case MAT_FACTOR_CHOLESKY:
 84:     MatCholeskyFactorSymbolic(F,A[0],perm,&info);
 85:     break;
 86:   default:
 87:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype);
 88:   }

 90:   /* Compute numeric factors using same F, then solve */
 91:   for (k = 0; k < num_numfac; k++) {
 92:     switch (facttype) {
 93:     case MAT_FACTOR_LU:
 94:     case MAT_FACTOR_ILU:
 95:       MatLUFactorNumeric(F,A[k],&info);
 96:       break;
 97:     case MAT_FACTOR_ICC:
 98:     case MAT_FACTOR_CHOLESKY:
 99:       MatCholeskyFactorNumeric(F,A[k],&info);
100:       break;
101:     default:
102:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype);
103:     }

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

108:     /* Check the residual */
109:     MatMult(A[k],x,u);
110:     VecAXPY(u,-1.0,b);
111:     VecNorm(u,NORM_INFINITY,&norm);
112:     if (norm > tol) {
113:       PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the %s numfact and solve: residual %g\n",k,factortype,(double)norm);
114:     }
115:   }

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

131: /*TEST

133:    test:

135:    test:
136:       suffix: 2
137:       args: -mat_solver_type superlu
138:       requires: superlu

140:    test:
141:       suffix: 3
142:       nsize: 2
143:       requires: mumps
144:       args: -mat_solver_type mumps

146:    test:
147:       suffix: 4
148:       args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
149:       requires: cuda

151: TEST*/