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;
15: pc->failedreason = PC_NOERROR;
16: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
18: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
19: if (dir->hdr.inplace) {
20: MatFactorType ftype;
22: MatGetFactorType(pc->pmat, &ftype);
23: if (ftype == MAT_FACTOR_NONE) {
24: MatQRFactor(pc->pmat,dir->col,&((PC_Factor*)dir)->info);
25: MatFactorGetError(pc->pmat,&err);
26: if (err) { /* Factor() fails */
27: pc->failedreason = (PCFailedReason)err;
28: return 0;
29: }
30: }
31: ((PC_Factor*)dir)->fact = pc->pmat;
32: } else {
33: MatInfo info;
35: if (!pc->setupcalled) {
36: if (!((PC_Factor*)dir)->fact) {
37: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_QR,&((PC_Factor*)dir)->fact);
38: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
39: }
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: PetscFree(pc->data);
89: return 0;
90: }
92: static PetscErrorCode PCApply_QR(PC pc,Vec x,Vec y)
93: {
94: PC_QR *dir = (PC_QR*)pc->data;
95: Mat fact;
97: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
98: MatSolve(fact,x,y);
99: return 0;
100: }
102: static PetscErrorCode PCMatApply_QR(PC pc,Mat X,Mat Y)
103: {
104: PC_QR *dir = (PC_QR*)pc->data;
105: Mat fact;
107: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
108: MatMatSolve(fact,X,Y);
109: return 0;
110: }
112: static PetscErrorCode PCApplyTranspose_QR(PC pc,Vec x,Vec y)
113: {
114: PC_QR *dir = (PC_QR*)pc->data;
115: Mat fact;
117: fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
118: MatSolveTranspose(fact,x,y);
119: return 0;
120: }
122: /* -----------------------------------------------------------------------------------*/
124: /*MC
125: PCQR - Uses a direct solver, based on QR factorization, as a preconditioner
127: Level: beginner
129: Notes:
130: Usually this will compute an "exact" solution in one iteration and does
131: not need a Krylov method (i.e. you can use -ksp_type preonly, or
132: KSPSetType(ksp,KSPPREONLY) for the Krylov method
134: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
135: PCILU, PCLU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
136: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
137: PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
138: PCFactorReorderForNonzeroDiagonal()
139: M*/
141: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
142: {
143: PC_QR *dir;
145: PetscNewLog(pc,&dir);
146: pc->data = (void*)dir;
147: PCFactorInitialize(pc, MAT_FACTOR_QR);
149: dir->col = NULL;
150: pc->ops->reset = PCReset_QR;
151: pc->ops->destroy = PCDestroy_QR;
152: pc->ops->apply = PCApply_QR;
153: pc->ops->matapply = PCMatApply_QR;
154: pc->ops->applytranspose = PCApplyTranspose_QR;
155: pc->ops->setup = PCSetUp_QR;
156: pc->ops->setfromoptions = PCSetFromOptions_Factor;
157: pc->ops->view = PCView_Factor;
158: pc->ops->applyrichardson = NULL;
159: return 0;
160: }