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: }