Actual source code: iguess.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/kspimpl.h>

  4: /* ---------------------------------------Method 1------------------------------------------------------------*/
  5: typedef struct {
  6:   PetscInt    method;   /* 1 or 2 */
  7:   PetscInt    curl;     /* Current number of basis vectors */
  8:   PetscInt    maxl;     /* Maximum number of basis vectors */
  9:   PetscInt    refcnt;
 10:   PetscBool   monitor;
 11:   Mat         mat;
 12:   KSP         ksp;
 13:   PetscScalar *alpha;   /* */
 14:   Vec         *xtilde;  /* Saved x vectors */
 15:   Vec         *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:   PetscMalloc1(maxl,&itg->alpha);
 31:   PetscLogObjectMemory((PetscObject)ksp,sizeof(KSPFischerGuess_Method1) + maxl*sizeof(PetscScalar));
 32:   KSPGetVecs(ksp,maxl,&itg->xtilde,0,NULL);
 33:   PetscLogObjectParents(ksp,maxl,itg->xtilde);
 34:   KSPGetVecs(ksp,maxl,&itg->btilde,0,NULL);
 35:   PetscLogObjectParents(ksp,maxl,itg->btilde);
 36:   VecDuplicate(itg->xtilde[0],&itg->guess);
 37:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)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 ",(double)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:   PetscInt    curl;     /* Current number of basis vectors */
133:   PetscInt    maxl;     /* Maximum number of basis vectors */
134:   PetscInt    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:   PetscMalloc1(maxl,&itg->alpha);
154:   PetscLogObjectMemory((PetscObject)ksp,sizeof(KSPFischerGuess_Method2) + maxl*sizeof(PetscScalar));
155:   KSPGetVecs(ksp,maxl,&itg->xtilde,0,NULL);
156:   PetscLogObjectParents(ksp,maxl,itg->xtilde);
157:   VecDuplicate(itg->xtilde[0],&itg->Ax);
158:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)itg->Ax);
159:   VecDuplicate(itg->xtilde[0],&itg->guess);
160:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)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 ",(double)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: {

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(PetscObjectComm((PetscObject)ksp),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;

294:   KSPGetOperators(ksp,&(*itg)->mat,NULL);
295:   return(0);
296: }

300: PetscErrorCode  KSPFischerGuessSetFromOptions(KSPFischerGuess ITG)
301: {

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

311: PetscErrorCode  KSPFischerGuessDestroy(KSPFischerGuess *ITG)
312: {

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

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

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

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

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

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

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

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