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:   PetscFunctionBegin;
 10:   /* get work vectors needed by LCD */
 11:   PetscCall(KSPSetWorkVecs(ksp, 2));

 13:   PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->P));
 14:   PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q));
 15:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 40:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 41:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

 43:   lcd   = (KSP_LCD *)ksp->data;
 44:   X     = ksp->vec_sol;
 45:   B     = ksp->vec_rhs;
 46:   R     = ksp->work[0];
 47:   Z     = ksp->work[1];
 48:   max_k = lcd->restart;
 49:   mone  = -1;

 51:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

 53:   ksp->its = 0;
 54:   if (!ksp->guess_zero) {
 55:     PetscCall(KSP_MatMult(ksp, Amat, X, Z)); /*   z <- b - Ax       */
 56:     PetscCall(VecAYPX(Z, mone, B));
 57:   } else {
 58:     PetscCall(VecCopy(B, Z)); /*     z <- b (x is 0) */
 59:   }

 61:   PetscCall(KSP_PCApply(ksp, Z, R)); /*     r <- M^-1z         */
 62:   if (ksp->normtype != KSP_NORM_NONE) {
 63:     PetscCall(VecNorm(R, NORM_2, &rnorm));
 64:     KSPCheckNorm(ksp, rnorm);
 65:   }
 66:   PetscCall(KSPLogResidualHistory(ksp, rnorm));
 67:   PetscCall(KSPMonitor(ksp, 0, rnorm));
 68:   ksp->rnorm = rnorm;

 70:   /* test for convergence */
 71:   PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
 72:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 74:   PetscCall(VecCopy(R, lcd->P[0]));

 76:   while (!ksp->reason && ksp->its < ksp->max_it) {
 77:     it = 0;
 78:     PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it], Z));
 79:     PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it]));

 81:     while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
 82:       ksp->its++;
 83:       PetscCall(VecDot(lcd->P[it], R, &num));
 84:       PetscCall(VecDot(lcd->P[it], lcd->Q[it], &den));
 85:       KSPCheckDot(ksp, den);
 86:       alfa = num / den;
 87:       PetscCall(VecAXPY(X, alfa, lcd->P[it]));
 88:       PetscCall(VecAXPY(R, -alfa, lcd->Q[it]));
 89:       if (ksp->normtype != KSP_NORM_NONE) {
 90:         PetscCall(VecNorm(R, NORM_2, &rnorm));
 91:         KSPCheckNorm(ksp, rnorm);
 92:       }

 94:       ksp->rnorm = rnorm;
 95:       PetscCall(KSPLogResidualHistory(ksp, rnorm));
 96:       PetscCall(KSPMonitor(ksp, ksp->its, rnorm));
 97:       PetscCall((*ksp->converged)(ksp, ksp->its, rnorm, &ksp->reason, ksp->cnvP));

 99:       if (ksp->reason) break;

101:       PetscCall(VecCopy(R, lcd->P[it + 1]));
102:       PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it + 1], Z));
103:       PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it + 1]));

