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