Actual source code: lcd.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>
6: PetscErrorCode KSPSetUp_LCD(KSP ksp)
7: {
8: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
10: PetscInt restart = lcd->restart;
13: /* get work vectors needed by LCD */
14: KSPDefaultGetWork(ksp,2);
15:
16: VecDuplicateVecs(ksp->work[0],restart+1,&lcd->P);
17: VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q);
18: PetscLogObjectMemory(ksp,2*(restart+2)*sizeof(Vec));
19: return(0);
20: }
22: /* KSPSolve_LCD - This routine actually applies the left conjugate
23: direction method
25: Input Parameter:
26: . ksp - the Krylov space object that was set to use LCD, by, for
27: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPLCD);
29: Output Parameter:
30: . its - number of iterations used
32: */
35: PetscErrorCode KSPSolve_LCD(KSP ksp)
36: {
38: PetscInt it,j,max_k;
39: PetscScalar alfa, beta, num, den, mone;
40: PetscReal rnorm;
41: Vec X,B,R,Z;
42: KSP_LCD *lcd;
43: Mat Amat,Pmat;
44: MatStructure pflag;
45: PetscBool diagonalscale;
48: PCGetDiagonalScale(ksp->pc,&diagonalscale);
49: if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
51: lcd = (KSP_LCD*)ksp->data;
52: X = ksp->vec_sol;
53: B = ksp->vec_rhs;
54: R = ksp->work[0];
55: Z = ksp->work[1];
56: max_k = lcd->restart;
57: mone = -1;
59: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
61: ksp->its = 0;
62: if (!ksp->guess_zero) {
63: KSP_MatMult(ksp,Amat,X,Z); /* z <- b - Ax */
64: VecAYPX(Z,mone,B);
65: } else {
66: VecCopy(B,Z); /* z <- b (x is 0) */
67: }
68:
69: KSP_PCApply(ksp,Z,R); /* r <- M^-1z */
70: VecNorm(R,NORM_2,&rnorm);
71: KSPLogResidualHistory(ksp,rnorm);
72: KSPMonitor(ksp,0,rnorm);
73: ksp->rnorm = rnorm;
74:
75: /* test for convergence */
76: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
77: if (ksp->reason) return(0);
79: it = 0;
80: VecCopy(R,lcd->P[0]);
81:
82: while (!ksp->reason && ksp->its < ksp->max_it) {
83: it = 0;
84: KSP_MatMult(ksp,Amat,lcd->P[it],Z);
85: KSP_PCApply(ksp,Z,lcd->Q[it]);
86:
87: while(!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
88: ksp->its++;
89: VecDot(lcd->P[it],R,&num);
90: VecDot(lcd->P[it],lcd->Q[it], &den);
91: alfa = num/den;
92: VecAXPY(X,alfa,lcd->P[it]);
93: VecAXPY(R,-alfa,lcd->Q[it]);
94: VecNorm(R,NORM_2,&rnorm);
96: ksp->rnorm = rnorm;
97: KSPLogResidualHistory(ksp,rnorm);
98: KSPMonitor(ksp,ksp->its,rnorm);
99: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
100:
101: if (ksp->reason) break;
102:
103: VecCopy(R,lcd->P[it+1]);
104: KSP_MatMult(ksp,Amat,lcd->P[it+1],Z);
105: KSP_PCApply(ksp,Z,lcd->Q[it+1]);
106:
107: for( j = 0; j <= it; j++) {
108: VecDot(lcd->P[j],lcd->Q[it+1],&num);
109: VecDot(lcd->P[j],lcd->Q[j],&den);
110: beta = - num/den;
111: VecAXPY(lcd->P[it+1],beta,lcd->P[j]);
112: VecAXPY(lcd->Q[it+1],beta,lcd->Q[j]);
113: }
114: it++;
115: }
116: VecCopy(lcd->P[it],lcd->P[0]);
117: }
118: if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
119: VecCopy(X,ksp->vec_sol);
120:
121: return(0);
122: }
123: /*
124: KSPDestroy_LCD - Frees all memory space used by the Krylov method
126: */
129: PetscErrorCode KSPReset_LCD(KSP ksp)
130: {
131: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
135: if (lcd->P) { VecDestroyVecs(lcd->restart+1,&lcd->P);}
136: if (lcd->Q) { VecDestroyVecs(lcd->restart+1,&lcd->Q);}
137: return(0);
138: }
143: PetscErrorCode KSPDestroy_LCD(KSP ksp)
144: {
148: KSPReset_LCD(ksp);
149: PetscFree(ksp->data);
150: return(0);
151: }
153: /*
154: KSPView_LCD - Prints information about the current Krylov method being used
156: Currently this only prints information to a file (or stdout) about the
157: symmetry of the problem. If your Krylov method has special options or
158: flags that information should be printed here.
160: */
163: PetscErrorCode KSPView_LCD(KSP ksp,PetscViewer viewer)
164: {
166: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
168: PetscBool iascii;
171: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
172: if (iascii) {
173: PetscViewerASCIIPrintf(viewer," LCD: restart=%d\n",lcd->restart);
174: PetscViewerASCIIPrintf(viewer," LCD: happy breakdown tolerance %g\n",lcd->haptol);
175: } else {
176: SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for KSP LCD",((PetscObject)viewer)->type_name);
177: }
178: return(0);
179: }
181: /*
182: KSPSetFromOptions_LCD - Checks the options database for options related to the
183: LCD method.
184: */
187: PetscErrorCode KSPSetFromOptions_LCD(KSP ksp)
188: {
190: PetscBool flg;
191: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
192:
194: PetscOptionsHead("KSP LCD options");
195: PetscOptionsInt("-ksp_lcd_restart","Number of vectors conjugate","KSPLCDSetRestart",lcd->restart,&lcd->restart,&flg);
196: if(flg && lcd->restart < 1) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
197: PetscOptionsReal("-ksp_lcd_haptol","Tolerance for exact convergence (happy ending)","KSPLCDSetHapTol",lcd->haptol,&lcd->haptol,&flg);
198: if (flg && lcd->haptol < 0.0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
199: return(0);
200: }
202: /*MC
203: KSPLCD - Implements the LCD (left conjugate direction) method in PETSc.
205: Options Database Keys:
206: + -ksp_lcd_restart - number of vectors conjudate
207: - -ksp_lcd_haptol - tolerance for exact convergence (happing ending)
209: Level: beginner
211: Notes: Support only for left preconditioning
213: References:
214: - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
215: direction methods for real positive definite system. BIT Numerical
216: Mathematics, 44(1):189-207,2004.
217: - Y. Dai and J.Y. Yuan. Study on semi-conjugate direction methods for
218: non-symmetric systems. International Journal for Numerical Methods in
219: Engineering, 60:1383-1399,2004.
220: - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
221: algorithm for solving linear systems of equations arising from implicit
222: SUPG formulation of compressible flows. International Journal for
223: Numerical Methods in Engineering, 60:1513-1534,2004
224: - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
225: A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
226: element and finite difference solution of convection-diffusion
227: equations, Communications in Numerical Methods in Engineering, (Early
228: View).
230: Contributed by: Lucia Catabriga <luciac@ices.utexas.edu>
233: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
234: KSPCGSetType(), KSPLCDSetRestart(), KSPLCDSetHapTol()
236: M*/
238: EXTERN_C_BEGIN
241: PetscErrorCode KSPCreate_LCD(KSP ksp)
242: {
244: KSP_LCD *lcd;
247: PetscNewLog(ksp,KSP_LCD,&lcd);
248: ksp->data = (void*)lcd;
249: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
250: lcd->restart = 30;
251: lcd->haptol = 1.0e-30;
253: /*
254: Sets the functions that are associated with this data structure
255: (in C++ this is the same as defining virtual functions)
256: */
257: ksp->ops->setup = KSPSetUp_LCD;
258: ksp->ops->solve = KSPSolve_LCD;
259: ksp->ops->reset = KSPReset_LCD;
260: ksp->ops->destroy = KSPDestroy_LCD;
261: ksp->ops->view = KSPView_LCD;
262: ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
263: ksp->ops->buildsolution = KSPDefaultBuildSolution;
264: ksp->ops->buildresidual = KSPDefaultBuildResidual;
266: return(0);
267: }
268: EXTERN_C_END