Actual source code: lcd.c
petsc-3.7.5 2017-01-01
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: KSPSetWorkVecs(ksp,2);
16: VecDuplicateVecs(ksp->work[0],restart+1,&lcd->P);
17: VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q);
18: PetscLogObjectMemory((PetscObject)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: PetscBool diagonalscale;
47: PCGetDiagonalScale(ksp->pc,&diagonalscale);
48: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
50: lcd = (KSP_LCD*)ksp->data;
51: X = ksp->vec_sol;
52: B = ksp->vec_rhs;
53: R = ksp->work[0];
54: Z = ksp->work[1];
55: max_k = lcd->restart;
56: mone = -1;
58: PCGetOperators(ksp->pc,&Amat,&Pmat);
60: ksp->its = 0;
61: if (!ksp->guess_zero) {
62: KSP_MatMult(ksp,Amat,X,Z); /* z <- b - Ax */
63: VecAYPX(Z,mone,B);
64: } else {
65: VecCopy(B,Z); /* z <- b (x is 0) */
66: }
68: KSP_PCApply(ksp,Z,R); /* r <- M^-1z */
69: VecNorm(R,NORM_2,&rnorm);
70: KSPLogResidualHistory(ksp,rnorm);
71: KSPMonitor(ksp,0,rnorm);
72: ksp->rnorm = rnorm;
74: /* test for convergence */
75: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
76: if (ksp->reason) return(0);
78: VecCopy(R,lcd->P[0]);
80: while (!ksp->reason && ksp->its < ksp->max_it) {
81: it = 0;
82: KSP_MatMult(ksp,Amat,lcd->P[it],Z);
83: KSP_PCApply(ksp,Z,lcd->Q[it]);
85: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
86: ksp->its++;
87: VecDot(lcd->P[it],R,&num);
88: VecDot(lcd->P[it],lcd->Q[it], &den);
89: alfa = num/den;
90: VecAXPY(X,alfa,lcd->P[it]);
91: VecAXPY(R,-alfa,lcd->Q[it]);
92: VecNorm(R,NORM_2,&rnorm);
94: ksp->rnorm = rnorm;
95: KSPLogResidualHistory(ksp,rnorm);
96: KSPMonitor(ksp,ksp->its,rnorm);
97: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
99: if (ksp->reason) break;
101: VecCopy(R,lcd->P[it+1]);
102: KSP_MatMult(ksp,Amat,lcd->P[it+1],Z);
103: KSP_PCApply(ksp,Z,lcd->Q[it+1]);
105: for (j = 0; j <= it; j++) {
106: VecDot(lcd->P[j],lcd->Q[it+1],&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: */
126: PetscErrorCode KSPReset_LCD(KSP ksp)
127: {
128: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
132: if (lcd->P) { VecDestroyVecs(lcd->restart+1,&lcd->P);}
133: if (lcd->Q) { VecDestroyVecs(lcd->restart+1,&lcd->Q);}
134: return(0);
135: }
140: PetscErrorCode KSPDestroy_LCD(KSP ksp)
141: {
145: KSPReset_LCD(ksp);
146: PetscFree(ksp->data);
147: return(0);
148: }
150: /*
151: KSPView_LCD - Prints information about the current Krylov method being used
153: Currently this only prints information to a file (or stdout) about the
154: symmetry of the problem. If your Krylov method has special options or
155: flags that information should be printed here.
157: */
160: PetscErrorCode KSPView_LCD(KSP ksp,PetscViewer viewer)
161: {
163: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
165: PetscBool iascii;
168: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
169: if (iascii) {
170: PetscViewerASCIIPrintf(viewer," LCD: restart=%d\n",lcd->restart);
171: PetscViewerASCIIPrintf(viewer," LCD: happy breakdown tolerance %g\n",lcd->haptol);
172: }
173: return(0);
174: }
176: /*
177: KSPSetFromOptions_LCD - Checks the options database for options related to the
178: LCD method.
179: */
182: PetscErrorCode KSPSetFromOptions_LCD(PetscOptionItems *PetscOptionsObject,KSP ksp)
183: {
185: PetscBool flg;
186: KSP_LCD *lcd = (KSP_LCD*)ksp->data;
189: PetscOptionsHead(PetscOptionsObject,"KSP LCD options");
190: PetscOptionsInt("-ksp_lcd_restart","Number of vectors conjugate","KSPLCDSetRestart",lcd->restart,&lcd->restart,&flg);
191: if (flg && lcd->restart < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
192: PetscOptionsReal("-ksp_lcd_haptol","Tolerance for exact convergence (happy ending)","KSPLCDSetHapTol",lcd->haptol,&lcd->haptol,&flg);
193: if (flg && lcd->haptol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
194: return(0);
195: }
197: /*MC
198: KSPLCD - Implements the LCD (left conjugate direction) method in PETSc.
200: Options Database Keys:
201: + -ksp_lcd_restart - number of vectors conjudate
202: - -ksp_lcd_haptol - tolerance for exact convergence (happing ending)
204: Level: beginner
206: Notes: Support only for left preconditioning
208: References:
209: + 1. - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
210: direction methods for real positive definite system. BIT Numerical
211: Mathematics, 44(1),2004.
212: . 2. - Y. Dai and J.Y. Yuan. Study on semiconjugate direction methods for
213: nonsymmetric systems. International Journal for Numerical Methods in
214: Engineering, 60, 2004.
215: . 3. - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
216: algorithm for solving linear systems of equations arising from implicit
217: SUPG formulation of compressible flows. International Journal for
218: Numerical Methods in Engineering, 60, 2004
219: - 4. - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
220: A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
221: element and finite difference solution of convection diffusion
222: equations, Communications in Numerical Methods in Engineering, (Early
223: View).
225: Contributed by: Lucia Catabriga <luciac@ices.utexas.edu>
228: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
229: KSPCGSetType(), KSPLCDSetRestart(), KSPLCDSetHapTol()
231: M*/
235: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
236: {
238: KSP_LCD *lcd;
241: PetscNewLog(ksp,&lcd);
242: ksp->data = (void*)lcd;
243: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
244: lcd->restart = 30;
245: lcd->haptol = 1.0e-30;
247: /*
248: Sets the functions that are associated with this data structure
249: (in C++ this is the same as defining virtual functions)
250: */
251: ksp->ops->setup = KSPSetUp_LCD;
252: ksp->ops->solve = KSPSolve_LCD;
253: ksp->ops->reset = KSPReset_LCD;
254: ksp->ops->destroy = KSPDestroy_LCD;
255: ksp->ops->view = KSPView_LCD;
256: ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
257: ksp->ops->buildsolution = KSPBuildSolutionDefault;
258: ksp->ops->buildresidual = KSPBuildResidualDefault;
259: return(0);
260: }