Actual source code: iguess.c

petsc-3.3-p7 2013-05-11
  2: #include <petsc-private/kspimpl.h>

  4: /* ---------------------------------------Method 1------------------------------------------------------------*/
  5: typedef struct {
  6:     PetscInt    method,    /* 1 or 2 */
  7:                 curl,     /* Current number of basis vectors */
  8:                 maxl,     /* Maximum number of basis vectors */
  9:                 refcnt;
 10:     PetscBool   monitor;
 11:     Mat         mat;
 12:     KSP         ksp;
 13:     PetscScalar *alpha;   /* */
 14:     Vec         *xtilde,  /* Saved x vectors */
 15:                 *btilde;  /* Saved b vectors */
 16:     Vec         guess;
 17: } KSPFischerGuess_Method1;


 22: PetscErrorCode  KSPFischerGuessCreate_Method1(KSP ksp,int  maxl,KSPFischerGuess_Method1 **ITG)
 23: {
 24:   KSPFischerGuess_Method1 *itg;
 25:   PetscErrorCode          ierr;

 29:   PetscMalloc(sizeof(KSPFischerGuess_Method1),&itg);
 30:   PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);
 31:   PetscLogObjectMemory(ksp,sizeof(KSPFischerGuess_Method1) + maxl*sizeof(PetscScalar));
 32:   KSPGetVecs(ksp,maxl,&itg->xtilde,0,PETSC_NULL);
 33:   PetscLogObjectParents(ksp,maxl,itg->xtilde);
 34:   KSPGetVecs(ksp,maxl,&itg->btilde,0,PETSC_NULL);
 35:   PetscLogObjectParents(ksp,maxl,itg->btilde);
 36:   VecDuplicate(itg->xtilde[0],&itg->guess);
 37:   PetscLogObjectParent(ksp,itg->guess);
 38:   *ITG = itg;
 39:   return(0);
 40: }


 45: PetscErrorCode  KSPFischerGuessDestroy_Method1(KSPFischerGuess_Method1 *itg)
 46: {

 50:   PetscFree(itg->alpha);
 51:   VecDestroyVecs(itg->maxl,&itg->btilde);
 52:   VecDestroyVecs(itg->maxl,&itg->xtilde);
 53:   VecDestroy(&itg->guess);
 54:   PetscFree(itg);
 55:   return(0);
 56: }


 59: /*
 60:         Given a basis generated already this computes a new guess x from the new right hand side b
 61:      Figures out the components of b in each btilde direction and adds them to x
 62: */
 65: PetscErrorCode  KSPFischerGuessFormGuess_Method1(KSPFischerGuess_Method1 *itg,Vec b,Vec x)
 66: {
 68:   PetscInt       i;

 73:   VecSet(x,0.0);
 74:   VecMDot(b,itg->curl,itg->btilde,itg->alpha);
 75:   if (itg->monitor) {
 76:     PetscPrintf(((PetscObject)itg->ksp)->comm,"KSPFischerGuess alphas = ");
 77:     for (i=0; i<itg->curl; i++ ){
 78:       PetscPrintf(((PetscObject)itg->ksp)->comm,"%G ",PetscAbsScalar(itg->alpha[i]));
 79:     }
 80:     PetscPrintf(((PetscObject)itg->ksp)->comm,"\n");
 81:   }
 82:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
 83:   VecCopy(x,itg->guess);
 84:   /* Note: do not change the b right hand side as is done in the publication */
 85:   return(0);
 86: }

 90: PetscErrorCode  KSPFischerGuessUpdate_Method1(KSPFischerGuess_Method1 *itg,Vec x)
 91: {
 92:   PetscReal      norm;
 94:   int            curl = itg->curl,i;

 99:   if (curl == itg->maxl) {
100:     KSP_MatMult(itg->ksp,itg->mat,x,itg->btilde[0]);
101:     VecNormalize(itg->btilde[0],&norm);
102:     VecCopy(x,itg->xtilde[0]);
103:     VecScale(itg->xtilde[0],1.0/norm);
104:     itg->curl = 1;
105:   } else {
106:     if (!curl) {
107:       VecCopy(x,itg->xtilde[curl]);
108:     } else {
109:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
110:     }

112:     KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->btilde[curl]);
113:     VecMDot(itg->btilde[curl],curl,itg->btilde,itg->alpha);
114:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
115:     VecMAXPY(itg->btilde[curl],curl,itg->alpha,itg->btilde);
116:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);

