Actual source code: lcd.c

  1: #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>

  3: static PetscErrorCode KSPSetUp_LCD(KSP ksp)
  4: {
  5:   KSP_LCD *lcd     = (KSP_LCD *)ksp->data;
  6:   PetscInt restart = lcd->restart;

  8:   PetscFunctionBegin;
  9:   /* get work vectors needed by LCD */
 10:   PetscCall(KSPSetWorkVecs(ksp, 2));

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

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

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

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

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

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

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

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

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

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

 98:       if (ksp->reason) break;

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

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

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

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

134: static PetscErrorCode KSPDestroy_LCD(KSP ksp)
135: {
136:   PetscFunctionBegin;
137:   PetscCall(KSPReset_LCD(ksp));
138:   PetscCall(PetscFree(ksp->data));
139:   PetscFunctionReturn(PETSC_SUCCESS);
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: static PetscErrorCode KSPView_LCD(KSP ksp, PetscViewer viewer)
151: {
152:   KSP_LCD  *lcd = (KSP_LCD *)ksp->data;
153:   PetscBool iascii;

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

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

173:   PetscFunctionBegin;
174:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP LCD options");
175:   PetscCall(PetscOptionsBoundedInt("-ksp_lcd_restart", "Number of vectors conjugate", "KSPLCDSetRestart", lcd->restart, &lcd->restart, &flg, 1));
176:   PetscCall(PetscOptionsBoundedReal("-ksp_lcd_haptol", "Tolerance for exact convergence (happy ending)", "KSPLCDSetHapTol", lcd->haptol, &lcd->haptol, &flg, 0.0));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*MC
181:    KSPLCD -  Implements the LCD (left conjugate direction) method

183:    Options Database Keys:
184: +  -ksp_lcd_restart - number of vectors conjugate
185: -  -ksp_lcd_haptol - tolerance for exact convergence (happy ending)

187:    Level: beginner

189:    Notes:
190:    Support only for left preconditioning

192:    See {cite}`yuan2004semi`, {cite}`dai2004study`, {cite}`catabriga2004evaluating`, and {cite}`catabriga2006performance`

194:    Contributed by:
195:    Lucia Catabriga <luciac@ices.utexas.edu>

197: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`,
198:           `KSPCGSetType()`, `KSPLCDSetRestart()`, `KSPLCDSetHapTol()`
199: M*/

201: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
202: {
203:   KSP_LCD *lcd;

205:   PetscFunctionBegin;
206:   PetscCall(PetscNew(&lcd));
207:   ksp->data = (void *)lcd;
208:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
209:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
210:   lcd->restart = 30;
211:   lcd->haptol  = 1.0e-30;

213:   /*
214:        Sets the functions that are associated with this data structure
215:        (in C++ this is the same as defining virtual functions)
216:   */
217:   ksp->ops->setup          = KSPSetUp_LCD;
218:   ksp->ops->solve          = KSPSolve_LCD;
219:   ksp->ops->reset          = KSPReset_LCD;
220:   ksp->ops->destroy        = KSPDestroy_LCD;
221:   ksp->ops->view           = KSPView_LCD;
222:   ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
223:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
224:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }