Actual source code: iguess.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/ksp/kspimpl.h
      4: /*
      5:   This code inplements Paul Fischer's initial guess code for situations where
      6:   a linear system is solved repeatedly
      7:  */

      9: typedef struct {
     10:     int         curl,     /* Current number of basis vectors */
     11:                 maxl;     /* Maximum number of basis vectors */
     12:     PetscScalar *alpha;   /* */
     13:     Vec         *xtilde,  /* Saved x vectors */
     14:                 *btilde;  /* Saved b vectors */
     15: } KSPIGUESS;

     19: PetscErrorCode PETSCKSP_DLLEXPORT KSPGuessCreate(KSP ksp,int  maxl,void **ITG)
     20: {
     21:   KSPIGUESS *itg;

     24:   *ITG = 0;
     27:   PetscMalloc(sizeof(KSPIGUESS),&itg);
     28:   itg->curl = 0;
     29:   itg->maxl = maxl;
     30:   PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);
     31:   PetscLogObjectMemory(ksp,sizeof(KSPIGUESS) + maxl*sizeof(PetscScalar));
     32:   KSPGetVecs(ksp,maxl,&itg->xtilde);
     33:   PetscLogObjectParents(ksp,maxl,itg->xtilde);
     34:   KSPGetVecs(ksp,maxl,&itg->btilde);
     35:   PetscLogObjectParents(ksp,maxl,itg->btilde);
     36:   *ITG = (void*)itg;
     37:   return(0);
     38: }

     42: PetscErrorCode PETSCKSP_DLLEXPORT KSPGuessDestroy(KSP ksp,KSPIGUESS *itg)
     43: {

     48:   PetscFree(itg->alpha);
     49:   VecDestroyVecs(itg->btilde,itg->maxl);
     50:   VecDestroyVecs(itg->xtilde,itg->maxl);
     51:   PetscFree(itg);
     52:   return(0);
     53: }

     57: PetscErrorCode PETSCKSP_DLLEXPORT KSPGuessFormB(KSP ksp,KSPIGUESS *itg,Vec b)
     58: {
     60:   int         i;
     61:   PetscScalar tmp;

     67:   for (i=1; i<=itg->curl; i++) {
     68:     VecDot(itg->btilde[i-1],b,&(itg->alpha[i-1]));
     69:     tmp = -itg->alpha[i-1];
     70:     VecAXPY(b,tmp,itg->btilde[i-1]);
     71:   }
     72:   return(0);
     73: }

     77: PetscErrorCode PETSCKSP_DLLEXPORT KSPGuessFormX(KSP ksp,KSPIGUESS *itg,Vec x)
     78: {
     80:   int i;

     86:   VecCopy(x,itg->xtilde[itg->curl]);
     87:   for (i=1; i<=itg->curl; i++) {
     88:     VecAXPY(x,itg->alpha[i-1],itg->xtilde[i-1]);
     89:   }
     90:   return(0);
     91: }

     95: PetscErrorCode PETSCKSP_DLLEXPORT KSPGuessUpdate(KSP ksp,Vec x,KSPIGUESS *itg)
     96: {
     97:   PetscReal    normax,norm;
     98:   PetscScalar  tmp;
     99:   MatStructure pflag;
    101:   int          curl = itg->curl,i;
    102:   Mat          Amat,Pmat;

    108:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
    109:   if (curl == itg->maxl) {
    110:     KSP_MatMult(ksp,Amat,x,itg->btilde[0]);
    111:     VecNorm(itg->btilde[0],NORM_2,&normax);
    112:     tmp = 1.0/normax; VecScale(itg->btilde[0],tmp);
    113:     /* VCOPY(ksp->vc,x,itg->xtilde[0]); */
    114:     VecScale(itg->xtilde[0],tmp);
    115:   } else {
    116:     KSP_MatMult(ksp,Amat,itg->xtilde[curl],itg->btilde[curl]);
    117:     for (i=1; i<=curl; i++) {
    118:       VecDot(itg->btilde[curl],itg->btilde[i-1],itg->alpha+i-1);
    119:     }
    120:     for (i=1; i<=curl; i++) {
    121:       tmp  = -itg->alpha[i-1];
    122:       VecAXPY(itg->btilde[curl],tmp,itg->btilde[i-1]);
    123:       VecAXPY(itg->xtilde[curl],itg->alpha[i-1],itg->xtilde[i-1]);
    124:     }
    125:     VecNorm(itg->btilde[curl],NORM_2,&norm);
    126:     tmp = 1.0/norm; VecScale(itg->btilde[curl],tmp);
    127:     VecNorm(itg->xtilde[curl],NORM_2,&norm);
    128:     VecScale(itg->xtilde[curl],tmp);
    129:     itg->curl++;
    130:   }
    131:   return(0);
    132: }