Actual source code: lcd.c

petsc-3.4.4 2014-03-13
  2: #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>

  6: PetscErrorCode KSPSetUp_LCD(KSP ksp)
  7: {
  8:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;
 10:   PetscInt       restart = lcd->restart;

 13:   /* get work vectors needed by LCD */
 14:   KSPSetWorkVecs(ksp,2);

 16:   VecDuplicateVecs(ksp->work[0],restart+1,&lcd->P);
 17:   VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q);
 18:   PetscLogObjectMemory(ksp,2*(restart+2)*sizeof(Vec));
 19:   return(0);
 20: }

 22: /*     KSPSolve_LCD - This routine actually applies the left conjugate
 23:     direction method

 25:    Input Parameter:
 26: .     ksp - the Krylov space object that was set to use LCD, by, for
 27:             example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPLCD);

 29:    Output Parameter:
 30: .     its - number of iterations used

 32: */
 35: PetscErrorCode  KSPSolve_LCD(KSP ksp)
 36: {
 38:   PetscInt       it,j,max_k;
 39:   PetscScalar    alfa, beta, num, den, mone;
 40:   PetscReal      rnorm;
 41:   Vec            X,B,R,Z;
 42:   KSP_LCD        *lcd;
 43:   Mat            Amat,Pmat;
 44:   MatStructure   pflag;
 45:   PetscBool      diagonalscale;

 48:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
 49:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

 51:   lcd   = (KSP_LCD*)ksp->data;
 52:   X     = ksp->vec_sol;
 53:   B     = ksp->vec_rhs;
 54:   R     = ksp->work[0];
 55:   Z     = ksp->work[1];
 56:   max_k = lcd->restart;
 57:   mone  = -1;

 59:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);

 61:   ksp->its = 0;
 62:   if (!ksp->guess_zero) {
 63:     KSP_MatMult(ksp,Amat,X,Z);             /*   z <- b - Ax       */
 64:     VecAYPX(Z,mone,B);
 65:   } else {
 66:     VecCopy(B,Z);                         /*     z <- b (x is 0) */
 67:   }

 69:   KSP_PCApply(ksp,Z,R);                   /*     r <- M^-1z         */
 70:   VecNorm(R,NORM_2,&rnorm);
 71:   KSPLogResidualHistory(ksp,rnorm);
 72:   KSPMonitor(ksp,0,rnorm);
 73:   ksp->rnorm = rnorm;

 75:   /* test for convergence */
 76:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
 77:   if (ksp->reason) return(0);

 79:   it = 0;
 80:   VecCopy(R,lcd->P[0]);

 82:   while (!ksp->reason && ksp->its < ksp->max_it) {
 83:     it   = 0;
 84:     KSP_MatMult(ksp,Amat,lcd->P[it],Z);
 85:     KSP_PCApply(ksp,Z,lcd->Q[it]);

 87:     while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
 88:       ksp->its++;
 89:       VecDot(lcd->P[it],R,&num);
 90:       VecDot(lcd->P[it],lcd->Q[it], &den);
 91:       alfa = num/den;
 92:       VecAXPY(X,alfa,lcd->P[it]);
 93:       VecAXPY(R,-alfa,lcd->Q[it]);
 94:       VecNorm(R,NORM_2,&rnorm);

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

101:       if (ksp->reason) break;

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

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

125: */
128: PetscErrorCode KSPReset_LCD(KSP ksp)
129: {
130:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;

134:   if (lcd->P) { VecDestroyVecs(lcd->restart+1,&lcd->P);}
135:   if (lcd->Q) { VecDestroyVecs(lcd->restart+1,&lcd->Q);}
136:   return(0);
137: }


142: PetscErrorCode KSPDestroy_LCD(KSP ksp)
143: {

147:   KSPReset_LCD(ksp);
148:   PetscFree(ksp->data);
149:   return(0);
150: }

152: /*
153:      KSPView_LCD - Prints information about the current Krylov method being used

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

159: */
162: PetscErrorCode KSPView_LCD(KSP ksp,PetscViewer viewer)
163: {

165:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;
167:   PetscBool      iascii;

170:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
171:   if (iascii) {
172:     PetscViewerASCIIPrintf(viewer,"  LCD: restart=%d\n",lcd->restart);
173:     PetscViewerASCIIPrintf(viewer,"  LCD: happy breakdown tolerance %g\n",lcd->haptol);
174:   }
175:   return(0);
176: }

178: /*
179:     KSPSetFromOptions_LCD - Checks the options database for options related to the
180:                             LCD method.
181: */
184: PetscErrorCode KSPSetFromOptions_LCD(KSP ksp)
185: {
187:   PetscBool      flg;
188:   KSP_LCD        *lcd = (KSP_LCD*)ksp->data;

191:   PetscOptionsHead("KSP LCD options");
192:   PetscOptionsInt("-ksp_lcd_restart","Number of vectors conjugate","KSPLCDSetRestart",lcd->restart,&lcd->restart,&flg);
193:   if (flg && lcd->restart < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
194:   PetscOptionsReal("-ksp_lcd_haptol","Tolerance for exact convergence (happy ending)","KSPLCDSetHapTol",lcd->haptol,&lcd->haptol,&flg);
195:   if (flg && lcd->haptol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
196:   return(0);
197: }

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

202:    Options Database Keys:
203: +   -ksp_lcd_restart - number of vectors conjudate
204: -   -ksp_lcd_haptol - tolerance for exact convergence (happing ending)

206:    Level: beginner

208:     Notes: Support only for left preconditioning

210:     References:
211:    - J.Y. Yuan, G.H.Golub, R.J. Plemmons, and W.A.G. Cecilio. Semiconjugate
212:      direction methods for real positive definite system. BIT Numerical
213:      Mathematics, 44(1):189-207,2004.
214:    - Y. Dai and J.Y. Yuan. Study on semi-conjugate direction methods for
215:      non-symmetric systems. International Journal for Numerical Methods in
216:      Engineering, 60:1383-1399,2004.
217:    - L. Catabriga, A.L.G.A. Coutinho, and L.P.Franca. Evaluating the LCD
218:      algorithm for solving linear systems of equations arising from implicit
219:      SUPG formulation of compressible flows. International Journal for
220:      Numerical Methods in Engineering, 60:1513-1534,2004
221:    - L. Catabriga, A. M. P. Valli, B. Z. Melotti, L. M. Pessoa,
222:      A. L. G. A. Coutinho, Performance of LCD iterative method in the finite
223:      element and finite difference solution of convection-diffusion
224:      equations,  Communications in Numerical Methods in Engineering, (Early
225:      View).

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


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

233: M*/

237: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
238: {
240:   KSP_LCD        *lcd;

243:   PetscNewLog(ksp,KSP_LCD,&lcd);
244:   ksp->data    = (void*)lcd;
245:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
246:   lcd->restart = 30;
247:   lcd->haptol  = 1.0e-30;

249:   /*
250:        Sets the functions that are associated with this data structure
251:        (in C++ this is the same as defining virtual functions)
252:   */
253:   ksp->ops->setup          = KSPSetUp_LCD;
254:   ksp->ops->solve          = KSPSolve_LCD;
255:   ksp->ops->reset          = KSPReset_LCD;
256:   ksp->ops->destroy        = KSPDestroy_LCD;
257:   ksp->ops->view           = KSPView_LCD;
258:   ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
259:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
260:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
261:   return(0);
262: }