Actual source code: lcd.c
2: #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>
4: PetscErrorCode KSPSetUp_LCD(KSP ksp)
5: {
6: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
7: PetscInt restart = lcd->restart;
9: PetscFunctionBegin;
10: /* get work vectors needed by LCD */
11: PetscCall(KSPSetWorkVecs(ksp, 2));
13: PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->P));
14: PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q));
15: PetscFunctionReturn(PETSC_SUCCESS);
16: }
18: /* KSPSolve_LCD - This routine actually applies the left conjugate
19: direction method
21: Input Parameter:
22: . ksp - the Krylov space object that was set to use LCD, by, for
23: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPLCD);
25: Output Parameter:
26: . its - number of iterations used
28: */
29: PetscErrorCode KSPSolve_LCD(KSP ksp)
30: {
31: PetscInt it, j, max_k;
32: PetscScalar alfa, beta, num, den, mone;
33: PetscReal rnorm = 0.0;
34: Vec X, B, R, Z;
35: KSP_LCD *lcd;
36: Mat Amat, Pmat;
37: PetscBool diagonalscale;
39: PetscFunctionBegin;
40: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
41: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
43: lcd = (KSP_LCD *)ksp->data;
44: X = ksp->vec_sol;
45: B = ksp->vec_rhs;
46: R = ksp->work[0];
47: Z = ksp->work[1];
48: max_k = lcd->restart;
49: mone = -1;
51: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
53: ksp->its = 0;
54: if (!ksp->guess_zero) {
55: PetscCall(KSP_MatMult(ksp, Amat, X, Z)); /* z <- b - Ax */
56: PetscCall(VecAYPX(Z, mone, B));
57: } else {
58: PetscCall(VecCopy(B, Z)); /* z <- b (x is 0) */
59: }
61: PetscCall(KSP_PCApply(ksp, Z, R)); /* r <- M^-1z */
62: if (ksp->normtype != KSP_NORM_NONE) {
63: PetscCall(VecNorm(R, NORM_2, &rnorm));
64: KSPCheckNorm(ksp, rnorm);
65: }
66: PetscCall(KSPLogResidualHistory(ksp, rnorm));
67: PetscCall(KSPMonitor(ksp, 0, rnorm));
68: ksp->rnorm = rnorm;
70: /* test for convergence */
71: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
72: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
74: PetscCall(VecCopy(R, lcd->P[0]));
76: while (!ksp->reason && ksp->its < ksp->max_it) {
77: it = 0;
78: PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it], Z));
79: PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it]));
81: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
82: ksp->its++;
83: PetscCall(VecDot(lcd->P[it], R, &num));
84: PetscCall(VecDot(lcd->P[it], lcd->Q[it], &den));
85: KSPCheckDot(ksp, den);
86: alfa = num / den;
87: PetscCall(VecAXPY(X, alfa, lcd->P[it]));
88: PetscCall(VecAXPY(R, -alfa, lcd->Q[it]));
89: if (ksp->normtype != KSP_NORM_NONE) {
90: PetscCall(VecNorm(R, NORM_2, &rnorm));
91: KSPCheckNorm(ksp, rnorm);
92: }
94: ksp->rnorm = rnorm;
95: PetscCall(KSPLogResidualHistory(ksp, rnorm));
96: PetscCall(KSPMonitor(ksp, ksp->its, rnorm));
97: PetscCall((*ksp->converged)(ksp, ksp->its, rnorm, &ksp->reason, ksp->cnvP));
99: if (ksp->reason) break;
101: PetscCall(VecCopy(R, lcd->P[it + 1]));
102: PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it + 1], Z));
103: PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it + 1]));
105: for (j = 0; j <= it; j++) {
106: PetscCall(VecDot(lcd->P[j], lcd->Q[it + 1], &num));
107: KSPCheckDot(ksp, num);
108: PetscCall(VecDot(lcd->P[j], lcd->Q[j], &den));
109: beta = -num / den;
110: PetscCall(VecAXPY(lcd->P[it + 1], beta, lcd->P[j]));
111: PetscCall(VecAXPY(lcd->Q[it + 1], beta, lcd->Q[j]));
112: }
113: it++;
114: }
115: PetscCall(VecCopy(lcd->P[it], lcd->P[0]));
116: }
117: if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
118: PetscCall(VecCopy(X, ksp->vec_sol));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
121: /*
122: KSPDestroy_LCD - Frees all memory space used by the Krylov method
124: */
125: PetscErrorCode KSPReset_LCD(KSP ksp)
126: {
127: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
129: PetscFunctionBegin;
130: if (lcd->P) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->P));
131: if (lcd->Q) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->Q));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: PetscErrorCode KSPDestroy_LCD(KSP ksp)
136: {
137: PetscFunctionBegin;
138: PetscCall(KSPReset_LCD(ksp));
139: PetscCall(PetscFree(ksp->data));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*
144: KSPView_LCD - Prints information about the current Krylov method being used
146: Currently this only prints information to a file (or stdout) about the
147: symmetry of the problem. If your Krylov method has special options or
148: flags that information should be printed here.
150: */
151: PetscErrorCode KSPView_LCD(KSP ksp, PetscViewer viewer)
152: {
153: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
154: PetscBool iascii;
156: PetscFunctionBegin;
157: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
158: if (iascii) {
159: PetscCall(PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT "\n", lcd->restart));
160: PetscCall(PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance %g\n", (double)lcd->haptol));
161: }
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*
166: KSPSetFromOptions_LCD - Checks the options database for options related to the
167: LCD method.
168: */
169: PetscErrorCode KSPSetFromOptions_LCD(KSP ksp, PetscOptionItems *PetscOptionsObject)
170: {
171: PetscBool flg;
172: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
174: PetscFunctionBegin;
175: PetscOptionsHeadBegin(PetscOptionsObject, "KSP LCD options");
176: PetscCall(PetscOptionsInt("-ksp_lcd_restart", "Number of vectors conjugate", "KSPLCDSetRestart", lcd->restart, &lcd->restart, &flg));
177: PetscCheck(!flg || lcd->restart >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Restart must be positive");
178: PetscCall(PetscOptionsReal("-ksp_lcd_haptol", "Tolerance for exact convergence (happy ending)", "KSPLCDSetHapTol", lcd->haptol, &lcd->haptol, &flg));
179: PetscCheck(!flg || lcd->haptol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*MC
184: KSPLCD - Implements the LCD (left conjugate direction) method
186: Options Database Keys:
187: + -ksp_lcd_restart - number of vectors conjugate
188: - -ksp_lcd_haptol - tolerance for exact convergence (happy ending)
190: Level: beginner
192: Note:
193: Support only for left preconditioning
195: References:
196: + * - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
197: direction methods for real positive definite system. BIT Numerical
198: Mathematics, 44(1),2004.
199: . * - Y. Dai and J.Y. Yuan. Study on semiconjugate direction methods for
200: nonsymmetric systems. International Journal for Numerical Methods in
201: Engineering, 60, 2004.
202: . * - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
203: algorithm for solving linear systems of equations arising from implicit
204: SUPG formulation of compressible flows. International Journal for
205: Numerical Methods in Engineering, 60, 2004
206: - * - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
207: A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
208: element and finite difference solution of convection diffusion
209: equations, Communications in Numerical Methods in Engineering, (Early
210: View).
212: Contributed by:
213: Lucia Catabriga <luciac@ices.utexas.edu>
215: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
216: `KSPCGSetType()`, `KSPLCDSetRestart()`, `KSPLCDSetHapTol()`
217: M*/
219: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
220: {
221: KSP_LCD *lcd;
223: PetscFunctionBegin;
224: PetscCall(PetscNew(&lcd));
225: ksp->data = (void *)lcd;
226: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
227: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
228: lcd->restart = 30;
229: lcd->haptol = 1.0e-30;
231: /*
232: Sets the functions that are associated with this data structure
233: (in C++ this is the same as defining virtual functions)
234: */
235: ksp->ops->setup = KSPSetUp_LCD;
236: ksp->ops->solve = KSPSolve_LCD;
237: ksp->ops->reset = KSPReset_LCD;
238: ksp->ops->destroy = KSPDestroy_LCD;
239: ksp->ops->view = KSPView_LCD;
240: ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
241: ksp->ops->buildsolution = KSPBuildSolutionDefault;
242: ksp->ops->buildresidual = KSPBuildResidualDefault;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }