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;
 11:   PetscMPIInt    rank;
 12:   PetscScalar    value[3];
 13:   PetscReal      norm,tol=100*PETSC_MACHINE_EPSILON;
 14:   IS             perm,iperm;
 15:   MatFactorInfo  info;
 16:   MatFactorType  facttype = MAT_FACTOR_LU;
 17:   char           solvertype[64];
 18:   char           factortype[64];

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

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

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

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

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

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

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

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

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

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

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

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

132: /*TEST

134:    test:

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

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

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

152: TEST*/