118:     VecNormalize(itg->btilde[curl],&norm);
119:     if (norm) {
120:       VecScale(itg->xtilde[curl],1.0/norm);
121:       itg->curl++;
122:     } else {
123:       PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous");
124:     }
125:   }
126:   return(0);
127: }

129: /* ---------------------------------------Method 2------------------------------------------------------------*/
130: typedef struct {
131:     PetscInt    method,    /* 1 or 2 */
132:                 curl,     /* Current number of basis vectors */
133:                 maxl,     /* Maximum number of basis vectors */
134:                 refcnt;
135:     PetscBool   monitor;
136:     Mat         mat;
137:     KSP         ksp;
138:     PetscScalar *alpha;   /* */
139:     Vec         *xtilde;  /* Saved x vectors */
140:   Vec         Ax,guess;
141: } KSPFischerGuess_Method2;

145: PetscErrorCode  KSPFischerGuessCreate_Method2(KSP ksp,int  maxl,KSPFischerGuess_Method2 **ITG)
146: {
147:   KSPFischerGuess_Method2 *itg;
148:   PetscErrorCode          ierr;

152:   PetscMalloc(sizeof(KSPFischerGuess_Method2),&itg);
153:   PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);
154:   PetscLogObjectMemory(ksp,sizeof(KSPFischerGuess_Method2) + maxl*sizeof(PetscScalar));
155:   KSPGetVecs(ksp,maxl,&itg->xtilde,0,PETSC_NULL);
156:   PetscLogObjectParents(ksp,maxl,itg->xtilde);
157:   VecDuplicate(itg->xtilde[0],&itg->Ax);
158:   PetscLogObjectParent(ksp,itg->Ax);
159:   VecDuplicate(itg->xtilde[0],&itg->guess);
160:   PetscLogObjectParent(ksp,itg->guess);
161:   *ITG = itg;
162:   return(0);
163: }


168: PetscErrorCode  KSPFischerGuessDestroy_Method2(KSPFischerGuess_Method2 *itg)
169: {

173:   PetscFree(itg->alpha);
174:   VecDestroyVecs(itg->maxl,&itg->xtilde);
175:   VecDestroy(&itg->Ax);
176:   VecDestroy(&itg->guess);
177:   PetscFree(itg);
178:   return(0);
179: }


182: /*
183:         Given a basis generated already this computes a new guess x from the new right hand side b
184:      Figures out the components of b in each btilde direction and adds them to x
185: */
188: PetscErrorCode  KSPFischerGuessFormGuess_Method2(KSPFischerGuess_Method2 *itg,Vec b,Vec x)
189: {
191:   PetscInt       i;

196:   VecSet(x,0.0);
197:   VecMDot(b,itg->curl,itg->xtilde,itg->alpha);
198:   if (itg->monitor) {
199:     PetscPrintf(((PetscObject)itg->ksp)->comm,"KSPFischerGuess alphas = ");
200:     for (i=0; i<itg->curl; i++ ){
201:       PetscPrintf(((PetscObject)itg->ksp)->comm,"%G ",PetscAbsScalar(itg->alpha[i]));
202:     }
203:     PetscPrintf(((PetscObject)itg->ksp)->comm,"\n");
204:   }
205:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
206:   VecCopy(x,itg->guess);
207:   /* Note: do not change the b right hand side as is done in the publication */
208:   return(0);
209: }

213: PetscErrorCode  KSPFischerGuessUpdate_Method2(KSPFischerGuess_Method2 *itg,Vec x)
214: {
215:   PetscScalar    norm;
217:   int            curl = itg->curl,i;

222:   if (curl == itg->maxl) {
223:     KSP_MatMult(itg->ksp,itg->mat,x,itg->Ax); /* norm = sqrt(x'Ax) */
224:     VecDot(x,itg->Ax,&norm);
225:     VecCopy(x,itg->xtilde[0]);
226:     VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));
227:     itg->curl = 1;
228:   } else {
229:     if (!curl) {
230:       VecCopy(x,itg->xtilde[curl]);
231:     } else {
232:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
233:     }
234:     KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);
235:     VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);
236:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
237:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);

239:     KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
240:     VecDot(itg->xtilde[curl],itg->Ax,&norm);
241:     if (PetscAbsScalar(norm) != 0.0) {
242:       VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));
243:       itg->curl++;
244:     } else {
245:       PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous");
246:     }
247:   }
248:   return(0);
249: }