105:       for (j = 0; j <= it; j++) {
106:         PetscCall(VecDot(lcd->P[j], lcd->Q[it + 1], &num));
107:         KSPCheckDot(ksp, num);
108:         PetscCall(VecDot(lcd->P[j], lcd->Q[j], &den));
109:         beta = -num / den;
110:         PetscCall(VecAXPY(lcd->P[it + 1], beta, lcd->P[j]));
111:         PetscCall(VecAXPY(lcd->Q[it + 1], beta, lcd->Q[j]));
112:       }
113:       it++;
114:     }
115:     PetscCall(VecCopy(lcd->P[it], lcd->P[0]));
116:   }
117:   if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
118:   PetscCall(VecCopy(X, ksp->vec_sol));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }
121: /*
122:        KSPDestroy_LCD - Frees all memory space used by the Krylov method

124: */
125: PetscErrorCode KSPReset_LCD(KSP ksp)
126: {
127:   KSP_LCD *lcd = (KSP_LCD *)ksp->data;

129:   PetscFunctionBegin;
130:   if (lcd->P) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->P));
131:   if (lcd->Q) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->Q));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: PetscErrorCode KSPDestroy_LCD(KSP ksp)
136: {
137:   PetscFunctionBegin;
138:   PetscCall(KSPReset_LCD(ksp));
139:   PetscCall(PetscFree(ksp->data));
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: /*
144:      KSPView_LCD - Prints information about the current Krylov method being used

146:       Currently this only prints information to a file (or stdout) about the
147:       symmetry of the problem. If your Krylov method has special options or
148:       flags that information should be printed here.

150: */
151: PetscErrorCode KSPView_LCD(KSP ksp, PetscViewer viewer)
152: {
153:   KSP_LCD  *lcd = (KSP_LCD *)ksp->data;
154:   PetscBool iascii;

156:   PetscFunctionBegin;
157:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
158:   if (iascii) {
159:     PetscCall(PetscViewerASCIIPrintf(viewer, "  restart=%" PetscInt_FMT "\n", lcd->restart));
160:     PetscCall(PetscViewerASCIIPrintf(viewer, "  happy breakdown tolerance %g\n", (double)lcd->haptol));
161:   }
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

165: /*
166:     KSPSetFromOptions_LCD - Checks the options database for options related to the
167:                             LCD method.
168: */
169: PetscErrorCode KSPSetFromOptions_LCD(KSP ksp, PetscOptionItems *PetscOptionsObject)
170: {
171:   PetscBool flg;
172:   KSP_LCD  *lcd = (KSP_LCD *)ksp->data;

174:   PetscFunctionBegin;
175:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP LCD options");
176:   PetscCall(PetscOptionsInt("-ksp_lcd_restart", "Number of vectors conjugate", "KSPLCDSetRestart", lcd->restart, &lcd->restart, &flg));
177:   PetscCheck(!flg || lcd->restart >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Restart must be positive");
178:   PetscCall(PetscOptionsReal("-ksp_lcd_haptol", "Tolerance for exact convergence (happy ending)", "KSPLCDSetHapTol", lcd->haptol, &lcd->haptol, &flg));
179:   PetscCheck(!flg || lcd->haptol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*MC
184:      KSPLCD -  Implements the LCD (left conjugate direction) method

186:    Options Database Keys:
187: +   -ksp_lcd_restart - number of vectors conjugate
188: -   -ksp_lcd_haptol - tolerance for exact convergence (happy ending)

190:    Level: beginner

192:     Note:
193:     Support only for left preconditioning

195:     References:
196: +   * - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
197:      direction methods for real positive definite system. BIT Numerical
198:      Mathematics, 44(1),2004.
199: .   * - Y. Dai and J.Y. Yuan. Study on semiconjugate direction methods for
200:      nonsymmetric systems. International Journal for Numerical Methods in
201:      Engineering, 60, 2004.
202: .   * - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
203:      algorithm for solving linear systems of equations arising from implicit
204:      SUPG formulation of compressible flows. International Journal for
205:      Numerical Methods in Engineering, 60, 2004
206: -   * - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
207:      A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
208:      element and finite difference solution of convection diffusion
209:      equations,  Communications in Numerical Methods in Engineering, (Early
210:      View).

212:   Contributed by:
213:   Lucia Catabriga <luciac@ices.utexas.edu>

215: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
216:           `KSPCGSetType()`, `KSPLCDSetRestart()`, `KSPLCDSetHapTol()`
217: M*/

219: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
220: {
221:   KSP_LCD *lcd;

223:   PetscFunctionBegin;
224:   PetscCall(PetscNew(&lcd));
225:   ksp->data = (void *)lcd;
226:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
227:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
228:   lcd->restart = 30;
229:   lcd->haptol  = 1.0e-30;

231:   /*
232:        Sets the functions that are associated with this data structure
233:        (in C++ this is the same as defining virtual functions)
234:   */
235:   ksp->ops->setup          = KSPSetUp_LCD;
236:   ksp->ops->solve          = KSPSolve_LCD;
237:   ksp->ops->reset          = KSPReset_LCD;
238:   ksp->ops->destroy        = KSPDestroy_LCD;
239:   ksp->ops->view           = KSPView_LCD;
240:   ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
241:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
242:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }