1 static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n"; 2 3 #include <petscmat.h> 4 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]; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 23 /* Create and assemble matrices, all have same data structure */ 24 for (k=0; k<num_numfac; k++) { 25 PetscCall(MatCreate(PETSC_COMM_WORLD,&A[k])); 26 PetscCall(MatSetSizes(A[k],PETSC_DECIDE,PETSC_DECIDE,N,N)); 27 PetscCall(MatSetFromOptions(A[k])); 28 PetscCall(MatSetUp(A[k])); 29 PetscCall(MatGetOwnershipRange(A[k],&rstart,&rend)); 30 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 PetscCall(MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES)); 38 } else if (i == N-1) { 39 PetscCall(MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES)); 40 } else { 41 PetscCall(MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES)); 42 } 43 } 44 PetscCall(MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY)); 45 PetscCall(MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 47 } 48 49 /* Create vectors */ 50 PetscCall(MatCreateVecs(A[0],&x,&b)); 51 PetscCall(VecDuplicate(x,&u)); 52 53 /* Set rhs vector b */ 54 PetscCall(VecSet(b,1.0)); 55 56 /* Get a symbolic factor F from A[0] */ 57 PetscCall(PetscStrncpy(solvertype,"petsc",sizeof(solvertype))); 58 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,sizeof(solvertype),NULL)); 59 PetscCall(PetscOptionsGetEnum(NULL,NULL,"-mat_factor_type",MatFactorTypes,(PetscEnum*)&facttype,NULL)); 60 61 PetscCall(MatGetFactor(A[0],solvertype,facttype,&F)); 62 /* test mumps options */ 63 #if defined(PETSC_HAVE_MUMPS) 64 PetscCall(MatMumpsSetIcntl(F,7,5)); 65 #endif 66 PetscCall(PetscStrncpy(factortype,MatFactorTypes[facttype],sizeof(factortype))); 67 PetscCall(PetscStrtoupper(solvertype)); 68 PetscCall(PetscStrtoupper(factortype)); 69 PetscCall(PetscPrintf(PETSC_COMM_WORLD," %s %s:\n",solvertype,factortype)); 70 71 PetscCall(MatFactorInfoInitialize(&info)); 72 info.fill = 5.0; 73 PetscCall(MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm)); 74 switch (facttype) { 75 case MAT_FACTOR_LU: 76 PetscCall(MatLUFactorSymbolic(F,A[0],perm,iperm,&info)); 77 break; 78 case MAT_FACTOR_ILU: 79 PetscCall(MatILUFactorSymbolic(F,A[0],perm,iperm,&info)); 80 break; 81 case MAT_FACTOR_ICC: 82 PetscCall(MatICCFactorSymbolic(F,A[0],perm,&info)); 83 break; 84 case MAT_FACTOR_CHOLESKY: 85 PetscCall(MatCholeskyFactorSymbolic(F,A[0],perm,&info)); 86 break; 87 default: 88 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype); 89 } 90 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 PetscCall(MatLUFactorNumeric(F,A[k],&info)); 97 break; 98 case MAT_FACTOR_ICC: 99 case MAT_FACTOR_CHOLESKY: 100 PetscCall(MatCholeskyFactorNumeric(F,A[k],&info)); 101 break; 102 default: 103 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype); 104 } 105 106 /* Solve A[k] * x = b */ 107 PetscCall(MatSolve(F,b,x)); 108 109 /* Check the residual */ 110 PetscCall(MatMult(A[k],x,u)); 111 PetscCall(VecAXPY(u,-1.0,b)); 112 PetscCall(VecNorm(u,NORM_INFINITY,&norm)); 113 if (norm > tol) { 114 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the %s numfact and solve: residual %g\n",k,factortype,(double)norm)); 115 } 116 } 117 118 /* Free data structures */ 119 for (k=0; k<num_numfac; k++) { 120 PetscCall(MatDestroy(&A[k])); 121 } 122 PetscCall(MatDestroy(&F)); 123 PetscCall(ISDestroy(&perm)); 124 PetscCall(ISDestroy(&iperm)); 125 PetscCall(VecDestroy(&x)); 126 PetscCall(VecDestroy(&b)); 127 PetscCall(VecDestroy(&u)); 128 PetscCall(PetscFinalize()); 129 return 0; 130 } 131 132 /*TEST 133 134 test: 135 136 test: 137 suffix: 2 138 args: -mat_solver_type superlu 139 requires: superlu 140 141 test: 142 suffix: 3 143 nsize: 2 144 requires: mumps 145 args: -mat_solver_type mumps 146 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 151 152 TEST*/ 153