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 ipack,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 char solvertype[64]="petsc"; 17 PetscBool flg,flg_superlu,flg_mumps; 18 19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 20 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 21 22 /* Create and assemble matrices, all have same data structure */ 23 for (k=0; k<num_numfac; k++) { 24 ierr = MatCreate(PETSC_COMM_WORLD,&A[k]);CHKERRQ(ierr); 25 ierr = MatSetSizes(A[k],PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); 26 ierr = MatSetFromOptions(A[k]);CHKERRQ(ierr); 27 ierr = MatSetUp(A[k]);CHKERRQ(ierr); 28 ierr = MatGetOwnershipRange(A[k],&rstart,&rend);CHKERRQ(ierr); 29 30 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 31 for (i=rstart; i<rend; i++) { 32 col[0] = i-1; col[1] = i; col[2] = i+1; 33 if (i == 0) { 34 ierr = MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);CHKERRQ(ierr); 35 } else if (i == N-1) { 36 ierr = MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 37 } else { 38 ierr = MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 39 } 40 } 41 ierr = MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 42 ierr = MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 43 ierr = MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 44 } 45 46 /* Create vectors */ 47 ierr = MatCreateVecs(A[0],&x,&b);CHKERRQ(ierr); 48 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); 49 50 /* Set rhs vector b */ 51 ierr = VecSet(b,1.0);CHKERRQ(ierr); 52 53 /* Get a symbolic factor F from A[0] */ 54 ierr = PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,sizeof(solvertype),&flg);CHKERRQ(ierr); 55 ierr = PetscStrcmp(solvertype,"superlu",&flg_superlu);CHKERRQ(ierr); 56 ierr = PetscStrcmp(solvertype,"mumps",&flg_mumps);CHKERRQ(ierr); 57 58 ipack = 0; 59 if (flg_superlu) ipack = 1; 60 else { 61 if (flg_mumps) ipack = 2; 62 } 63 64 switch (ipack) { 65 case 1: 66 #if defined(PETSC_HAVE_SUPERLU) 67 ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");CHKERRQ(ierr); 68 ierr = MatGetFactor(A[0],MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 69 break; 70 #else 71 ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU is not installed, using PETSC LU\n");CHKERRQ(ierr); 72 #endif 73 case 2: 74 #if defined(PETSC_HAVE_MUMPS) 75 ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");CHKERRQ(ierr); 76 ierr = MatGetFactor(A[0],MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 77 { 78 /* test mumps options */ 79 PetscInt icntl_7 = 5; 80 ierr = MatMumpsSetIcntl(F,7,icntl_7);CHKERRQ(ierr); 81 } 82 break; 83 #else 84 ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS is not installed, use PETSC LU\n");CHKERRQ(ierr); 85 #endif 86 default: 87 ierr = PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");CHKERRQ(ierr); 88 ierr = MatGetFactor(A[0],MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 89 } 90 91 ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr); 92 info.fill = 5.0; 93 ierr = MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);CHKERRQ(ierr); 94 ierr = MatLUFactorSymbolic(F,A[0],perm,iperm,&info);CHKERRQ(ierr); 95 96 /* Compute numeric factors using same F, then solve */ 97 for (k = 0; k < num_numfac; k++) { 98 /* Get numeric factor of A[k] */ 99 ierr = MatLUFactorNumeric(F,A[k],&info);CHKERRQ(ierr); 100 101 /* Solve A[k] * x = b */ 102 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 103 104 /* Check the residual */ 105 ierr = MatMult(A[k],x,u);CHKERRQ(ierr); 106 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); 107 ierr = VecNorm(u,NORM_INFINITY,&norm);CHKERRQ(ierr); 108 if (norm > tol) { 109 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D-the LU numfact and solve: residual %g\n",k,(double)norm);CHKERRQ(ierr); 110 } 111 } 112 113 /* Free data structures */ 114 for (k=0; k<num_numfac; k++) { 115 ierr = MatDestroy(&A[k]);CHKERRQ(ierr); 116 } 117 ierr = MatDestroy(&F);CHKERRQ(ierr); 118 ierr = ISDestroy(&perm);CHKERRQ(ierr); 119 ierr = ISDestroy(&iperm);CHKERRQ(ierr); 120 ierr = VecDestroy(&x);CHKERRQ(ierr); 121 ierr = VecDestroy(&b);CHKERRQ(ierr); 122 ierr = VecDestroy(&u);CHKERRQ(ierr); 123 ierr = PetscFinalize(); 124 return ierr; 125 } 126 127 /*TEST 128 129 test: 130 131 test: 132 suffix: 2 133 args: -mat_solver_type superlu 134 requires: superlu 135 136 test: 137 suffix: 3 138 nsize: 2 139 requires: mumps 140 args: -mat_solver_type mumps 141 142 TEST*/ 143