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