Actual source code: fischer.c

      1: #include <petsc/private/kspimpl.h>

      3: typedef struct {
      4:   PetscInt    method;   /* 1 or 2 */
      5:   PetscInt    curl;     /* Current number of basis vectors */
      6:   PetscInt    maxl;     /* Maximum number of basis vectors */
      7:   PetscBool   monitor;
      8:   PetscScalar *alpha;   /* */
      9:   Vec         *xtilde;  /* Saved x vectors */
     10:   Vec         *btilde;  /* Saved b vectors, method 1 */
     11:   Vec         Ax;       /* method 2 */
     12:   Vec         guess;
     13: } KSPGuessFischer;

     15: static PetscErrorCode KSPGuessReset_Fischer(KSPGuess guess)
     16: {
     17:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
     18:   PetscLayout     Alay = NULL,vlay = NULL;
     19:   PetscBool       cong;
     20:   PetscErrorCode  ierr;

     23:   itg->curl = 0;
     24:   /* destroy vectors if the size of the linear system has changed */
     25:   if (guess->A) {
     26:     MatGetLayouts(guess->A,&Alay,NULL);
     27:   }
     28:   if (itg->xtilde) {
     29:     VecGetLayout(itg->xtilde[0],&vlay);
     30:   }
     31:   cong = PETSC_FALSE;
     32:   if (vlay && Alay) {
     33:     PetscLayoutCompare(Alay,vlay,&cong);
     34:   }
     35:   if (!cong) {
     36:     VecDestroyVecs(itg->maxl,&itg->btilde);
     37:     VecDestroyVecs(itg->maxl,&itg->xtilde);
     38:     VecDestroy(&itg->guess);
     39:     VecDestroy(&itg->Ax);
     40:   }
     41:   return(0);
     42: }

     44: static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
     45: {
     46:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
     47:   PetscErrorCode  ierr;

     50:   if (!itg->alpha) {
     51:     PetscMalloc1(itg->maxl,&itg->alpha);
     52:     PetscLogObjectMemory((PetscObject)guess,itg->maxl*sizeof(PetscScalar));
     53:   }
     54:   if (!itg->xtilde) {
     55:     KSPCreateVecs(guess->ksp,itg->maxl,&itg->xtilde,0,NULL);
     56:     PetscLogObjectParents(guess,itg->maxl,itg->xtilde);
     57:   }
     58:   if (!itg->btilde && itg->method == 1) {
     59:     KSPCreateVecs(guess->ksp,itg->maxl,&itg->btilde,0,NULL);
     60:     PetscLogObjectParents(guess,itg->maxl,itg->btilde);
     61:   }
     62:   if (!itg->Ax && itg->method == 2) {
     63:     VecDuplicate(itg->xtilde[0],&itg->Ax);
     64:     PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->Ax);
     65:   }
     66:   if (!itg->guess) {
     67:     VecDuplicate(itg->xtilde[0],&itg->guess);
     68:     PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->guess);
     69:   }
     70:   return(0);
     71: }

     73: static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
     74: {
     75:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
     76:   PetscErrorCode  ierr;

     79:   PetscFree(itg->alpha);
     80:   VecDestroyVecs(itg->maxl,&itg->btilde);
     81:   VecDestroyVecs(itg->maxl,&itg->xtilde);
     82:   VecDestroy(&itg->guess);
     83:   VecDestroy(&itg->Ax);
     84:   PetscFree(itg);
     85:   PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",NULL);
     86:   return(0);
     87: }

     89: /* Note: do not change the b right hand side as is done in the publication */
     90: static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess,Vec b,Vec x)
     91: {
     92:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
     93:   PetscErrorCode  ierr;
     94:   PetscInt        i;

     97:   VecSet(x,0.0);
     98:   VecMDot(b,itg->curl,itg->btilde,itg->alpha);
     99:   if (itg->monitor) {
    100:     PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas = ");
    101:     for (i=0; i<itg->curl; i++) {
    102:       PetscPrintf(((PetscObject)guess)->comm,"%g ",(double)PetscAbsScalar(itg->alpha[i]));
    103:     }
    104:     PetscPrintf(((PetscObject)guess)->comm,"\n");
    105:   }
    106:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
    107:   VecCopy(x,itg->guess);
    108:   return(0);
    109: }

    111: static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
    112: {
    113:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
    114:   PetscReal       norm;
    115:   PetscErrorCode  ierr;
    116:   int             curl = itg->curl,i;

    119:   if (curl == itg->maxl) {
    120:     KSP_MatMult(guess->ksp,guess->A,x,itg->btilde[0]);
    121:     /* VecCopy(b,itg->btilde[0]); */
    122:     VecNormalize(itg->btilde[0],&norm);
    123:     VecCopy(x,itg->xtilde[0]);
    124:     VecScale(itg->xtilde[0],1.0/norm);
    125:     itg->curl = 1;
    126:   } else {
    127:     if (!curl) {
    128:       VecCopy(x,itg->xtilde[curl]);
    129:     } else {
    130:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
    131:     }
    132:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->btilde[curl]);
    133:     VecMDot(itg->btilde[curl],curl,itg->btilde,itg->alpha);
    134:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
    135:     VecMAXPY(itg->btilde[curl],curl,itg->alpha,itg->btilde);
    136:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
    137:     VecNormalize(itg->btilde[curl],&norm);
    138:     if (norm) {
    139:       VecScale(itg->xtilde[curl],1.0/norm);
    140:       itg->curl++;
    141:     } else {
    142:       PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
    143:     }
    144:   }
    145:   return(0);
    146: }

    148: /*
    149:   Given a basis generated already this computes a new guess x from the new right hand side b
    150:   Figures out the components of b in each btilde direction and adds them to x
    151:   Note: do not change the b right hand side as is done in the publication
    152: */
    153: static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
    154: {
    155:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
    156:   PetscErrorCode  ierr;
    157:   PetscInt        i;

    160:   VecSet(x,0.0);
    161:   VecMDot(b,itg->curl,itg->xtilde,itg->alpha);
    162:   if (itg->monitor) {
    163:     PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas = ");
    164:     for (i=0; i<itg->curl; i++) {
    165:       PetscPrintf(((PetscObject)guess)->comm,"%g ",(double)PetscAbsScalar(itg->alpha[i]));
    166:     }
    167:     PetscPrintf(((PetscObject)guess)->comm,"\n");
    168:   }
    169:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
    170:   VecCopy(x,itg->guess);
    171:   return(0);
    172: }

    174: static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
    175: {
    176:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
    177:   PetscScalar     norm;
    178:   PetscErrorCode  ierr;
    179:   int             curl = itg->curl,i;

    182:   if (curl == itg->maxl) {
    183:     KSP_MatMult(guess->ksp,guess->A,x,itg->Ax); /* norm = sqrt(x'Ax) */
    184:     VecDot(x,itg->Ax,&norm);
    185:     VecCopy(x,itg->xtilde[0]);
    186:     VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));
    187:     itg->curl = 1;
    188:   } else {
    189:     if (!curl) {
    190:       VecCopy(x,itg->xtilde[curl]);
    191:     } else {
    192:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
    193:     }
    194:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax);
    195:     VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);
    196:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
    197:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);

    199:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
    200:     VecDot(itg->xtilde[curl],itg->Ax,&norm);
    201:     if (PetscAbsScalar(norm) != 0.0) {
    202:       VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));
    203:       itg->curl++;
    204:     } else {
    205:       PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
    206:     }
    207:   }
    208:   return(0);
    209: }

    211: static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
    212: {
    213:   KSPGuessFischer *ITG = (KSPGuessFischer *)guess->data;
    214:   PetscInt        nmax = 2, model[2];
    215:   PetscBool       flg;
    216:   PetscErrorCode  ierr;

    219:   model[0] = ITG->method;
    220:   model[1] = ITG->maxl;
    221:   PetscOptionsBegin(PetscObjectComm((PetscObject)guess),((PetscObject)guess)->prefix,"Fischer guess options","KSPGuess");
    222:   PetscOptionsIntArray("-ksp_guess_fischer_model","Model type and dimension of basis","KSPGuessFischerSetModel",model,&nmax,&flg);
    223:   if (flg) {
    224:     KSPGuessFischerSetModel(guess,model[0],model[1]);
    225:   }
    226:   PetscOptionsBool("-ksp_guess_fischer_monitor","Monitor the guess",NULL,ITG->monitor,&ITG->monitor,NULL);
    227:   PetscOptionsEnd();
    228:   return(0);
    229: }

    231: static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess,PetscViewer viewer)
    232: {
    233:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
    234:   PetscBool       isascii;
    235:   PetscErrorCode  ierr;

    238:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
    239:   if (isascii) {
    240:     PetscViewerASCIIPrintf(viewer,"Model %D, size %D\n",itg->method,itg->maxl);
    241:   }
    242:   return(0);
    243: }

    245: /*@
    246:    KSPGuessFischerSetModel - Use the Paul Fischer algorithm

    248:    Logically Collective on guess

    250:    Input Parameters:
    251: +  guess - the initial guess context
    252: .  model - use model 1, model 2 or any other number to turn it off
    253: -  size  - size of subspace used to generate initial guess

    255:     Options Database:
    256: .   -ksp_guess_fischer_model <model,size> - uses the Fischer initial guess generator for repeated linear solves

    258:    Level: advanced

    260: .seealso: KSPGuess, KSPGuessCreate(), KSPSetUseFischerGuess(), KSPSetGuess(), KSPGetGuess(), KSP
    261: @*/
    262: PetscErrorCode  KSPGuessFischerSetModel(KSPGuess guess,PetscInt model,PetscInt size)
    263: {

    270:   PetscTryMethod(guess,"KSPGuessFischerSetModel_C",(KSPGuess,PetscInt,PetscInt),(guess,model,size));
    271:   return(0);
    272: }

    274: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess,PetscInt model,PetscInt size)
    275: {
    276:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;

    280:   if (model == 1) {
    281:     guess->ops->update    = KSPGuessUpdate_Fischer_1;
    282:     guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
    283:   } else if (model == 2) {
    284:     guess->ops->update    = KSPGuessUpdate_Fischer_2;
    285:     guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
    286:   } else {
    287:     guess->ops->update    = NULL;
    288:     guess->ops->formguess = NULL;
    289:     itg->method           = 0;
    290:     return(0);
    291:   }
    292:   if (size != itg->maxl) {
    293:     PetscFree(itg->alpha);
    294:     VecDestroyVecs(itg->maxl,&itg->btilde);
    295:     VecDestroyVecs(itg->maxl,&itg->xtilde);
    296:     VecDestroy(&itg->guess);
    297:     VecDestroy(&itg->Ax);
    298:   }
    299:   itg->method = model;
    300:   itg->maxl   = size;
    301:   return(0);
    302: }

    304: /*
    305:     KSPGUESSFISCHER - Implements Paul Fischer's initial guess algorithm Method 1 and 2 for situations where
    306:     a linear system is solved repeatedly

    308:   References:
    309: .   1. -   https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf

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

    318:     These are not intended to be used directly, they are called by KSP automatically with the command line options -ksp_guess_type fischer -ksp_guess_fischer_model <int,int> or programmatically as
    319: .vb
    320:     KSPGetGuess(ksp,&guess);
    321:     KSPGuessSetType(guess,KSPGUESSFISCHER);
    322:     KSPGuessFischerSetModel(guess,model,basis);

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

    326:     Developer note: the option -ksp_fischer_guess <int,int> is still available for backward compatibility

    328:     Level: intermediate

    330: @*/
    331: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
    332: {
    333:   KSPGuessFischer *fischer;
    334:   PetscErrorCode  ierr;

    337:   PetscNewLog(guess,&fischer);
    338:   fischer->method = 1;  /* defaults to method 1 */
    339:   fischer->maxl   = 10;
    340:   guess->data     = fischer;

    342:   guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
    343:   guess->ops->destroy        = KSPGuessDestroy_Fischer;
    344:   guess->ops->setup          = KSPGuessSetUp_Fischer;
    345:   guess->ops->view           = KSPGuessView_Fischer;
    346:   guess->ops->reset          = KSPGuessReset_Fischer;
    347:   guess->ops->update         = KSPGuessUpdate_Fischer_1;
    348:   guess->ops->formguess      = KSPGuessFormGuess_Fischer_1;

    350:   PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",KSPGuessFischerSetModel_Fischer);
    351:   return(0);
    352: }