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