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: PCGetOptionsPrefix(pc, &prefix);
17: 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: MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
22: if (dir->hdr.inplace) {
23: MatFactorType ftype;
25: MatGetFactorType(pc->pmat, &ftype);
26: if (ftype == MAT_FACTOR_NONE) {
27: MatQRFactor(pc->pmat, dir->col, &((PC_Factor *)dir)->info);
28: MatFactorGetError(pc->pmat, &err);
29: if (err) { /* Factor() fails */
30: pc->failedreason = (PCFailedReason)err;
31: return 0;
32: }
33: }
34: ((PC_Factor *)dir)->fact = pc->pmat;
35: } else {
36: MatInfo info;
38: if (!pc->setupcalled) {
39: if (!((PC_Factor *)dir)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact); }
40: MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info);
41: 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: MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info);
45: MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
46: dir->hdr.actualfill = info.fill_ratio_needed;
47: } else {
48: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
49: }
50: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
51: if (err) { /* FactorSymbolic() fails */
52: pc->failedreason = (PCFailedReason)err;
53: return 0;
54: }
56: MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info);
57: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
58: if (err) { /* FactorNumeric() fails */
59: pc->failedreason = (PCFailedReason)err;
60: }
61: }
63: PCFactorGetMatSolverType(pc, &stype);
64: if (!stype) {
65: MatSolverType solverpackage;
66: MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage);
67: PCFactorSetMatSolverType(pc, solverpackage);
68: }
69: return 0;
70: }
72: static PetscErrorCode PCReset_QR(PC pc)
73: {
74: PC_QR *dir = (PC_QR *)pc->data;
76: if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) MatDestroy(&((PC_Factor *)dir)->fact);
77: ISDestroy(&dir->col);
78: return 0;
79: }
81: static PetscErrorCode PCDestroy_QR(PC pc)
82: {
83: PC_QR *dir = (PC_QR *)pc->data;
85: PCReset_QR(pc);
86: PetscFree(((PC_Factor *)dir)->ordering);
87: PetscFree(((PC_Factor *)dir)->solvertype);
88: PCFactorClearComposedFunctions(pc);
89: PetscFree(pc->data);
90: return 0;
91: }
93: static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y)
94: {
95: PC_QR *dir = (PC_QR *)pc->data;
96: Mat fact;
98: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
99: MatSolve(fact, x, y);
100: return 0;
101: }
103: static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y)
104: {
105: PC_QR *dir = (PC_QR *)pc->data;
106: Mat fact;
108: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
109: MatMatSolve(fact, X, Y);
110: return 0;
111: }
113: static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y)
114: {
115: PC_QR *dir = (PC_QR *)pc->data;
116: Mat fact;
118: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
119: MatSolveTranspose(fact, x, y);
120: return 0;
121: }
123: /*MC
124: PCQR - Uses a direct solver, based on QR factorization, as a preconditioner
126: Level: beginner
128: Note:
129: Usually this will compute an "exact" solution in one iteration and does
130: not need a Krylov method (i.e. you can use -ksp_type preonly, or
131: `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
133: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSVD`,
134: `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
135: `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
136: `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
137: `PCFactorReorderForNonzeroDiagonal()`
138: M*/
140: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
141: {
142: PC_QR *dir;
144: PetscNew(&dir);
145: pc->data = (void *)dir;
146: PCFactorInitialize(pc, MAT_FACTOR_QR);
148: dir->col = NULL;
149: pc->ops->reset = PCReset_QR;
150: pc->ops->destroy = PCDestroy_QR;
151: pc->ops->apply = PCApply_QR;
152: pc->ops->matapply = PCMatApply_QR;
153: pc->ops->applytranspose = PCApplyTranspose_QR;
154: pc->ops->setup = PCSetUp_QR;
155: pc->ops->setfromoptions = PCSetFromOptions_Factor;
156: pc->ops->view = PCView_Factor;
157: pc->ops->applyrichardson = NULL;
158: return 0;
159: }