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