1 /*
2 Defines a direct QR factorization preconditioner for any Mat implementation
3 Note: this need not be considered a preconditioner since it supplies
4 a direct solver.
5 */
6 #include <../src/ksp/pc/impls/factor/qr/qr.h> /*I "petscpc.h" I*/
7
PCSetUp_QR(PC pc)8 static PetscErrorCode PCSetUp_QR(PC pc)
9 {
10 PC_QR *dir = (PC_QR *)pc->data;
11 MatSolverType stype;
12 MatFactorError err;
13 const char *prefix;
14
15 PetscFunctionBegin;
16 PetscCall(PCGetOptionsPrefix(pc, &prefix));
17 PetscCall(MatSetOptionsPrefix(pc->pmat, prefix));
18 pc->failedreason = PC_NOERROR;
19 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
20
21 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
22 if (dir->hdr.inplace) {
23 MatFactorType ftype;
24
25 PetscCall(MatGetFactorType(pc->pmat, &ftype));
26 if (ftype == MAT_FACTOR_NONE) {
27 PetscCall(MatQRFactor(pc->pmat, dir->col, &((PC_Factor *)dir)->info));
28 PetscCall(MatFactorGetError(pc->pmat, &err));
29 if (err) { /* Factor() fails */
30 pc->failedreason = (PCFailedReason)err;
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33 }
34 ((PC_Factor *)dir)->fact = pc->pmat;
35 } else {
36 MatInfo info;
37
38 if (!pc->setupcalled) {
39 if (!((PC_Factor *)dir)->fact) PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact));
40 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info));
41 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
42 dir->hdr.actualfill = info.fill_ratio_needed;
43 } else if (pc->flag != SAME_NONZERO_PATTERN) {
44 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info));
45 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
46 dir->hdr.actualfill = info.fill_ratio_needed;
47 } else {
48 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
49 }
50 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
51 if (err) { /* FactorSymbolic() fails */
52 pc->failedreason = (PCFailedReason)err;
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55
56 PetscCall(MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
57 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
58 if (err) { /* FactorNumeric() fails */
59 pc->failedreason = (PCFailedReason)err;
60 }
61 }
62
63 PetscCall(PCFactorGetMatSolverType(pc, &stype));
64 if (!stype) {
65 MatSolverType solverpackage;
66 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
67 PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
68 }
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71
PCReset_QR(PC pc)72 static PetscErrorCode PCReset_QR(PC pc)
73 {
74 PC_QR *dir = (PC_QR *)pc->data;
75
76 PetscFunctionBegin;
77 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
78 PetscCall(ISDestroy(&dir->col));
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81
PCDestroy_QR(PC pc)82 static PetscErrorCode PCDestroy_QR(PC pc)
83 {
84 PC_QR *dir = (PC_QR *)pc->data;
85
86 PetscFunctionBegin;
87 PetscCall(PCReset_QR(pc));
88 PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
89 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
90 PetscCall(PCFactorClearComposedFunctions(pc));
91 PetscCall(PetscFree(pc->data));
92 PetscFunctionReturn(PETSC_SUCCESS);
93 }
94
PCApply_QR(PC pc,Vec x,Vec y)95 static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y)
96 {
97 PC_QR *dir = (PC_QR *)pc->data;
98 Mat fact;
99
100 PetscFunctionBegin;
101 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
102 PetscCall(MatSolve(fact, x, y));
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
PCMatApply_QR(PC pc,Mat X,Mat Y)106 static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y)
107 {
108 PC_QR *dir = (PC_QR *)pc->data;
109 Mat fact;
110
111 PetscFunctionBegin;
112 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
113 PetscCall(MatMatSolve(fact, X, Y));
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
PCApplyTranspose_QR(PC pc,Vec x,Vec y)117 static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y)
118 {
119 PC_QR *dir = (PC_QR *)pc->data;
120 Mat fact;
121
122 PetscFunctionBegin;
123 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
124 PetscCall(MatSolveTranspose(fact, x, y));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /*MC
129 PCQR - Uses a direct solver, based on QR factorization, as a preconditioner
130
131 Level: beginner
132
133 Notes:
134 This is usually used as a preconditioner for `KSPLSQR`
135
136 If the matrix used to construct the preconditioner is a `MATNORMAL` or `MATNORMALHERMITIAN` matrix, then the QR is performed on the
137 inner matrix to provide the preconditioner.
138
139 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSVD`,
140 `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
141 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
142 `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
143 `PCFactorReorderForNonzeroDiagonal()`, `KSPLSQR`
144 M*/
145
PCCreate_QR(PC pc)146 PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
147 {
148 PC_QR *dir;
149
150 PetscFunctionBegin;
151 PetscCall(PetscNew(&dir));
152 pc->data = (void *)dir;
153 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR));
154
155 dir->col = NULL;
156 pc->ops->reset = PCReset_QR;
157 pc->ops->destroy = PCDestroy_QR;
158 pc->ops->apply = PCApply_QR;
159 pc->ops->matapply = PCMatApply_QR;
160 pc->ops->applytranspose = PCApplyTranspose_QR;
161 pc->ops->setup = PCSetUp_QR;
162 pc->ops->setfromoptions = PCSetFromOptions_Factor;
163 pc->ops->view = PCView_Factor;
164 pc->ops->applyrichardson = NULL;
165 PetscFunctionReturn(PETSC_SUCCESS);
166 }
167