251: /* ---------------------------------------------------------------------------------------------------------*/
254: /*@C
255:     KSPFischerGuessCreate - Implements Paul Fischer's initial guess algorithm Method 1 and 2 for situations where
256:     a linear system is solved repeatedly 

258:   References:
259:       http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf

261:    Notes: the algorithm is different from the paper because we do not CHANGE the right hand side of the new 
262:     problem and solve the problem with an initial guess of zero, rather we solve the original new problem
263:     with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
264:     the original RHS.) But we use the xtilde = x - xguess as the new direction so that it is not
265:     mostly orthogonal to the previous solutions.

267:     These are not intended to be used directly, they are called by KSP automatically when the 
268:     KSP option KSPSetFischerGuess(KSP,PetscInt,PetscInt) or -ksp_guess_fischer <int,int>

270:     Method 2 is only for positive definite matrices, since it uses the A norm.

272:     This is not currently programmed as a PETSc class because there are only two methods; if more methods
273:     are introduced it should be changed. For example the Knoll guess should be included

275:     Level: advanced

277: @*/
278: PetscErrorCode  KSPFischerGuessCreate(KSP ksp,PetscInt method,PetscInt maxl,KSPFischerGuess *itg)
279: {
280:   PetscErrorCode  ierr;

283:   if (method == 1) {
284:     KSPFischerGuessCreate_Method1(ksp,maxl,(KSPFischerGuess_Method1 **)itg);
285:   } else if (method == 2) {
286:     KSPFischerGuessCreate_Method2(ksp,maxl,(KSPFischerGuess_Method2 **)itg);
287:   } else SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
288:   (*itg)->method  = method;
289:   (*itg)->curl    = 0;
290:   (*itg)->maxl    = maxl;
291:   (*itg)->ksp     = ksp;
292:   (*itg)->refcnt  = 1;
293:   KSPGetOperators(ksp,&(*itg)->mat,PETSC_NULL,PETSC_NULL);
294:   return(0);
295: }

299: PetscErrorCode  KSPFischerGuessSetFromOptions(KSPFischerGuess ITG)
300: {
301:   PetscErrorCode  ierr;

304:   PetscOptionsGetBool(((PetscObject)ITG->ksp)->prefix,"-ksp_fischer_guess_monitor",&ITG->monitor,PETSC_NULL);
305:   return(0);
306: }

310: PetscErrorCode  KSPFischerGuessDestroy(KSPFischerGuess *ITG)
311: {
312:   PetscErrorCode  ierr;

315:   if (!*ITG) return(0);
316:   if (--(*ITG)->refcnt) {*ITG = 0; return(0);}

318:   if ((*ITG)->method == 1) {
319:     KSPFischerGuessDestroy_Method1((KSPFischerGuess_Method1 *)*ITG);
320:   } else if ((*ITG)->method == 2) {
321:     KSPFischerGuessDestroy_Method2((KSPFischerGuess_Method2 *)*ITG);
322:   } else SETERRQ(((PetscObject)(*ITG)->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
323:   *ITG = PETSC_NULL;
324:   return(0);
325: }

329: PetscErrorCode  KSPFischerGuessUpdate(KSPFischerGuess itg,Vec x)
330: {
331:   PetscErrorCode  ierr;

334:   if (itg->method == 1) {
335:     KSPFischerGuessUpdate_Method1((KSPFischerGuess_Method1 *)itg,x);
336:   } else if (itg->method == 2) {
337:     KSPFischerGuessUpdate_Method2((KSPFischerGuess_Method2 *)itg,x);
338:   } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
339:   return(0);
340: }

344: PetscErrorCode  KSPFischerGuessFormGuess(KSPFischerGuess itg,Vec b,Vec x)
345: {
346:   PetscErrorCode  ierr;

349:   if (itg->method == 1) {
350:     KSPFischerGuessFormGuess_Method1((KSPFischerGuess_Method1 *)itg,b,x);
351:   } else if (itg->method == 2) {
352:     KSPFischerGuessFormGuess_Method2((KSPFischerGuess_Method2 *)itg,b,x);
353:   } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
354:   return(0);
355: }

359: /*
360:     KSPFischerGuessReset - This is called whenever KSPSetOperators() is called to tell the
361:       initial guess object that the matrix is changed and so the initial guess object
362:       must restart from scratch building the subspace where the guess is computed from.
363: */
364: PetscErrorCode  KSPFischerGuessReset(KSPFischerGuess itg)
365: {

369:   itg->curl = 0;
370:   KSPGetOperators(itg->ksp,&itg->mat,PETSC_NULL,PETSC_NULL);
371:   return(0);
372: }