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: PetscLogObjectMemory((PetscObject)ksp,2*(restart+2)*sizeof(Vec));
15: return 0;
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: PCGetDiagonalScale(ksp->pc,&diagonalscale);
42: lcd = (KSP_LCD*)ksp->data;
43: X = ksp->vec_sol;
44: B = ksp->vec_rhs;
45: R = ksp->work[0];
46: Z = ksp->work[1];
47: max_k = lcd->restart;
48: mone = -1;
50: PCGetOperators(ksp->pc,&Amat,&Pmat);
52: ksp->its = 0;
53: if (!ksp->guess_zero) {
54: KSP_MatMult(ksp,Amat,X,Z); /* z <- b - Ax */
55: VecAYPX(Z,mone,B);
56: } else {
57: VecCopy(B,Z); /* z <- b (x is 0) */
58: }
60: KSP_PCApply(ksp,Z,R); /* r <- M^-1z */
61: if (ksp->normtype != KSP_NORM_NONE) {
62: VecNorm(R,NORM_2,&rnorm);
63: KSPCheckNorm(ksp,rnorm);
64: }
65: KSPLogResidualHistory(ksp,rnorm);
66: KSPMonitor(ksp,0,rnorm);
67: ksp->rnorm = rnorm;
69: /* test for convergence */
70: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
71: if (ksp->reason) return 0;
73: VecCopy(R,lcd->P[0]);
75: while (!ksp->reason && ksp->its < ksp->max_it) {
76: it = 0;
77: KSP_MatMult(ksp,Amat,lcd->P[it],Z);
78: KSP_PCApply(ksp,Z,lcd->Q[it]);
80: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
81: ksp->its++;
82: VecDot(lcd->P[it],R,&num);
83: VecDot(lcd->P[it],lcd->Q[it], &den);
84: KSPCheckDot(ksp,den);
85: alfa = num/den;
86: VecAXPY(X,alfa,lcd->P[it]);
87: VecAXPY(R,-alfa,lcd->Q[it]);
88: if (ksp->normtype != KSP_NORM_NONE) {
89: VecNorm(R,NORM_2,&rnorm);
90: KSPCheckNorm(ksp,rnorm);
91: }
93: ksp->rnorm = rnorm;
94: KSPLogResidualHistory(ksp,rnorm);
95: KSPMonitor(ksp,ksp->its,rnorm);
96: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
98: if (ksp->reason) break;
100: VecCopy(R,lcd->P[it+1]);
101: KSP_MatMult(ksp,Amat,lcd->P[it+1],Z);
102: KSP_PCApply(ksp,Z,lcd->Q[it+1]);
104: for (j = 0; j <= it; j++) {
105: VecDot(lcd->P[j],lcd->Q[it+1],&num);
106: KSPCheckDot(ksp,num);
107: VecDot(lcd->P[j],lcd->Q[j],&den);
108: beta = -num/den;
109: VecAXPY(lcd->P[it+1],beta,lcd->P[j]);
110: VecAXPY(lcd->Q[it+1],beta,lcd->Q[j]);
111: }
112: it++;
113: }
114: VecCopy(lcd->P[it],lcd->P[0]);
115: }
116: if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
117: VecCopy(X,ksp->vec_sol);
118: return 0;
119: }
120: /*
121: KSPDestroy_LCD - Frees all memory space used by the Krylov method
123: */
124: PetscErrorCode KSPReset_LCD(KSP ksp)
125: {
126: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
128: if (lcd->P) VecDestroyVecs(lcd->restart+1,&lcd->P);
129: if (lcd->Q) VecDestroyVecs(lcd->restart+1,&lcd->Q);
130: return 0;
131: }
133: PetscErrorCode KSPDestroy_LCD(KSP ksp)
134: {
135: KSPReset_LCD(ksp);
136: PetscFree(ksp->data);
137: return 0;
138: }
140: /*
141: KSPView_LCD - Prints information about the current Krylov method being used
143: Currently this only prints information to a file (or stdout) about the
144: symmetry of the problem. If your Krylov method has special options or
145: flags that information should be printed here.
147: */
148: PetscErrorCode KSPView_LCD(KSP ksp,PetscViewer viewer)
149: {
151: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
152: PetscBool iascii;
154: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
155: if (iascii) {
156: PetscViewerASCIIPrintf(viewer," restart=%d\n",lcd->restart);
157: PetscViewerASCIIPrintf(viewer," happy breakdown tolerance %g\n",lcd->haptol);
158: }
159: return 0;
160: }
162: /*
163: KSPSetFromOptions_LCD - Checks the options database for options related to the
164: LCD method.
165: */
166: PetscErrorCode KSPSetFromOptions_LCD(PetscOptionItems *PetscOptionsObject,KSP ksp)
167: {
168: PetscBool flg;
169: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
171: PetscOptionsHead(PetscOptionsObject,"KSP LCD options");
172: PetscOptionsInt("-ksp_lcd_restart","Number of vectors conjugate","KSPLCDSetRestart",lcd->restart,&lcd->restart,&flg);
174: PetscOptionsReal("-ksp_lcd_haptol","Tolerance for exact convergence (happy ending)","KSPLCDSetHapTol",lcd->haptol,&lcd->haptol,&flg);
176: return 0;
177: }
179: /*MC
180: KSPLCD - Implements the LCD (left conjugate direction) method in PETSc.
182: Options Database Keys:
183: + -ksp_lcd_restart - number of vectors conjudate
184: - -ksp_lcd_haptol - tolerance for exact convergence (happing ending)
186: Level: beginner
188: Notes:
189: Support only for left preconditioning
191: References:
192: + * - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
193: direction methods for real positive definite system. BIT Numerical
194: Mathematics, 44(1),2004.
195: . * - Y. Dai and J.Y. Yuan. Study on semiconjugate direction methods for
196: nonsymmetric systems. International Journal for Numerical Methods in
197: Engineering, 60, 2004.
198: . * - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
199: algorithm for solving linear systems of equations arising from implicit
200: SUPG formulation of compressible flows. International Journal for
201: Numerical Methods in Engineering, 60, 2004
202: - * - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
203: A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
204: element and finite difference solution of convection diffusion
205: equations, Communications in Numerical Methods in Engineering, (Early
206: View).
208: Contributed by: Lucia Catabriga <luciac@ices.utexas.edu>
210: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
211: KSPCGSetType(), KSPLCDSetRestart(), KSPLCDSetHapTol()
213: M*/
215: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
216: {
217: KSP_LCD *lcd;
219: PetscNewLog(ksp,&lcd);
220: ksp->data = (void*)lcd;
221: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
222: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
223: lcd->restart = 30;
224: lcd->haptol = 1.0e-30;
226: /*
227: Sets the functions that are associated with this data structure
228: (in C++ this is the same as defining virtual functions)
229: */
230: ksp->ops->setup = KSPSetUp_LCD;
231: ksp->ops->solve = KSPSolve_LCD;
232: ksp->ops->reset = KSPReset_LCD;
233: ksp->ops->destroy = KSPDestroy_LCD;
234: ksp->ops->view = KSPView_LCD;
235: ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
236: ksp->ops->buildsolution = KSPBuildSolutionDefault;
237: ksp->ops->buildresidual = KSPBuildResidualDefault;
238: return 0;
239: }