Actual source code: fgmres.c

      1: #define PETSCKSP_DLL

      3: /*
      4:     This file implements FGMRES (a Generalized Minimal Residual) method.
      5:     Reference:  Saad, 1993.

      7:     Preconditioning:  It the preconditioner is constant then this fgmres
      8:     code is equivalent to RIGHT-PRECONDITIONED GMRES.

     10:     Restarts:  Restarts are basically solves with x0 not equal to zero.
     11:
     12:        Contributed by Allison Baker

     14: */

     16:  #include src/ksp/ksp/impls/gmres/fgmres/fgmresp.h
     17: #define FGMRES_DELTA_DIRECTIONS 10
     18: #define FGMRES_DEFAULT_MAXK     30
     19: static PetscErrorCode FGMRESGetNewVectors(KSP,PetscInt);
     20: static PetscErrorCode FGMRESUpdateHessenberg(KSP,PetscInt,PetscTruth,PetscReal *);
     21: static PetscErrorCode BuildFgmresSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);

     23: EXTERN PetscErrorCode KSPView_GMRES(KSP,PetscViewer);
     24: /*

     26:     KSPSetUp_FGMRES - Sets up the workspace needed by fgmres.

     28:     This is called once, usually automatically by KSPSolveQ() or KSPSetUp(),
     29:     but can be called directly by KSPSetUp().

     31: */
     34: PetscErrorCode    KSPSetUp_FGMRES(KSP ksp)
     35: {
     36:   PetscInt       size,hh,hes,rs,cc;
     38:   PetscInt       max_k,k;
     39:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)ksp->data;

     42:   if (ksp->pc_side == PC_SYMMETRIC) {
     43:     SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPFGMRES");
     44:   } else if (ksp->pc_side == PC_LEFT) {
     45:     SETERRQ(PETSC_ERR_SUP,"no left preconditioning for KSPFGMRES");
     46:   }
     47:   max_k         = fgmres->max_k;
     48:   hh            = (max_k + 2) * (max_k + 1);
     49:   hes           = (max_k + 1) * (max_k + 1);
     50:   rs            = (max_k + 2);
     51:   cc            = (max_k + 1);  /* SS and CC are the same size */
     52:   size          = (hh + hes + rs + 2*cc) * sizeof(PetscScalar);

     54:   /* Allocate space and set pointers to beginning */
     55:   PetscMalloc(size,&fgmres->hh_origin);
     56:   PetscMemzero(fgmres->hh_origin,size);
     57:   PetscLogObjectMemory(ksp,size); /* HH - modified (by plane rotations) hessenburg */
     58:   fgmres->hes_origin = fgmres->hh_origin + hh;     /* HES - unmodified hessenburg */
     59:   fgmres->rs_origin  = fgmres->hes_origin + hes;   /* RS - the right-hand-side of the
     60:                                                       Hessenberg system */
     61:   fgmres->cc_origin  = fgmres->rs_origin + rs;     /* CC - cosines for rotations */
     62:   fgmres->ss_origin  = fgmres->cc_origin + cc;     /* SS - sines for rotations */

     64:   if (ksp->calc_sings) {
     65:     /* Allocate workspace to hold Hessenberg matrix needed by Eispack */
     66:     size = (max_k + 3)*(max_k + 9)*sizeof(PetscScalar);
     67:     PetscMalloc(size,&fgmres->Rsvd);
     68:     PetscMalloc(5*(max_k+2)*sizeof(PetscReal),&fgmres->Dsvd);
     69:     PetscLogObjectMemory(ksp,size+5*(max_k+2)*sizeof(PetscReal));
     70:   }

     72:   /* Allocate array to hold pointers to user vectors.  Note that we need
     73:    4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
     74:   PetscMalloc((VEC_OFFSET+2+max_k)*sizeof(void*),&fgmres->vecs);
     75:   fgmres->vecs_allocated = VEC_OFFSET + 2 + max_k;
     76:   PetscMalloc((VEC_OFFSET+2+max_k)*sizeof(void*),&fgmres->user_work);
     77:   PetscMalloc((VEC_OFFSET+2+max_k)*sizeof(PetscInt),&fgmres->mwork_alloc);
     78:   PetscLogObjectMemory(ksp,(VEC_OFFSET+2+max_k)*(2*sizeof(void*)+sizeof(PetscInt)));

     80:   /* New for FGMRES - Allocate array to hold pointers to preconditioned
     81:      vectors - same sizes as user vectors above */
     82:   PetscMalloc((VEC_OFFSET+2+max_k)*sizeof(void*),&fgmres->prevecs);
     83:   PetscMalloc((VEC_OFFSET+2+max_k)*sizeof(void*),&fgmres->prevecs_user_work);
     84:   PetscLogObjectMemory(ksp,(VEC_OFFSET+2+max_k)*(2*sizeof(void*)));

     87:   /* if q_preallocate = 0 then only allocate one "chunck" of space (for
     88:      5 vectors) - additional will then be allocated from FGMREScycle()
     89:      as needed.  Otherwise, allocate all of the space that could be needed */
     90:   if (fgmres->q_preallocate) {
     91:     fgmres->vv_allocated   = VEC_OFFSET + 2 + max_k;
     92:   } else {
     93:     fgmres->vv_allocated    = 5;
     94:   }

     96:   /* space for work vectors */
     97:   KSPGetVecs(ksp,fgmres->vv_allocated,&fgmres->user_work[0]);
     98:   PetscLogObjectParents(ksp,fgmres->vv_allocated,fgmres->user_work[0]);
     99:   for (k=0; k < fgmres->vv_allocated; k++) {
    100:     fgmres->vecs[k] = fgmres->user_work[0][k];
    101:   }

    103:   /* space for preconditioned vectors */
    104:   KSPGetVecs(ksp,fgmres->vv_allocated,&fgmres->prevecs_user_work[0]);
    105:   PetscLogObjectParents(ksp,fgmres->vv_allocated,fgmres->prevecs_user_work[0]);
    106:   for (k=0; k < fgmres->vv_allocated; k++) {
    107:     fgmres->prevecs[k] = fgmres->prevecs_user_work[0][k];
    108:   }

    110:   /* specify how many work vectors have been allocated in this
    111:      chunck" (the first one) */
    112:   fgmres->mwork_alloc[0] = fgmres->vv_allocated;
    113:   fgmres->nwork_alloc    = 1;

    115:   return(0);
    116: }

    118: /*
    119:     FGMRESResidual - This routine computes the initial residual (NOT PRECONDITIONED)
    120: */
    123: static PetscErrorCode FGMRESResidual(KSP ksp)
    124: {
    125:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)(ksp->data);
    126:   PetscScalar    mone = -1.0;
    127:   Mat            Amat,Pmat;
    128:   MatStructure   pflag;

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

    134:   /* put A*x into VEC_TEMP */
    135:   MatMult(Amat,ksp->vec_sol,VEC_TEMP);
    136:   /* now put residual (-A*x + f) into vec_vv(0) */
    137:   VecWAXPY(VEC_VV(0),mone,VEC_TEMP,ksp->vec_rhs);
    138:   return(0);
    139: }

    141: /*

    143:     FGMRESCycle - Run fgmres, possibly with restart.  Return residual
    144:                   history if requested.

    146:     input parameters:
    147: .         fgmres  - structure containing parameters and work areas

    149:     output parameters:
    150: .        itcount - number of iterations used.  If null, ignored.
    151: .        converged - 0 if not converged

    153:
    154:     Notes:
    155:     On entry, the value in vector VEC_VV(0) should be
    156:     the initial residual.

    159:  */
    162: PetscErrorCode FGMREScycle(PetscInt *itcount,KSP ksp)
    163: {

    165:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)(ksp->data);
    166:   PetscReal      res_norm;
    167:   PetscReal      hapbnd,tt;
    168:   PetscScalar    zero = 0.0;
    169:   PetscScalar    tmp;
    170:   PetscTruth     hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
    172:   PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
    173:   PetscInt       max_k = fgmres->max_k; /* max # of directions Krylov space */
    174:   Mat            Amat,Pmat;
    175:   MatStructure   pflag;

    179:   /* Number of pseudo iterations since last restart is the number
    180:      of prestart directions */
    181:   loc_it = 0;

    183:   /* note: (fgmres->it) is always set one less than (loc_it) It is used in
    184:      KSPBUILDSolution_FGMRES, where it is passed to BuildFGmresSoln.
    185:      Note that when BuildFGmresSoln is called from this function,
    186:      (loc_it -1) is passed, so the two are equivalent */
    187:   fgmres->it = (loc_it - 1);

    189:   /* initial residual is in VEC_VV(0)  - compute its norm*/
    190:   VecNorm(VEC_VV(0),NORM_2,&res_norm);

    192:   /* first entry in right-hand-side of hessenberg system is just
    193:      the initial residual norm */
    194:   *RS(0) = res_norm;

    196:   /* FYI: AMS calls are for memory snooper */
    197:   PetscObjectTakeAccess(ksp);
    198:   ksp->rnorm = res_norm;
    199:   PetscObjectGrantAccess(ksp);
    200:   KSPLogResidualHistory(ksp,res_norm);

    202:   /* check for the convergence - maybe the current guess is good enough */
    203:   (*ksp->converged)(ksp,ksp->its,res_norm,&ksp->reason,ksp->cnvP);
    204:   if (ksp->reason) {
    205:     if (itcount) *itcount = 0;
    206:     return(0);
    207:   }

    209:   /* scale VEC_VV (the initial residual) */
    210:   tmp = 1.0/res_norm; VecScale(VEC_VV(0),tmp);

    214:
    215:   /* MAIN ITERATION LOOP BEGINNING*/
    216:   /* keep iterating until we have converged OR generated the max number
    217:      of directions OR reached the max number of iterations for the method */
    218:   while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
    219:     KSPLogResidualHistory(ksp,res_norm);
    220:     fgmres->it = (loc_it - 1);
    221:     KSPMonitor(ksp,ksp->its,res_norm);

    223:     /* see if more space is needed for work vectors */
    224:     if (fgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
    225:       FGMRESGetNewVectors(ksp,loc_it+1);
    226:       /* (loc_it+1) is passed in as number of the first vector that should
    227:          be allocated */
    228:     }

    230:     /* CHANGE THE PRECONDITIONER? */
    231:     /* ModifyPC is the callback function that can be used to
    232:        change the PC or its attributes before its applied */
    233:     (*fgmres->modifypc)(ksp,ksp->its,loc_it,res_norm,fgmres->modifyctx);
    234:
    235:
    236:     /* apply PRECONDITIONER to direction vector and store with
    237:        preconditioned vectors in prevec */
    238:     KSP_PCApply(ksp,VEC_VV(loc_it),PREVEC(loc_it));
    239:
    240:     PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
    241:     /* Multiply preconditioned vector by operator - put in VEC_VV(loc_it+1) */
    242:     MatMult(Amat,PREVEC(loc_it),VEC_VV(1+loc_it));

    244:
    245:     /* update hessenberg matrix and do Gram-Schmidt - new direction is in
    246:        VEC_VV(1+loc_it)*/
    247:     (*fgmres->orthog)(ksp,loc_it);

    249:     /* new entry in hessenburg is the 2-norm of our new direction */
    250:     VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);
    251:     *HH(loc_it+1,loc_it)   = tt;
    252:     *HES(loc_it+1,loc_it)  = tt;

    254:     /* Happy Breakdown Check */
    255:     hapbnd  = PetscAbsScalar((tt) / *RS(loc_it));
    256:     /* RS(loc_it) contains the res_norm from the last iteration  */
    257:     hapbnd = PetscMin(fgmres->haptol,hapbnd);
    258:     if (tt > hapbnd) {
    259:         tmp = 1.0/tt;
    260:         /* scale new direction by its norm */
    261:         VecScale(VEC_VV(loc_it+1),tmp);
    262:     } else {
    263:         /* This happens when the solution is exactly reached. */
    264:         /* So there is no new direction... */
    265:           VecSet(VEC_TEMP,zero); /* set VEC_TEMP to 0 */
    266:           hapend = PETSC_TRUE;
    267:     }
    268:     /* note that for FGMRES we could get HES(loc_it+1, loc_it)  = 0 and the
    269:        current solution would not be exact if HES was singular.  Note that
    270:        HH non-singular implies that HES is no singular, and HES is guaranteed
    271:        to be nonsingular when PREVECS are linearly independent and A is
    272:        nonsingular (in GMRES, the nonsingularity of A implies the nonsingularity
    273:        of HES). So we should really add a check to verify that HES is nonsingular.*/

    275:
    276:     /* Now apply rotations to new col of hessenberg (and right side of system),
    277:        calculate new rotation, and get new residual norm at the same time*/
    278:     FGMRESUpdateHessenberg(ksp,loc_it,hapend,&res_norm);
    279:     if (ksp->reason) break;

    281:     loc_it++;
    282:     fgmres->it  = (loc_it-1);  /* Add this here in case it has converged */
    283:
    284:     PetscObjectTakeAccess(ksp);
    285:     ksp->its++;
    286:     ksp->rnorm = res_norm;
    287:     PetscObjectGrantAccess(ksp);

    289:     (*ksp->converged)(ksp,ksp->its,res_norm,&ksp->reason,ksp->cnvP);

    291:     /* Catch error in happy breakdown and signal convergence and break from loop */
    292:     if (hapend) {
    293:       if (!ksp->reason) {
    294:         SETERRQ(0,"You reached the happy break down,but convergence was not indicated.");
    295:       }
    296:       break;
    297:     }
    298:   }
    299:   /* END OF ITERATION LOOP */

    301:   KSPLogResidualHistory(ksp,res_norm);

    303:   /*
    304:      Monitor if we know that we will not return for a restart */
    305:   if (ksp->reason || ksp->its >= ksp->max_it) {
    306:     KSPMonitor(ksp,ksp->its,res_norm);
    307:   }

    309:   if (itcount) *itcount    = loc_it;

    311:   /*
    312:     Down here we have to solve for the "best" coefficients of the Krylov
    313:     columns, add the solution values together, and possibly unwind the
    314:     preconditioning from the solution
    315:    */
    316:
    317:   /* Form the solution (or the solution so far) */
    318:   /* Note: must pass in (loc_it-1) for iteration count so that BuildFgmresSoln
    319:      properly navigates */

    321:   BuildFgmresSoln(RS(0),ksp->vec_sol,ksp->vec_sol,ksp,loc_it-1);

    323:   return(0);
    324: }

    326: /*
    327:     KSPSolve_FGMRES - This routine applies the FGMRES method.

    330:    Input Parameter:
    331: .     ksp - the Krylov space object that was set to use fgmres

    333:    Output Parameter:
    334: .     outits - number of iterations used

    336: */

    340: PetscErrorCode KSPSolve_FGMRES(KSP ksp)
    341: {
    343:   PetscInt       cycle_its = 0; /* iterations done in a call to FGMREScycle */
    344:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)ksp->data;
    345:   PetscTruth     diagonalscale;

    348:   PCDiagonalScale(ksp->pc,&diagonalscale);
    349:   if (diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",ksp->type_name);

    351:   PetscObjectTakeAccess(ksp);
    352:   ksp->its = 0;
    353:   PetscObjectGrantAccess(ksp);

    355:   /* Compute the initial (NOT preconditioned) residual */
    356:   if (!ksp->guess_zero) {
    357:     FGMRESResidual(ksp);
    358:   } else { /* guess is 0 so residual is F (which is in ksp->vec_rhs) */
    359:     VecCopy(ksp->vec_rhs,VEC_VV(0));
    360:   }
    361:   /* now the residual is in VEC_VV(0) - which is what
    362:      FGMREScycle expects... */
    363:
    364:   FGMREScycle(&cycle_its,ksp);
    365:   while (!ksp->reason) {
    366:     FGMRESResidual(ksp);
    367:     if (ksp->its >= ksp->max_it) break;
    368:     FGMREScycle(&cycle_its,ksp);
    369:   }
    370:   /* mark lack of convergence */
    371:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;

    373:   return(0);
    374: }

    376: /*

    378:    KSPDestroy_FGMRES - Frees all memory space used by the Krylov method.

    380: */
    383: PetscErrorCode KSPDestroy_FGMRES(KSP ksp)
    384: {
    385:   KSP_FGMRES     *fgmres = (KSP_FGMRES*)ksp->data;
    387:   PetscInt       i;

    390:   /* Free the Hessenberg matrices */
    391:   if (fgmres->hh_origin) {PetscFree(fgmres->hh_origin);}

    393:   /* Free pointers to user variables */
    394:   if (fgmres->vecs) {PetscFree(fgmres->vecs);}
    395:   if (fgmres->prevecs) {PetscFree (fgmres->prevecs);}

    397:   /* free work vectors */
    398:   for (i=0; i < fgmres->nwork_alloc; i++) {
    399:     VecDestroyVecs(fgmres->user_work[i],fgmres->mwork_alloc[i]);
    400:   }
    401:   if (fgmres->user_work)  {PetscFree(fgmres->user_work);}

    403:   for (i=0; i < fgmres->nwork_alloc; i++) {
    404:     VecDestroyVecs(fgmres->prevecs_user_work[i],fgmres->mwork_alloc[i]);
    405:   }
    406:   if (fgmres->prevecs_user_work) {PetscFree(fgmres->prevecs_user_work);}

    408:   if (fgmres->mwork_alloc) {PetscFree(fgmres->mwork_alloc);}
    409:   if (fgmres->nrs) {PetscFree(fgmres->nrs);}
    410:   if (fgmres->sol_temp) {VecDestroy(fgmres->sol_temp);}
    411:   if (fgmres->Rsvd) {PetscFree(fgmres->Rsvd);}
    412:   if (fgmres->Dsvd) {PetscFree(fgmres->Dsvd);}
    413:   if (fgmres->modifydestroy) {
    414:     (*fgmres->modifydestroy)(fgmres->modifyctx);
    415:   }
    416:   PetscFree(fgmres);
    417:   return(0);
    418: }

    420: /*
    421:     BuildFgmresSoln - create the solution from the starting vector and the
    422:                       current iterates.

    424:     Input parameters:
    425:         nrs - work area of size it + 1.
    426:         vguess  - index of initial guess
    427:         vdest - index of result.  Note that vguess may == vdest (replace
    428:                 guess with the solution).
    429:         it - HH upper triangular part is a block of size (it+1) x (it+1)

    431:      This is an internal routine that knows about the FGMRES internals.
    432:  */
    435: static PetscErrorCode BuildFgmresSoln(PetscScalar* nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
    436: {
    437:   PetscScalar    tt,zero = 0.0,one = 1.0;
    439:   PetscInt       ii,k,j;
    440:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)(ksp->data);

    443:   /* Solve for solution vector that minimizes the residual */

    445:   /* If it is < 0, no fgmres steps have been performed */
    446:   if (it < 0) {
    447:     if (vdest != vguess) {
    448:       VecCopy(vguess,vdest);
    449:     }
    450:     return(0);
    451:   }

    453:   /* so fgmres steps HAVE been performed */

    455:   /* solve the upper triangular system - RS is the right side and HH is
    456:      the upper triangular matrix  - put soln in nrs */
    457:   nrs[it] = *RS(it) / *HH(it,it);
    458:   for (ii=1; ii<=it; ii++) {
    459:     k   = it - ii;
    460:     tt  = *RS(k);
    461:     for (j=k+1; j<=it; j++) tt  = tt - *HH(k,j) * nrs[j];
    462:     nrs[k]   = tt / *HH(k,k);
    463:   }

    465:   /* Accumulate the correction to the soln of the preconditioned prob. in
    466:      VEC_TEMP - note that we use the preconditioned vectors  */
    467:   VecSet(VEC_TEMP,zero); /* set VEC_TEMP components to 0 */
    468:   VecMAXPY(VEC_TEMP,it+1,nrs,&PREVEC(0));

    470:   /* put updated solution into vdest.*/
    471:   if (vdest != vguess) {
    472:     VecCopy(VEC_TEMP,vdest);
    473:     VecAXPY(vdest,one,vguess);
    474:   } else  {/* replace guess with solution */
    475:     VecAXPY(vdest,one,VEC_TEMP);
    476:   }
    477:   return(0);
    478: }

    480: /*

    482:     FGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
    483:                             Return new residual.

    485:     input parameters:

    487: .        ksp -    Krylov space object
    488: .         it  -    plane rotations are applied to the (it+1)th column of the
    489:                   modified hessenberg (i.e. HH(:,it))
    490: .        hapend - PETSC_FALSE not happy breakdown ending.

    492:     output parameters:
    493: .        res - the new residual
    494:
    495:  */
    498: static PetscErrorCode FGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscTruth hapend,PetscReal *res)
    499: {
    500:   PetscScalar   *hh,*cc,*ss,tt;
    501:   PetscInt      j;
    502:   KSP_FGMRES    *fgmres = (KSP_FGMRES *)(ksp->data);

    505:   hh  = HH(0,it);  /* pointer to beginning of column to update - so
    506:                       incrementing hh "steps down" the (it+1)th col of HH*/
    507:   cc  = CC(0);     /* beginning of cosine rotations */
    508:   ss  = SS(0);     /* beginning of sine rotations */

    510:   /* Apply all the previously computed plane rotations to the new column
    511:      of the Hessenberg matrix */
    512:   /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta),
    513:      and some refs have [c   s ; -conj(s)  c] (don't be confused!) */

    515:   for (j=1; j<=it; j++) {
    516:     tt  = *hh;
    517: #if defined(PETSC_USE_COMPLEX)
    518:     *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
    519: #else
    520:     *hh = *cc * tt + *ss * *(hh+1);
    521: #endif
    522:     hh++;
    523:     *hh = *cc++ * *hh - (*ss++ * tt);
    524:     /* hh, cc, and ss have all been incremented one by end of loop */
    525:   }

    527:   /*
    528:     compute the new plane rotation, and apply it to:
    529:      1) the right-hand-side of the Hessenberg system (RS)
    530:         note: it affects RS(it) and RS(it+1)
    531:      2) the new column of the Hessenberg matrix
    532:         note: it affects HH(it,it) which is currently pointed to
    533:         by hh and HH(it+1, it) (*(hh+1))
    534:     thus obtaining the updated value of the residual...
    535:   */

    537:   /* compute new plane rotation */

    539:   if (!hapend) {
    540: #if defined(PETSC_USE_COMPLEX)
    541:     tt        = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
    542: #else
    543:     tt        = PetscSqrtScalar(*hh * *hh + *(hh+1) * *(hh+1));
    544: #endif
    545:     if (tt == 0.0) {
    546:       ksp->reason = KSP_DIVERGED_NULL;
    547:       return(0);
    548:     }

    550:     *cc       = *hh / tt;   /* new cosine value */
    551:     *ss       = *(hh+1) / tt;  /* new sine value */

    553:     /* apply to 1) and 2) */
    554:     *RS(it+1) = - (*ss * *RS(it));
    555: #if defined(PETSC_USE_COMPLEX)
    556:     *RS(it)   = PetscConj(*cc) * *RS(it);
    557:     *hh       = PetscConj(*cc) * *hh + *ss * *(hh+1);
    558: #else
    559:     *RS(it)   = *cc * *RS(it);
    560:     *hh       = *cc * *hh + *ss * *(hh+1);
    561: #endif

    563:     /* residual is the last element (it+1) of right-hand side! */
    564:     *res      = PetscAbsScalar(*RS(it+1));

    566:   } else { /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
    567:             another rotation matrix (so RH doesn't change).  The new residual is
    568:             always the new sine term times the residual from last time (RS(it)),
    569:             but now the new sine rotation would be zero...so the residual should
    570:             be zero...so we will multiply "zero" by the last residual.  This might
    571:             not be exactly what we want to do here -could just return "zero". */
    572:
    573:     *res = 0.0;
    574:   }
    575:   return(0);
    576: }

    578: /*

    580:    FGMRESGetNewVectors - This routine allocates more work vectors, starting from
    581:                          VEC_VV(it), and more preconditioned work vectors, starting
    582:                          from PREVEC(i).

    584: */
    587: static PetscErrorCode FGMRESGetNewVectors(KSP ksp,PetscInt it)
    588: {
    589:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)ksp->data;
    590:   PetscInt       nwork = fgmres->nwork_alloc; /* number of work vector chunks allocated */
    591:   PetscInt       nalloc;                      /* number to allocate */
    593:   PetscInt       k;
    594:
    596:   nalloc = fgmres->delta_allocate; /* number of vectors to allocate
    597:                                       in a single chunk */

    599:   /* Adjust the number to allocate to make sure that we don't exceed the
    600:      number of available slots (fgmres->vecs_allocated)*/
    601:   if (it + VEC_OFFSET + nalloc >= fgmres->vecs_allocated){
    602:     nalloc = fgmres->vecs_allocated - it - VEC_OFFSET;
    603:   }
    604:   if (!nalloc) return(0);

    606:   fgmres->vv_allocated += nalloc; /* vv_allocated is the number of vectors allocated */

    608:   /* work vectors */
    609:   KSPGetVecs(ksp,nalloc,&fgmres->user_work[nwork]);
    610:   PetscLogObjectParents(ksp,nalloc,fgmres->user_work[nwork]);
    611:   for (k=0; k < nalloc; k++) {
    612:     fgmres->vecs[it+VEC_OFFSET+k] = fgmres->user_work[nwork][k];
    613:   }
    614:   /* specify size of chunk allocated */
    615:   fgmres->mwork_alloc[nwork] = nalloc;

    617:   /* preconditioned vectors */
    618:   KSPGetVecs(ksp,nalloc,&fgmres->prevecs_user_work[nwork]);
    619:   PetscLogObjectParents(ksp,nalloc,fgmres->prevecs_user_work[nwork]);
    620:   for (k=0; k < nalloc; k++) {
    621:     fgmres->prevecs[it+VEC_OFFSET+k] = fgmres->prevecs_user_work[nwork][k];
    622:   }

    624:   /* increment the number of work vector chunks */
    625:   fgmres->nwork_alloc++;
    626:   return(0);
    627: }

    629: /*

    631:    KSPBuildSolution_FGMRES

    633:      Input Parameter:
    634: .     ksp - the Krylov space object
    635: .     ptr-

    637:    Output Parameter:
    638: .     result - the solution

    640:    Note: this calls BuildFgmresSoln - the same function that FGMREScycle
    641:    calls directly.

    643: */
    646: PetscErrorCode KSPBuildSolution_FGMRES(KSP ksp,Vec ptr,Vec *result)
    647: {
    648:   KSP_FGMRES     *fgmres = (KSP_FGMRES *)ksp->data;

    652:   if (!ptr) {
    653:     if (!fgmres->sol_temp) {
    654:       VecDuplicate(ksp->vec_sol,&fgmres->sol_temp);
    655:       PetscLogObjectParent(ksp,fgmres->sol_temp);
    656:     }
    657:     ptr = fgmres->sol_temp;
    658:   }
    659:   if (!fgmres->nrs) {
    660:     /* allocate the work area */
    661:     PetscMalloc(fgmres->max_k*sizeof(PetscScalar),&fgmres->nrs);
    662:     PetscLogObjectMemory(ksp,fgmres->max_k*sizeof(PetscScalar));
    663:   }
    664:
    665:   BuildFgmresSoln(fgmres->nrs,ksp->vec_sol,ptr,ksp,fgmres->it);
    666:   *result = ptr;
    667:
    668:   return(0);
    669: }

    675: PetscErrorCode KSPSetFromOptions_FGMRES(KSP ksp)
    676: {
    678:   PetscTruth     flg;

    681:   KSPSetFromOptions_GMRES(ksp);
    682:   PetscOptionsHead("KSP flexible GMRES Options");
    683:     PetscOptionsTruthGroupBegin("-ksp_fgmres_modifypcnochange","do not vary the preconditioner","KSPFGMRESSetModifyPC",&flg);
    684:     if (flg) {KSPFGMRESSetModifyPC(ksp,KSPFGMRESModifyPCNoChange,0,0);}
    685:     PetscOptionsTruthGroupEnd("-ksp_fgmres_modifypcksp","vary the KSP based preconditioner","KSPFGMRESSetModifyPC",&flg);
    686:     if (flg) {KSPFGMRESSetModifyPC(ksp,KSPFGMRESModifyPCKSP,0,0);}
    687:   PetscOptionsTail();
    688:   return(0);
    689: }

    691: EXTERN PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP,PetscReal *,PetscReal *);
    692: EXTERN PetscErrorCode KSPComputeEigenvalues_GMRES(KSP,PetscInt,PetscReal *,PetscReal *,PetscInt *);

    695: typedef PetscErrorCode (*FCN2)(void*);
    699: PetscErrorCode PETSCKSP_DLLEXPORT KSPFGMRESSetModifyPC_FGMRES(KSP ksp,FCN1 fcn,void *ctx,FCN2 d)
    700: {
    703:   ((KSP_FGMRES *)ksp->data)->modifypc      = fcn;
    704:   ((KSP_FGMRES *)ksp->data)->modifydestroy = d;
    705:   ((KSP_FGMRES *)ksp->data)->modifyctx     = ctx;
    706:   return(0);
    707: }

    711: EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESSetPreAllocateVectors_GMRES(KSP);
    712: EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESSetRestart_GMRES(KSP,PetscInt);
    713: EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESSetOrthogonalization_GMRES(KSP,PetscErrorCode (*)(KSP,PetscInt));

    718: PetscErrorCode KSPDestroy_FGMRES_Internal(KSP ksp)
    719: {
    720:   KSP_FGMRES     *gmres = (KSP_FGMRES*)ksp->data;
    722:   PetscInt       i;

    725:   /* Free the Hessenberg matrix */
    726:   if (gmres->hh_origin) {
    727:     PetscFree(gmres->hh_origin);
    728:     gmres->hh_origin = 0;
    729:   }

    731:   /* Free the pointer to user variables */
    732:   if (gmres->vecs) {
    733:     PetscFree(gmres->vecs);
    734:     gmres->vecs = 0;
    735:   }
    736:   if (gmres->prevecs) {
    737:     PetscFree (gmres->prevecs);
    738:     gmres->prevecs = 0;
    739:   }

    741:   /* free work vectors */
    742:   for (i=0; i<gmres->nwork_alloc; i++) {
    743:     VecDestroyVecs(gmres->user_work[i],gmres->mwork_alloc[i]);
    744:     VecDestroyVecs(gmres->prevecs_user_work[i],gmres->mwork_alloc[i]);
    745:   }
    746:   if (gmres->user_work)  {
    747:     PetscFree(gmres->user_work);
    748:     gmres->user_work = 0;
    749:   }
    750:   if (gmres->prevecs_user_work) {
    751:     PetscFree(gmres->prevecs_user_work);
    752:     gmres->prevecs_user_work = 0;
    753:   }
    754:   if (gmres->mwork_alloc) {
    755:     PetscFree(gmres->mwork_alloc);
    756:     gmres->mwork_alloc = 0;
    757:   }
    758:   if (gmres->nrs) {
    759:     PetscFree(gmres->nrs);
    760:     gmres->nrs = 0;
    761:   }
    762:   if (gmres->sol_temp) {
    763:     VecDestroy(gmres->sol_temp);
    764:     gmres->sol_temp = 0;
    765:   }
    766:   if (gmres->Rsvd) {
    767:     PetscFree(gmres->Rsvd);
    768:     gmres->Rsvd = 0;
    769:   }
    770:   if (gmres->Dsvd) {
    771:     PetscFree(gmres->Dsvd);
    772:     gmres->Dsvd = 0;
    773:   }
    774:   if (gmres->modifydestroy) {
    775:     (*gmres->modifydestroy)(gmres->modifyctx);
    776:   }

    778:   gmres->vv_allocated   = 0;
    779:   gmres->vecs_allocated = 0;
    780:   gmres->sol_temp       = 0;
    781:   gmres->nwork_alloc    = 0;
    782:   return(0);
    783: }

    788: PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESSetRestart_FGMRES(KSP ksp,PetscInt max_k)
    789: {
    790:   KSP_FGMRES     *gmres = (KSP_FGMRES *)ksp->data;

    794:   if (max_k < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
    795:   if (!ksp->setupcalled) {
    796:     gmres->max_k = max_k;
    797:   } else if (gmres->max_k != max_k) {
    798:      gmres->max_k = max_k;
    799:      ksp->setupcalled = 0;
    800:      /* free the data structures, then create them again */
    801:      KSPDestroy_FGMRES_Internal(ksp);
    802:   }
    803:   return(0);
    804: }

    808: EXTERN PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESSetCGSRefinementType_GMRES(KSP,KSPGMRESCGSRefinementType);

    811: /*MC
    812:      KSPFGMRES - Implements the Flexible Generalized Minimal Residual method.
    813:                 developed by Saad with restart

    816:    Options Database Keys:
    817: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
    818: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
    819: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
    820:                              vectors are allocated as needed)
    821: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
    822: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
    823: .   -ksp_gmres_cgs_refinement_type <never,ifneeded,always> - determine if iterative refinement is used to increase the
    824:                                    stability of the classical Gram-Schmidt  orthogonalization.
    825: .   -ksp_gmres_krylov_monitor - plot the Krylov space generated
    826: .   -ksp_fgmres_modifypcnochange - do not change the preconditioner between iterations
    827: -   -ksp_fgmres_modifypcksp - modify the preconditioner using KSPFGMRESModifyPCKSP()

    829:    Level: beginner

    831:     Notes: See KSPFGMRESSetModifyPC() for how to vary the preconditioner between iterations
    832:            This object is subclassed off of KSPGMRES

    834: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPLGMRES,
    835:            KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization()
    836:            KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
    837:            KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESKrylovMonitor(), KSPFGMRESSetModifyPC(),
    838:            KSPFGMRESModifyPCKSP()

    840: M*/

    845: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_FGMRES(KSP ksp)
    846: {
    847:   KSP_FGMRES     *fgmres;

    851:   PetscNew(KSP_FGMRES,&fgmres);
    852:   PetscLogObjectMemory(ksp,sizeof(KSP_FGMRES));
    853:   ksp->data                              = (void*)fgmres;
    854:   ksp->ops->buildsolution                = KSPBuildSolution_FGMRES;

    856:   ksp->ops->setup                        = KSPSetUp_FGMRES;
    857:   ksp->ops->solve                        = KSPSolve_FGMRES;
    858:   ksp->ops->destroy                      = KSPDestroy_FGMRES;
    859:   ksp->ops->view                         = KSPView_GMRES;
    860:   ksp->ops->setfromoptions               = KSPSetFromOptions_FGMRES;
    861:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
    862:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

    864:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",
    865:                                     "KSPGMRESSetPreAllocateVectors_GMRES",
    866:                                      KSPGMRESSetPreAllocateVectors_GMRES);
    867:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",
    868:                                     "KSPGMRESSetOrthogonalization_GMRES",
    869:                                      KSPGMRESSetOrthogonalization_GMRES);
    870:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGMRESSetRestart_C",
    871:                                     "KSPGMRESSetRestart_FGMRES",
    872:                                      KSPGMRESSetRestart_FGMRES);
    873:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPFGMRESSetModifyPC_C",
    874:                                     "KSPFGMRESSetModifyPC_FGMRES",
    875:                                      KSPFGMRESSetModifyPC_FGMRES);
    876:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",
    877:                                     "KSPGMRESSetCGSRefinementType_GMRES",
    878:                                      KSPGMRESSetCGSRefinementType_GMRES);

    881:   fgmres->haptol              = 1.0e-30;
    882:   fgmres->q_preallocate       = 0;
    883:   fgmres->delta_allocate      = FGMRES_DELTA_DIRECTIONS;
    884:   fgmres->orthog              = KSPGMRESClassicalGramSchmidtOrthogonalization;
    885:   fgmres->nrs                 = 0;
    886:   fgmres->sol_temp            = 0;
    887:   fgmres->max_k               = FGMRES_DEFAULT_MAXK;
    888:   fgmres->Rsvd                = 0;
    889:   fgmres->modifypc            = KSPFGMRESModifyPCNoChange;
    890:   fgmres->modifyctx           = PETSC_NULL;
    891:   fgmres->modifydestroy       = PETSC_NULL;
    892:   fgmres->cgstype             = KSP_GMRES_CGS_REFINE_NEVER;
    893:   /*
    894:         This is not great since it changes this without explicit request from the user
    895:      but there is no left preconditioning in the FGMRES
    896:   */
    897:   PetscLogInfo((ksp,"KSPCreate_FGMRES: WARNING! Setting PC_SIDE for FGMRES to right!\n"));
    898:   ksp->pc_side                = PC_RIGHT;

    900:   return(0);
    901: }