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:   PetscLogObjectMemory((PetscObject)ksp,2*(restart+2)*sizeof(Vec));
 15:   return 0;
 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:   PCGetDiagonalScale(ksp->pc,&diagonalscale);

 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:   PCGetOperators(ksp->pc,&Amat,&Pmat);

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

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

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

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

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

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

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

 98:       if (ksp->reason) break;

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

104:       for (j = 0; j <= it; j++) {
105:         VecDot(lcd->P[j],lcd->Q[it+1],&num);
106:         KSPCheckDot(ksp,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: */
124: PetscErrorCode KSPReset_LCD(KSP ksp)
125: {
126:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;

128:   if (lcd->P) VecDestroyVecs(lcd->restart+1,&lcd->P);
129:   if (lcd->Q) VecDestroyVecs(lcd->restart+1,&lcd->Q);
130:   return 0;
131: }

133: PetscErrorCode KSPDestroy_LCD(KSP ksp)
134: {
135:   KSPReset_LCD(ksp);
136:   PetscFree(ksp->data);
137:   return 0;
138: }

140: /*
141:      KSPView_LCD - Prints information about the current Krylov method being used

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

147: */
148: PetscErrorCode KSPView_LCD(KSP ksp,PetscViewer viewer)
149: {

151:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;
152:   PetscBool      iascii;

154:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
155:   if (iascii) {
156:     PetscViewerASCIIPrintf(viewer,"  restart=%d\n",lcd->restart);
157:     PetscViewerASCIIPrintf(viewer,"  happy breakdown tolerance %g\n",lcd->haptol);
158:   }
159:   return 0;
160: }

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

171:   PetscOptionsHead(PetscOptionsObject,"KSP LCD options");
172:   PetscOptionsInt("-ksp_lcd_restart","Number of vectors conjugate","KSPLCDSetRestart",lcd->restart,&lcd->restart,&flg);
174:   PetscOptionsReal("-ksp_lcd_haptol","Tolerance for exact convergence (happy ending)","KSPLCDSetHapTol",lcd->haptol,&lcd->haptol,&flg);
176:   return 0;
177: }

179: /*MC
180:      KSPLCD -  Implements the LCD (left conjugate direction) method in PETSc.

182:    Options Database Keys:
183: +   -ksp_lcd_restart - number of vectors conjudate
184: -   -ksp_lcd_haptol - tolerance for exact convergence (happing ending)

186:    Level: beginner

188:     Notes:
189:     Support only for left preconditioning

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

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

210: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
211:            KSPCGSetType(), KSPLCDSetRestart(), KSPLCDSetHapTol()

213: M*/

215: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
216: {
217:   KSP_LCD        *lcd;

219:   PetscNewLog(ksp,&lcd);
220:   ksp->data    = (void*)lcd;
221:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
222:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
223:   lcd->restart = 30;
224:   lcd->haptol  = 1.0e-30;

226:   /*
227:        Sets the functions that are associated with this data structure
228:        (in C++ this is the same as defining virtual functions)
229:   */
230:   ksp->ops->setup          = KSPSetUp_LCD;
231:   ksp->ops->solve          = KSPSolve_LCD;
232:   ksp->ops->reset          = KSPReset_LCD;
233:   ksp->ops->destroy        = KSPDestroy_LCD;
234:   ksp->ops->view           = KSPView_LCD;
235:   ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
236:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
237:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
238:   return 0;
239: }