Actual source code: qr.c
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>
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;
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;
21: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
22: if (dir->hdr.inplace) {
23: MatFactorType ftype;
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;
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: }
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: }
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: }
72: static PetscErrorCode PCReset_QR(PC pc)
73: {
74: PC_QR *dir = (PC_QR *)pc->data;
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: }
82: static PetscErrorCode PCDestroy_QR(PC pc)
83: {
84: PC_QR *dir = (PC_QR *)pc->data;
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: }
95: static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y)
96: {
97: PC_QR *dir = (PC_QR *)pc->data;
98: Mat fact;
100: PetscFunctionBegin;
101: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
102: PetscCall(MatSolve(fact, x, y));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y)
107: {
108: PC_QR *dir = (PC_QR *)pc->data;
109: Mat fact;
111: PetscFunctionBegin;
112: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
113: PetscCall(MatMatSolve(fact, X, Y));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y)
118: {
119: PC_QR *dir = (PC_QR *)pc->data;
120: Mat fact;
122: PetscFunctionBegin;
123: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
124: PetscCall(MatSolveTranspose(fact, x, y));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*MC
129: PCQR - Uses a direct solver, based on QR factorization, as a preconditioner
131: Level: beginner
133: Note:
134: Usually this will compute an "exact" solution in one iteration and does
135: not need a Krylov method (i.e. you can use -ksp_type preonly, or
136: `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
138: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSVD`,
139: `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
140: `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
141: `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
142: `PCFactorReorderForNonzeroDiagonal()`
143: M*/
145: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
146: {
147: PC_QR *dir;
149: PetscFunctionBegin;
150: PetscCall(PetscNew(&dir));
151: pc->data = (void *)dir;
152: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR));
154: dir->col = NULL;
155: pc->ops->reset = PCReset_QR;
156: pc->ops->destroy = PCDestroy_QR;
157: pc->ops->apply = PCApply_QR;
158: pc->ops->matapply = PCMatApply_QR;
159: pc->ops->applytranspose = PCApplyTranspose_QR;
160: pc->ops->setup = PCSetUp_QR;
161: pc->ops->setfromoptions = PCSetFromOptions_Factor;
162: pc->ops->view = PCView_Factor;
163: pc->ops->applyrichardson = NULL;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }