Actual source code: qr.c
2: /*
3: Defines a direct QR factorization preconditioner for any Mat implementation
4: Note: this need not be considered a preconditioner since it supplies
5: a direct solver.
6: */
7: #include <../src/ksp/pc/impls/factor/qr/qr.h>
9: static PetscErrorCode PCSetUp_QR(PC pc)
10: {
11: PC_QR *dir = (PC_QR *)pc->data;
12: MatSolverType stype;
13: MatFactorError err;
14: const char *prefix;
16: PetscFunctionBegin;
17: PetscCall(PCGetOptionsPrefix(pc, &prefix));
18: PetscCall(MatSetOptionsPrefix(pc->pmat, prefix));
19: pc->failedreason = PC_NOERROR;
20: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
22: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
23: if (dir->hdr.inplace) {
24: MatFactorType ftype;
26: PetscCall(MatGetFactorType(pc->pmat, &ftype));
27: if (ftype == MAT_FACTOR_NONE) {
28: PetscCall(MatQRFactor(pc->pmat, dir->col, &((PC_Factor *)dir)->info));
29: PetscCall(MatFactorGetError(pc->pmat, &err));
30: if (err) { /* Factor() fails */
31: pc->failedreason = (PCFailedReason)err;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: }
35: ((PC_Factor *)dir)->fact = pc->pmat;
36: } else {
37: MatInfo info;
39: if (!pc->setupcalled) {
40: if (!((PC_Factor *)dir)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact)); }
41: PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info));
42: PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
43: dir->hdr.actualfill = info.fill_ratio_needed;
44: } else if (pc->flag != SAME_NONZERO_PATTERN) {
45: PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info));
46: PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
47: dir->hdr.actualfill = info.fill_ratio_needed;
48: } else {
49: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
50: }
51: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
52: if (err) { /* FactorSymbolic() fails */
53: pc->failedreason = (PCFailedReason)err;
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: PetscCall(MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
58: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
59: if (err) { /* FactorNumeric() fails */
60: pc->failedreason = (PCFailedReason)err;
61: }
62: }
64: PetscCall(PCFactorGetMatSolverType(pc, &stype));
65: if (!stype) {
66: MatSolverType solverpackage;
67: PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
68: PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
69: }
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode PCReset_QR(PC pc)
74: {
75: PC_QR *dir = (PC_QR *)pc->data;
77: PetscFunctionBegin;
78: if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
79: PetscCall(ISDestroy(&dir->col));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode PCDestroy_QR(PC pc)
84: {
85: PC_QR *dir = (PC_QR *)pc->data;
87: PetscFunctionBegin;
88: PetscCall(PCReset_QR(pc));
89: PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
90: PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
91: PetscCall(PCFactorClearComposedFunctions(pc));
92: PetscCall(PetscFree(pc->data));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y)
97: {
98: PC_QR *dir = (PC_QR *)pc->data;
99: Mat fact;
101: PetscFunctionBegin;
102: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
103: PetscCall(MatSolve(fact, x, y));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y)
108: {
109: PC_QR *dir = (PC_QR *)pc->data;
110: Mat fact;
112: PetscFunctionBegin;
113: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
114: PetscCall(MatMatSolve(fact, X, Y));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y)
119: {
120: PC_QR *dir = (PC_QR *)pc->data;
121: Mat fact;
123: PetscFunctionBegin;
124: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
125: PetscCall(MatSolveTranspose(fact, x, y));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*MC
130: PCQR - Uses a direct solver, based on QR factorization, as a preconditioner
132: Level: beginner
134: Note:
135: Usually this will compute an "exact" solution in one iteration and does
136: not need a Krylov method (i.e. you can use -ksp_type preonly, or
137: `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
139: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSVD`,
140: `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
141: `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
142: `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
143: `PCFactorReorderForNonzeroDiagonal()`
144: M*/
146: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
147: {
148: PC_QR *dir;
150: PetscFunctionBegin;
151: PetscCall(PetscNew(&dir));
152: pc->data = (void *)dir;
153: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR));
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: }