xref: /petsc/src/ksp/pc/impls/factor/cholesky/cholesky.c (revision c7393fdb8164c15777eb367667474a1e4da0d712)
1dba47a55SKris Buschelman #define PETSCKSP_DLL
29b54502bSHong Zhang 
39b54502bSHong Zhang /*
49b54502bSHong Zhang    Defines a direct factorization preconditioner for any Mat implementation
59b54502bSHong Zhang    Note: this need not be consided a preconditioner since it supplies
69b54502bSHong Zhang          a direct solver.
79b54502bSHong Zhang */
86356e834SBarry Smith #include "private/pcimpl.h"                /*I "petscpc.h" I*/
99b54502bSHong Zhang 
109b54502bSHong Zhang typedef struct {
119b54502bSHong Zhang   Mat             fact;             /* factored matrix */
129b54502bSHong Zhang   PetscReal       actualfill;       /* actual fill in factor */
139b54502bSHong Zhang   PetscTruth      inplace;          /* flag indicating in-place factorization */
149b54502bSHong Zhang   IS              row,col;          /* index sets used for reordering */
159b54502bSHong Zhang   MatOrderingType ordering;         /* matrix ordering */
169b54502bSHong Zhang   PetscTruth      reuseordering;    /* reuses previous reordering computed */
179b54502bSHong Zhang   PetscTruth      reusefill;        /* reuse fill from previous Cholesky */
189b54502bSHong Zhang   MatFactorInfo   info;
19*c7393fdbSBarry Smith   MatSolverPackage   solvertype;
209b54502bSHong Zhang } PC_Cholesky;
219b54502bSHong Zhang 
229b54502bSHong Zhang EXTERN_C_BEGIN
239b54502bSHong Zhang #undef __FUNCT__
24*c7393fdbSBarry Smith #define __FUNCT__ "PCFactorSetMatSolverPackage_LU"
25*c7393fdbSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatSolverPackage_Cholesky(PC pc,const MatSolverPackage stype)
26c62fe420SBarry Smith {
27c62fe420SBarry Smith   PetscErrorCode ierr;
28c62fe420SBarry Smith   PC_Cholesky    *choleksy = (PC_Cholesky*)pc->data;
29c62fe420SBarry Smith 
30c62fe420SBarry Smith   PetscFunctionBegin;
31c62fe420SBarry Smith   ierr = PetscStrallocpy(stype,&choleksy->solvertype);CHKERRQ(ierr);
32c62fe420SBarry Smith   PetscFunctionReturn(0);
33c62fe420SBarry Smith }
34c62fe420SBarry Smith EXTERN_C_END
35c62fe420SBarry Smith 
36c62fe420SBarry Smith EXTERN_C_BEGIN
37c62fe420SBarry Smith #undef __FUNCT__
38afaefe49SHong Zhang #define __FUNCT__ "PCFactorSetZeroPivot_Cholesky"
39dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_Cholesky(PC pc,PetscReal z)
40afaefe49SHong Zhang {
41afaefe49SHong Zhang   PC_Cholesky *ch;
42afaefe49SHong Zhang 
43afaefe49SHong Zhang   PetscFunctionBegin;
44afaefe49SHong Zhang   ch                 = (PC_Cholesky*)pc->data;
45afaefe49SHong Zhang   ch->info.zeropivot = z;
46afaefe49SHong Zhang   PetscFunctionReturn(0);
47afaefe49SHong Zhang }
48afaefe49SHong Zhang EXTERN_C_END
49afaefe49SHong Zhang 
50afaefe49SHong Zhang EXTERN_C_BEGIN
51afaefe49SHong Zhang #undef __FUNCT__
52afaefe49SHong Zhang #define __FUNCT__ "PCFactorSetShiftNonzero_Cholesky"
53dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftNonzero_Cholesky(PC pc,PetscReal shift)
54afaefe49SHong Zhang {
55afaefe49SHong Zhang   PC_Cholesky *dir;
56afaefe49SHong Zhang 
57afaefe49SHong Zhang   PetscFunctionBegin;
58afaefe49SHong Zhang   dir = (PC_Cholesky*)pc->data;
59afaefe49SHong Zhang   if (shift == (PetscReal) PETSC_DECIDE) {
60afaefe49SHong Zhang     dir->info.shiftnz = 1.e-12;
61afaefe49SHong Zhang   } else {
62afaefe49SHong Zhang     dir->info.shiftnz = shift;
63afaefe49SHong Zhang   }
64afaefe49SHong Zhang   PetscFunctionReturn(0);
65afaefe49SHong Zhang }
66afaefe49SHong Zhang EXTERN_C_END
67afaefe49SHong Zhang 
68afaefe49SHong Zhang EXTERN_C_BEGIN
69afaefe49SHong Zhang #undef __FUNCT__
70afaefe49SHong Zhang #define __FUNCT__ "PCFactorSetShiftPd_Cholesky"
71dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftPd_Cholesky(PC pc,PetscTruth shift)
72afaefe49SHong Zhang {
73afaefe49SHong Zhang   PC_Cholesky *dir;
74afaefe49SHong Zhang 
75afaefe49SHong Zhang   PetscFunctionBegin;
76afaefe49SHong Zhang   dir = (PC_Cholesky*)pc->data;
77fbf22428SSatish Balay   if (shift) {
78fbf22428SSatish Balay     dir->info.shift_fraction = 0.0;
79fbf22428SSatish Balay     dir->info.shiftpd = 1.0;
80fbf22428SSatish Balay   } else {
81fbf22428SSatish Balay     dir->info.shiftpd = 0.0;
82fbf22428SSatish Balay   }
83afaefe49SHong Zhang   PetscFunctionReturn(0);
84afaefe49SHong Zhang }
85afaefe49SHong Zhang EXTERN_C_END
86afaefe49SHong Zhang 
87afaefe49SHong Zhang EXTERN_C_BEGIN
88afaefe49SHong Zhang #undef __FUNCT__
892401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseOrdering_Cholesky"
902401956bSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
919b54502bSHong Zhang {
929b54502bSHong Zhang   PC_Cholesky *lu;
939b54502bSHong Zhang 
949b54502bSHong Zhang   PetscFunctionBegin;
959b54502bSHong Zhang   lu               = (PC_Cholesky*)pc->data;
969b54502bSHong Zhang   lu->reuseordering = flag;
979b54502bSHong Zhang   PetscFunctionReturn(0);
989b54502bSHong Zhang }
999b54502bSHong Zhang EXTERN_C_END
1009b54502bSHong Zhang 
1019b54502bSHong Zhang EXTERN_C_BEGIN
1029b54502bSHong Zhang #undef __FUNCT__
1032401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseFill_Cholesky"
1042401956bSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill_Cholesky(PC pc,PetscTruth flag)
1059b54502bSHong Zhang {
1069b54502bSHong Zhang   PC_Cholesky *lu;
1079b54502bSHong Zhang 
1089b54502bSHong Zhang   PetscFunctionBegin;
1099b54502bSHong Zhang   lu = (PC_Cholesky*)pc->data;
1109b54502bSHong Zhang   lu->reusefill = flag;
1119b54502bSHong Zhang   PetscFunctionReturn(0);
1129b54502bSHong Zhang }
1139b54502bSHong Zhang EXTERN_C_END
1149b54502bSHong Zhang 
1159b54502bSHong Zhang #undef __FUNCT__
1169b54502bSHong Zhang #define __FUNCT__ "PCSetFromOptions_Cholesky"
1179b54502bSHong Zhang static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
1189b54502bSHong Zhang {
1199b54502bSHong Zhang   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
1209b54502bSHong Zhang   PetscErrorCode ierr;
1219b54502bSHong Zhang   PetscTruth     flg;
1225c9eb25fSBarry Smith   char           tname[256], solvertype[64];
1239b54502bSHong Zhang   PetscFList     ordlist;
1249b54502bSHong Zhang 
1259b54502bSHong Zhang   PetscFunctionBegin;
1269b54502bSHong Zhang   ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);
1279b54502bSHong Zhang   ierr = PetscOptionsHead("Cholesky options");CHKERRQ(ierr);
1282401956bSBarry Smith     ierr = PetscOptionsName("-pc_factor_in_place","Form Cholesky in the same memory as the matrix","PCFactorSetUseInPlace",&flg);CHKERRQ(ierr);
1299b54502bSHong Zhang     if (flg) {
1302401956bSBarry Smith       ierr = PCFactorSetUseInPlace(pc);CHKERRQ(ierr);
1319b54502bSHong Zhang     }
13255ba2a51SBarry Smith     ierr = PetscOptionsReal("-pc_factor_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCFactorSetFill",lu->info.fill,&lu->info.fill,0);CHKERRQ(ierr);
1339b54502bSHong Zhang 
1342401956bSBarry Smith     ierr = PetscOptionsName("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",&flg);CHKERRQ(ierr);
1359b54502bSHong Zhang     if (flg) {
1362401956bSBarry Smith       ierr = PCFactorSetReuseFill(pc,PETSC_TRUE);CHKERRQ(ierr);
1379b54502bSHong Zhang     }
1382401956bSBarry Smith     ierr = PetscOptionsName("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",&flg);CHKERRQ(ierr);
1399b54502bSHong Zhang     if (flg) {
1402401956bSBarry Smith       ierr = PCFactorSetReuseOrdering(pc,PETSC_TRUE);CHKERRQ(ierr);
1419b54502bSHong Zhang     }
1429b54502bSHong Zhang 
1439b54502bSHong Zhang     ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr);
144e5a9bf91SBarry Smith     ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCFactorSetMatOrderingType",ordlist,lu->ordering,tname,256,&flg);CHKERRQ(ierr);
1459b54502bSHong Zhang     if (flg) {
146e5a9bf91SBarry Smith       ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr);
1479b54502bSHong Zhang     }
1485c9eb25fSBarry Smith 
1495c9eb25fSBarry Smith     /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
150*c7393fdbSBarry Smith     ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific Cholesky solver to use","MatGetFactor",lu->solvertype,solvertype,64,&flg);CHKERRQ(ierr);
1515c9eb25fSBarry Smith     if (flg) {
152*c7393fdbSBarry Smith       ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr);
1535c9eb25fSBarry Smith     }
1545c9eb25fSBarry Smith 
1559f95998fSHong Zhang     ierr = PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);CHKERRQ(ierr);
1569b54502bSHong Zhang     if (flg) {
157afaefe49SHong Zhang       ierr = PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);CHKERRQ(ierr);
1589b54502bSHong Zhang     }
1599f95998fSHong Zhang     ierr = PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",lu->info.shiftnz,&lu->info.shiftnz,0);CHKERRQ(ierr);
1609f95998fSHong Zhang     ierr = PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);CHKERRQ(ierr);
1619b54502bSHong Zhang     if (flg) {
162afaefe49SHong Zhang       ierr = PCFactorSetShiftPd(pc,PETSC_TRUE);CHKERRQ(ierr);
1639b54502bSHong Zhang     }
164ee45ca4aSHong Zhang     ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",lu->info.zeropivot,&lu->info.zeropivot,0);CHKERRQ(ierr);
1659b54502bSHong Zhang 
1669b54502bSHong Zhang   ierr = PetscOptionsTail();CHKERRQ(ierr);
1679b54502bSHong Zhang   PetscFunctionReturn(0);
1689b54502bSHong Zhang }
1699b54502bSHong Zhang 
1709b54502bSHong Zhang #undef __FUNCT__
1719b54502bSHong Zhang #define __FUNCT__ "PCView_Cholesky"
1729b54502bSHong Zhang static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
1739b54502bSHong Zhang {
1749b54502bSHong Zhang   PC_Cholesky    *lu = (PC_Cholesky*)pc->data;
1759b54502bSHong Zhang   PetscErrorCode ierr;
1769b54502bSHong Zhang   PetscTruth     iascii,isstring;
1779b54502bSHong Zhang 
1789b54502bSHong Zhang   PetscFunctionBegin;
1799b54502bSHong Zhang   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
1809b54502bSHong Zhang   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
1819b54502bSHong Zhang   if (iascii) {
1829b54502bSHong Zhang 
1839b54502bSHong Zhang     if (lu->inplace) {ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorization\n");CHKERRQ(ierr);}
1849b54502bSHong Zhang     else             {ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorization\n");CHKERRQ(ierr);}
1859b54502bSHong Zhang     ierr = PetscViewerASCIIPrintf(viewer,"    matrix ordering: %s\n",lu->ordering);CHKERRQ(ierr);
1869b54502bSHong Zhang     if (lu->fact) {
187f3a39becSBarry Smith       ierr = PetscViewerASCIIPrintf(viewer,"  Cholesky: factor fill ratio needed %G\n",lu->actualfill);CHKERRQ(ierr);
188f3a39becSBarry Smith       ierr = PetscViewerASCIIPrintf(viewer,"       Factored matrix follows\n");CHKERRQ(ierr);
189f3a39becSBarry Smith       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
190f3a39becSBarry Smith       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
191f3a39becSBarry Smith       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
192f3a39becSBarry Smith       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
1939b54502bSHong Zhang       ierr = MatView(lu->fact,viewer);CHKERRQ(ierr);
1949b54502bSHong Zhang       ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
195f3a39becSBarry Smith       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
196f3a39becSBarry Smith       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
197f3a39becSBarry Smith       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1989b54502bSHong Zhang     }
1999b54502bSHong Zhang     if (lu->reusefill)    {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");CHKERRQ(ierr);}
2009b54502bSHong Zhang     if (lu->reuseordering) {ierr = PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");CHKERRQ(ierr);}
2019b54502bSHong Zhang   } else if (isstring) {
2029b54502bSHong Zhang     ierr = PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);CHKERRQ(ierr);CHKERRQ(ierr);
2039b54502bSHong Zhang   } else {
2042401956bSBarry Smith     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name);
2059b54502bSHong Zhang   }
2069b54502bSHong Zhang   PetscFunctionReturn(0);
2079b54502bSHong Zhang }
2089b54502bSHong Zhang 
2099b54502bSHong Zhang #undef __FUNCT__
210a4fd02acSBarry Smith #define __FUNCT__ "PCFactorGetMatrix_Cholesky"
211a4fd02acSBarry Smith static PetscErrorCode PCFactorGetMatrix_Cholesky(PC pc,Mat *mat)
2129b54502bSHong Zhang {
2139b54502bSHong Zhang   PC_Cholesky *dir = (PC_Cholesky*)pc->data;
2149b54502bSHong Zhang 
2159b54502bSHong Zhang   PetscFunctionBegin;
2169b54502bSHong Zhang   if (!dir->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
2179b54502bSHong Zhang   *mat = dir->fact;
2189b54502bSHong Zhang   PetscFunctionReturn(0);
2199b54502bSHong Zhang }
2209b54502bSHong Zhang 
2219b54502bSHong Zhang #undef __FUNCT__
2229b54502bSHong Zhang #define __FUNCT__ "PCSetUp_Cholesky"
2239b54502bSHong Zhang static PetscErrorCode PCSetUp_Cholesky(PC pc)
2249b54502bSHong Zhang {
2259b54502bSHong Zhang   PetscErrorCode ierr;
2269b54502bSHong Zhang   PetscTruth     flg;
2279b54502bSHong Zhang   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
2289b54502bSHong Zhang 
2299b54502bSHong Zhang   PetscFunctionBegin;
2309b54502bSHong Zhang   if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
2319b54502bSHong Zhang 
2329b54502bSHong Zhang   if (dir->inplace) {
2339b54502bSHong Zhang     if (dir->row && dir->col && (dir->row != dir->col)) {
2349b54502bSHong Zhang       ierr = ISDestroy(dir->row);CHKERRQ(ierr);
2359b54502bSHong Zhang       dir->row = 0;
2369b54502bSHong Zhang     }
2379b54502bSHong Zhang     if (dir->col) {
2389b54502bSHong Zhang       ierr = ISDestroy(dir->col);CHKERRQ(ierr);
2399b54502bSHong Zhang       dir->col = 0;
2409b54502bSHong Zhang     }
2419b54502bSHong Zhang     ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
2429b54502bSHong Zhang     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
2439b54502bSHong Zhang       ierr = ISDestroy(dir->col);CHKERRQ(ierr);
2449b54502bSHong Zhang       dir->col=0;
2459b54502bSHong Zhang     }
246efee365bSSatish Balay     if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
2479b54502bSHong Zhang     ierr = MatCholeskyFactor(pc->pmat,dir->row,&dir->info);CHKERRQ(ierr);
2489b54502bSHong Zhang     dir->fact = pc->pmat;
2499b54502bSHong Zhang   } else {
2509b54502bSHong Zhang     MatInfo info;
2519b54502bSHong Zhang     if (!pc->setupcalled) {
2529b54502bSHong Zhang       ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
2539bfd6278SHong Zhang       /* check if dir->row == dir->col */
2549bfd6278SHong Zhang       ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr);
2559bfd6278SHong Zhang       if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
2569bfd6278SHong Zhang       ierr = ISDestroy(dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */
2579b54502bSHong Zhang       dir->col=0;
2589bfd6278SHong Zhang 
2597adad957SLisandro Dalcin       ierr = PetscOptionsHasName(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);CHKERRQ(ierr);
2609b54502bSHong Zhang       if (flg) {
2619b54502bSHong Zhang         PetscReal tol = 1.e-10;
2627adad957SLisandro Dalcin         ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr);
2639b54502bSHong Zhang         ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr);
2649b54502bSHong Zhang       }
265efee365bSSatish Balay       if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
2665c9eb25fSBarry Smith       ierr = MatGetFactor(pc->pmat,dir->solvertype,MAT_FACTOR_CHOLESKY,&dir->fact);CHKERRQ(ierr);
2679b54502bSHong Zhang       ierr = MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);CHKERRQ(ierr);
2689b54502bSHong Zhang       ierr = MatGetInfo(dir->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2699b54502bSHong Zhang       dir->actualfill = info.fill_ratio_needed;
27052e6d16bSBarry Smith       ierr = PetscLogObjectParent(pc,dir->fact);CHKERRQ(ierr);
2719b54502bSHong Zhang     } else if (pc->flag != SAME_NONZERO_PATTERN) {
2729b54502bSHong Zhang       if (!dir->reuseordering) {
2739b54502bSHong Zhang         if (dir->row && dir->col && (dir->row != dir->col)) {
2749b54502bSHong Zhang           ierr = ISDestroy(dir->row);CHKERRQ(ierr);
2759b54502bSHong Zhang           dir->row = 0;
2769b54502bSHong Zhang         }
2779b54502bSHong Zhang         if (dir->col) {
2789b54502bSHong Zhang           ierr = ISDestroy(dir->col);CHKERRQ(ierr);
2799b54502bSHong Zhang           dir->col =0;
2809b54502bSHong Zhang         }
2819b54502bSHong Zhang         ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
2829b54502bSHong Zhang         if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
2839b54502bSHong Zhang           ierr = ISDestroy(dir->col);CHKERRQ(ierr);
2849b54502bSHong Zhang           dir->col=0;
2859b54502bSHong Zhang         }
2867adad957SLisandro Dalcin         ierr = PetscOptionsHasName(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);CHKERRQ(ierr);
2879b54502bSHong Zhang         if (flg) {
2889b54502bSHong Zhang           PetscReal tol = 1.e-10;
2897adad957SLisandro Dalcin           ierr = PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);CHKERRQ(ierr);
2909b54502bSHong Zhang           ierr = MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);CHKERRQ(ierr);
2919b54502bSHong Zhang         }
292efee365bSSatish Balay         if (dir->row) {ierr = PetscLogObjectParent(pc,dir->row);CHKERRQ(ierr);}
2939b54502bSHong Zhang       }
2949b54502bSHong Zhang       ierr = MatDestroy(dir->fact);CHKERRQ(ierr);
2955c9eb25fSBarry Smith       ierr = MatGetFactor(pc->pmat,dir->solvertype,MAT_FACTOR_CHOLESKY,&dir->fact);CHKERRQ(ierr);
2969b54502bSHong Zhang       ierr = MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);CHKERRQ(ierr);
2979b54502bSHong Zhang       ierr = MatGetInfo(dir->fact,MAT_LOCAL,&info);CHKERRQ(ierr);
2989b54502bSHong Zhang       dir->actualfill = info.fill_ratio_needed;
29952e6d16bSBarry Smith       ierr = PetscLogObjectParent(pc,dir->fact);CHKERRQ(ierr);
3009b54502bSHong Zhang     }
3019b54502bSHong Zhang     ierr = MatCholeskyFactorNumeric(pc->pmat,&dir->info,&dir->fact);CHKERRQ(ierr);
3029b54502bSHong Zhang   }
3039b54502bSHong Zhang   PetscFunctionReturn(0);
3049b54502bSHong Zhang }
3059b54502bSHong Zhang 
3069b54502bSHong Zhang #undef __FUNCT__
3079b54502bSHong Zhang #define __FUNCT__ "PCDestroy_Cholesky"
3089b54502bSHong Zhang static PetscErrorCode PCDestroy_Cholesky(PC pc)
3099b54502bSHong Zhang {
3109b54502bSHong Zhang   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
3119b54502bSHong Zhang   PetscErrorCode ierr;
3129b54502bSHong Zhang 
3139b54502bSHong Zhang   PetscFunctionBegin;
3149b54502bSHong Zhang   if (!dir->inplace && dir->fact) {ierr = MatDestroy(dir->fact);CHKERRQ(ierr);}
3159b54502bSHong Zhang   if (dir->row) {ierr = ISDestroy(dir->row);CHKERRQ(ierr);}
3169b54502bSHong Zhang   if (dir->col) {ierr = ISDestroy(dir->col);CHKERRQ(ierr);}
3179b54502bSHong Zhang   ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr);
3185c9eb25fSBarry Smith   ierr = PetscStrfree(dir->solvertype);CHKERRQ(ierr);
3199b54502bSHong Zhang   ierr = PetscFree(dir);CHKERRQ(ierr);
3209b54502bSHong Zhang   PetscFunctionReturn(0);
3219b54502bSHong Zhang }
3229b54502bSHong Zhang 
3239b54502bSHong Zhang #undef __FUNCT__
3249b54502bSHong Zhang #define __FUNCT__ "PCApply_Cholesky"
3259b54502bSHong Zhang static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
3269b54502bSHong Zhang {
3279b54502bSHong Zhang   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
3289b54502bSHong Zhang   PetscErrorCode ierr;
3299b54502bSHong Zhang 
3309b54502bSHong Zhang   PetscFunctionBegin;
3319b54502bSHong Zhang   if (dir->inplace) {ierr = MatSolve(pc->pmat,x,y);CHKERRQ(ierr);}
3329b54502bSHong Zhang   else              {ierr = MatSolve(dir->fact,x,y);CHKERRQ(ierr);}
3339b54502bSHong Zhang   PetscFunctionReturn(0);
3349b54502bSHong Zhang }
3359b54502bSHong Zhang 
3369b54502bSHong Zhang #undef __FUNCT__
3379b54502bSHong Zhang #define __FUNCT__ "PCApplyTranspose_Cholesky"
3389b54502bSHong Zhang static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
3399b54502bSHong Zhang {
3409b54502bSHong Zhang   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
3419b54502bSHong Zhang   PetscErrorCode ierr;
3429b54502bSHong Zhang 
3439b54502bSHong Zhang   PetscFunctionBegin;
3449b54502bSHong Zhang   if (dir->inplace) {ierr = MatSolveTranspose(pc->pmat,x,y);CHKERRQ(ierr);}
3459b54502bSHong Zhang   else              {ierr = MatSolveTranspose(dir->fact,x,y);CHKERRQ(ierr);}
3469b54502bSHong Zhang   PetscFunctionReturn(0);
3479b54502bSHong Zhang }
3489b54502bSHong Zhang 
3499b54502bSHong Zhang /* -----------------------------------------------------------------------------------*/
3509b54502bSHong Zhang 
3519b54502bSHong Zhang EXTERN_C_BEGIN
3529b54502bSHong Zhang #undef __FUNCT__
35355ba2a51SBarry Smith #define __FUNCT__ "PCFactorSetFill_Cholesky"
35455ba2a51SBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetFill_Cholesky(PC pc,PetscReal fill)
3559b54502bSHong Zhang {
3569b54502bSHong Zhang   PC_Cholesky *dir;
3579b54502bSHong Zhang 
3589b54502bSHong Zhang   PetscFunctionBegin;
3599b54502bSHong Zhang   dir = (PC_Cholesky*)pc->data;
3609b54502bSHong Zhang   dir->info.fill = fill;
3619b54502bSHong Zhang   PetscFunctionReturn(0);
3629b54502bSHong Zhang }
3639b54502bSHong Zhang EXTERN_C_END
3649b54502bSHong Zhang 
3659b54502bSHong Zhang EXTERN_C_BEGIN
3669b54502bSHong Zhang #undef __FUNCT__
3672401956bSBarry Smith #define __FUNCT__ "PCFactorSetUseInPlace_Cholesky"
3682401956bSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetUseInPlace_Cholesky(PC pc)
3699b54502bSHong Zhang {
3709b54502bSHong Zhang   PC_Cholesky *dir;
3719b54502bSHong Zhang 
3729b54502bSHong Zhang   PetscFunctionBegin;
3739b54502bSHong Zhang   dir = (PC_Cholesky*)pc->data;
3749b54502bSHong Zhang   dir->inplace = PETSC_TRUE;
3759b54502bSHong Zhang   PetscFunctionReturn(0);
3769b54502bSHong Zhang }
3779b54502bSHong Zhang EXTERN_C_END
3789b54502bSHong Zhang 
3799b54502bSHong Zhang EXTERN_C_BEGIN
3809b54502bSHong Zhang #undef __FUNCT__
381e5a9bf91SBarry Smith #define __FUNCT__ "PCFactorSetMatOrderingType_Cholesky"
382e36ff184SSatish Balay PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatOrderingType_Cholesky(PC pc,const MatOrderingType ordering)
3839b54502bSHong Zhang {
3849b54502bSHong Zhang   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
3859b54502bSHong Zhang   PetscErrorCode ierr;
3869b54502bSHong Zhang 
3879b54502bSHong Zhang   PetscFunctionBegin;
3889b54502bSHong Zhang   ierr = PetscStrfree(dir->ordering);CHKERRQ(ierr);
3899b54502bSHong Zhang   ierr = PetscStrallocpy(ordering,&dir->ordering);CHKERRQ(ierr);
3909b54502bSHong Zhang   PetscFunctionReturn(0);
3919b54502bSHong Zhang }
3929b54502bSHong Zhang EXTERN_C_END
3939b54502bSHong Zhang 
3949b54502bSHong Zhang /* -----------------------------------------------------------------------------------*/
3959b54502bSHong Zhang 
3969b54502bSHong Zhang #undef __FUNCT__
3972401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseOrdering"
3989b54502bSHong Zhang /*@
3992401956bSBarry Smith    PCFactorSetReuseOrdering - When similar matrices are factored, this
4009b54502bSHong Zhang    causes the ordering computed in the first factor to be used for all
4019b54502bSHong Zhang    following factors.
4029b54502bSHong Zhang 
4039b54502bSHong Zhang    Collective on PC
4049b54502bSHong Zhang 
4059b54502bSHong Zhang    Input Parameters:
4069b54502bSHong Zhang +  pc - the preconditioner context
4079b54502bSHong Zhang -  flag - PETSC_TRUE to reuse else PETSC_FALSE
4089b54502bSHong Zhang 
4099b54502bSHong Zhang    Options Database Key:
4102401956bSBarry Smith .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
4119b54502bSHong Zhang 
4129b54502bSHong Zhang    Level: intermediate
4139b54502bSHong Zhang 
4149b54502bSHong Zhang .keywords: PC, levels, reordering, factorization, incomplete, LU
4159b54502bSHong Zhang 
4162401956bSBarry Smith .seealso: PCFactorSetReuseFill()
4179b54502bSHong Zhang @*/
4182401956bSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseOrdering(PC pc,PetscTruth flag)
4199b54502bSHong Zhang {
4209b54502bSHong Zhang   PetscErrorCode ierr,(*f)(PC,PetscTruth);
4219b54502bSHong Zhang 
4229b54502bSHong Zhang   PetscFunctionBegin;
4239b54502bSHong Zhang   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
4242401956bSBarry Smith   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);CHKERRQ(ierr);
4259b54502bSHong Zhang   if (f) {
4269b54502bSHong Zhang     ierr = (*f)(pc,flag);CHKERRQ(ierr);
4279b54502bSHong Zhang   }
4289b54502bSHong Zhang   PetscFunctionReturn(0);
4299b54502bSHong Zhang }
4309b54502bSHong Zhang 
4319b54502bSHong Zhang #undef __FUNCT__
4322401956bSBarry Smith #define __FUNCT__ "PCFactorSetReuseFill"
4339b54502bSHong Zhang /*@
4342401956bSBarry Smith    PCFactorSetReuseFill - When matrices with same different nonzero structure are factored,
4352401956bSBarry Smith    this causes later ones to use the fill ratio computed in the initial factorization.
4369b54502bSHong Zhang 
4379b54502bSHong Zhang    Collective on PC
4389b54502bSHong Zhang 
4399b54502bSHong Zhang    Input Parameters:
4409b54502bSHong Zhang +  pc - the preconditioner context
4419b54502bSHong Zhang -  flag - PETSC_TRUE to reuse else PETSC_FALSE
4429b54502bSHong Zhang 
4439b54502bSHong Zhang    Options Database Key:
4442401956bSBarry Smith .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
4459b54502bSHong Zhang 
4469b54502bSHong Zhang    Level: intermediate
4479b54502bSHong Zhang 
4489b54502bSHong Zhang .keywords: PC, levels, reordering, factorization, incomplete, Cholesky
4499b54502bSHong Zhang 
4502401956bSBarry Smith .seealso: PCFactorSetReuseOrdering()
4519b54502bSHong Zhang @*/
4522401956bSBarry Smith PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetReuseFill(PC pc,PetscTruth flag)
4539b54502bSHong Zhang {
4549b54502bSHong Zhang   PetscErrorCode ierr,(*f)(PC,PetscTruth);
4559b54502bSHong Zhang 
4569b54502bSHong Zhang   PetscFunctionBegin;
4579b54502bSHong Zhang   PetscValidHeaderSpecific(pc,PC_COOKIE,2);
4582401956bSBarry Smith   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseFill_C",(void (**)(void))&f);CHKERRQ(ierr);
4599b54502bSHong Zhang   if (f) {
4609b54502bSHong Zhang     ierr = (*f)(pc,flag);CHKERRQ(ierr);
4619b54502bSHong Zhang   }
4629b54502bSHong Zhang   PetscFunctionReturn(0);
4639b54502bSHong Zhang }
4649b54502bSHong Zhang 
4659b54502bSHong Zhang /*MC
46696fc60bcSBarry Smith    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
4679b54502bSHong Zhang 
4689b54502bSHong Zhang    Options Database Keys:
4692401956bSBarry Smith +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
470*c7393fdbSBarry Smith .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles
4712401956bSBarry Smith .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
47255ba2a51SBarry Smith .  -pc_factor_fill <fill> - Sets fill amount
4732401956bSBarry Smith .  -pc_factor_in_place - Activates in-place factorization
4742401956bSBarry Smith .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
475f251bdbdSHong Zhang .  -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
476f251bdbdSHong Zhang -  -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
477f251bdbdSHong Zhang    is optional with PETSC_TRUE being the default
4789b54502bSHong Zhang 
4799b54502bSHong Zhang    Notes: Not all options work for all matrix formats
4809b54502bSHong Zhang 
4819b54502bSHong Zhang    Level: beginner
4829b54502bSHong Zhang 
4839b54502bSHong Zhang    Concepts: Cholesky factorization, direct solver
4849b54502bSHong Zhang 
4859b54502bSHong Zhang    Notes: Usually this will compute an "exact" solution in one iteration and does
4869b54502bSHong Zhang           not need a Krylov method (i.e. you can use -ksp_type preonly, or
4879b54502bSHong Zhang           KSPSetType(ksp,KSPPREONLY) for the Krylov method
4889b54502bSHong Zhang 
4899b54502bSHong Zhang .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
490a4fd02acSBarry Smith            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
49155ba2a51SBarry Smith            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
492e5a9bf91SBarry Smith 	   PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd()
4939b54502bSHong Zhang 
4949b54502bSHong Zhang M*/
4959b54502bSHong Zhang 
4969b54502bSHong Zhang EXTERN_C_BEGIN
4979b54502bSHong Zhang #undef __FUNCT__
4989b54502bSHong Zhang #define __FUNCT__ "PCCreate_Cholesky"
499dba47a55SKris Buschelman PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Cholesky(PC pc)
5009b54502bSHong Zhang {
5019b54502bSHong Zhang   PetscErrorCode ierr;
5029b54502bSHong Zhang   PC_Cholesky    *dir;
5039b54502bSHong Zhang 
5049b54502bSHong Zhang   PetscFunctionBegin;
50538f2d2fdSLisandro Dalcin   ierr = PetscNewLog(pc,PC_Cholesky,&dir);CHKERRQ(ierr);
5069b54502bSHong Zhang 
5079b54502bSHong Zhang   dir->fact                   = 0;
5089b54502bSHong Zhang   dir->inplace                = PETSC_FALSE;
5099b54502bSHong Zhang   ierr = MatFactorInfoInitialize(&dir->info);CHKERRQ(ierr);
5109b54502bSHong Zhang   dir->info.fill              = 5.0;
5110a29876aSHong Zhang   dir->info.shiftnz           = 0.0;
512fbf22428SSatish Balay   dir->info.shiftpd           = 0.0; /* false */
5139b54502bSHong Zhang   dir->info.shift_fraction    = 0.0;
5149b54502bSHong Zhang   dir->info.pivotinblocks     = 1.0;
5159b54502bSHong Zhang   dir->col                    = 0;
5169b54502bSHong Zhang   dir->row                    = 0;
5179b54502bSHong Zhang   ierr = PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);CHKERRQ(ierr);
5185c9eb25fSBarry Smith   ierr = PetscStrallocpy("petsc",&dir->solvertype);CHKERRQ(ierr);
5199b54502bSHong Zhang   dir->reusefill        = PETSC_FALSE;
5209b54502bSHong Zhang   dir->reuseordering    = PETSC_FALSE;
5219b54502bSHong Zhang   pc->data              = (void*)dir;
5229b54502bSHong Zhang 
5239b54502bSHong Zhang   pc->ops->destroy           = PCDestroy_Cholesky;
5249b54502bSHong Zhang   pc->ops->apply             = PCApply_Cholesky;
5259b54502bSHong Zhang   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
5269b54502bSHong Zhang   pc->ops->setup             = PCSetUp_Cholesky;
5279b54502bSHong Zhang   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
5289b54502bSHong Zhang   pc->ops->view              = PCView_Cholesky;
5299b54502bSHong Zhang   pc->ops->applyrichardson   = 0;
530a4fd02acSBarry Smith   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Cholesky;
5319b54502bSHong Zhang 
532*c7393fdbSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Cholesky",
533*c7393fdbSBarry Smith                     PCFactorSetMatSolverPackage_Cholesky);CHKERRQ(ierr);
534afaefe49SHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Cholesky",
535afaefe49SHong Zhang                     PCFactorSetZeroPivot_Cholesky);CHKERRQ(ierr);
536afaefe49SHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Cholesky",
537afaefe49SHong Zhang                     PCFactorSetShiftNonzero_Cholesky);CHKERRQ(ierr);
538afaefe49SHong Zhang   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Cholesky",
539afaefe49SHong Zhang                     PCFactorSetShiftPd_Cholesky);CHKERRQ(ierr);
540afaefe49SHong Zhang 
54155ba2a51SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Cholesky",
54255ba2a51SBarry Smith                     PCFactorSetFill_Cholesky);CHKERRQ(ierr);
5432401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky",
5442401956bSBarry Smith                     PCFactorSetUseInPlace_Cholesky);CHKERRQ(ierr);
545e5a9bf91SBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Cholesky",
546e5a9bf91SBarry Smith                     PCFactorSetMatOrderingType_Cholesky);CHKERRQ(ierr);
5472401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky",
5482401956bSBarry Smith                     PCFactorSetReuseOrdering_Cholesky);CHKERRQ(ierr);
5492401956bSBarry Smith   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky",
5502401956bSBarry Smith                     PCFactorSetReuseFill_Cholesky);CHKERRQ(ierr);
5519b54502bSHong Zhang   PetscFunctionReturn(0);
5529b54502bSHong Zhang }
5539b54502bSHong Zhang EXTERN_C_END
554