2: #include <../src/ksp/ksp/impls/gmres/lgmres/lgmresimpl.h> 4: #define LGMRES_DELTA_DIRECTIONS 10 5: #define LGMRES_DEFAULT_MAXK 30 6: #define LGMRES_DEFAULT_AUGDIM 2 /*default number of augmentation vectors */ 7: static PetscErrorCode KSPLGMRESGetNewVectors(KSP,PetscInt);
8: static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
9: static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
11: PetscErrorCode KSPLGMRESSetAugDim(KSP ksp, PetscInt dim) 12: {
16: PetscTryMethod((ksp),"KSPLGMRESSetAugDim_C",(KSP,PetscInt),(ksp,dim));
17: return(0);
18: }
20: PetscErrorCode KSPLGMRESSetConstant(KSP ksp) 21: {
25: PetscTryMethod((ksp),"KSPLGMRESSetConstant_C",(KSP),(ksp));
26: return(0);
27: }
29: /*
30: KSPSetUp_LGMRES - Sets up the workspace needed by lgmres.
32: This is called once, usually automatically by KSPSolve() or KSPSetUp(),
33: but can be called directly by KSPSetUp().
35: */
36: PetscErrorCode KSPSetUp_LGMRES(KSP ksp) 37: {
39: PetscInt max_k,k, aug_dim;
40: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
43: max_k = lgmres->max_k;
44: aug_dim = lgmres->aug_dim;
45: KSPSetUp_GMRES(ksp);
47: /* need array of pointers to augvecs*/
48: PetscMalloc1(2*aug_dim + AUG_OFFSET,&lgmres->augvecs);
50: lgmres->aug_vecs_allocated = 2 *aug_dim + AUG_OFFSET;
52: PetscMalloc1(2*aug_dim + AUG_OFFSET,&lgmres->augvecs_user_work);
53: PetscMalloc1(aug_dim,&lgmres->aug_order);
54: PetscLogObjectMemory((PetscObject)ksp,(aug_dim)*(4*sizeof(void*) + sizeof(PetscInt)) + AUG_OFFSET*2*sizeof(void*));
56: /* for now we will preallocate the augvecs - because aug_dim << restart
57: ... also keep in mind that we need to keep augvecs from cycle to cycle*/
58: lgmres->aug_vv_allocated = 2* aug_dim + AUG_OFFSET;
59: lgmres->augwork_alloc = 2* aug_dim + AUG_OFFSET;
61: KSPCreateVecs(ksp,lgmres->aug_vv_allocated,&lgmres->augvecs_user_work[0],0,NULL);
62: PetscMalloc1(max_k+1,&lgmres->hwork);
63: PetscLogObjectParents(ksp,lgmres->aug_vv_allocated,lgmres->augvecs_user_work[0]);
64: for (k=0; k<lgmres->aug_vv_allocated; k++) {
65: lgmres->augvecs[k] = lgmres->augvecs_user_work[0][k];
66: }
67: return(0);
68: }
71: /*
73: KSPLGMRESCycle - Run lgmres, possibly with restart. Return residual
74: history if requested.
76: input parameters:
77: . lgmres - structure containing parameters and work areas
79: output parameters:
80: . nres - residuals (from preconditioned system) at each step.
81: If restarting, consider passing nres+it. If null,
82: ignored
83: . itcount - number of iterations used. nres[0] to nres[itcount]
84: are defined. If null, ignored. If null, ignored.
85: . converged - 0 if not converged
88: Notes:
89: On entry, the value in vector VEC_VV(0) should be
90: the initial residual.
93: */
94: PetscErrorCode KSPLGMRESCycle(PetscInt *itcount,KSP ksp) 95: {
96: KSP_LGMRES *lgmres = (KSP_LGMRES*)(ksp->data);
97: PetscReal res_norm, res;
98: PetscReal hapbnd, tt;
99: PetscScalar tmp;
100: PetscBool hapend = PETSC_FALSE; /* indicates happy breakdown ending */
102: PetscInt loc_it; /* local count of # of dir. in Krylov space */
103: PetscInt max_k = lgmres->max_k; /* max approx space size */
104: PetscInt max_it = ksp->max_it; /* max # of overall iterations for the method */
106: /* LGMRES_MOD - new variables*/
107: PetscInt aug_dim = lgmres->aug_dim;
108: PetscInt spot = 0;
109: PetscInt order = 0;
110: PetscInt it_arnoldi; /* number of arnoldi steps to take */
111: PetscInt it_total; /* total number of its to take (=approx space size)*/
112: PetscInt ii, jj;
113: PetscReal tmp_norm;
114: PetscScalar inv_tmp_norm;
115: PetscScalar *avec;
118: /* Number of pseudo iterations since last restart is the number
119: of prestart directions */
120: loc_it = 0;
122: /* LGMRES_MOD: determine number of arnoldi steps to take */
123: /* if approx_constant then we keep the space the same size even if
124: we don't have the full number of aug vectors yet*/
125: if (lgmres->approx_constant) it_arnoldi = max_k - lgmres->aug_ct;
126: else it_arnoldi = max_k - aug_dim;
128: it_total = it_arnoldi + lgmres->aug_ct;
130: /* initial residual is in VEC_VV(0) - compute its norm*/
131: VecNorm(VEC_VV(0),NORM_2,&res_norm);
132: KSPCheckNorm(ksp,res_norm);
133: res = res_norm;
135: /* first entry in right-hand-side of hessenberg system is just
136: the initial residual norm */
137: *GRS(0) = res_norm;
139: /* check for the convergence */
140: if (!res) {
141: if (itcount) *itcount = 0;
142: ksp->reason = KSP_CONVERGED_ATOL;
143: PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
144: return(0);
145: }
147: /* scale VEC_VV (the initial residual) */
148: tmp = 1.0/res_norm; VecScale(VEC_VV(0),tmp);
150: ksp->rnorm = res;
153: /* note: (lgmres->it) is always set one less than (loc_it) It is used in
154: KSPBUILDSolution_LGMRES, where it is passed to KSPLGMRESBuildSoln.
155: Note that when KSPLGMRESBuildSoln is called from this function,
156: (loc_it -1) is passed, so the two are equivalent */
157: lgmres->it = (loc_it - 1);
160: /* MAIN ITERATION LOOP BEGINNING*/
163: /* keep iterating until we have converged OR generated the max number
164: of directions OR reached the max number of iterations for the method */
165: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
167: while (!ksp->reason && loc_it < it_total && ksp->its < max_it) { /* LGMRES_MOD: changed to it_total */
168: KSPLogResidualHistory(ksp,res);
169: lgmres->it = (loc_it - 1);
170: KSPMonitor(ksp,ksp->its,res);
172: /* see if more space is needed for work vectors */
173: if (lgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
174: KSPLGMRESGetNewVectors(ksp,loc_it+1);
175: /* (loc_it+1) is passed in as number of the first vector that should
176: be allocated */
177: }
179: /*LGMRES_MOD: decide whether this is an arnoldi step or an aug step */
180: if (loc_it < it_arnoldi) { /* Arnoldi */
181: KSP_PCApplyBAorAB(ksp,VEC_VV(loc_it),VEC_VV(1+loc_it),VEC_TEMP_MATOP);
182: } else { /*aug step */
183: order = loc_it - it_arnoldi + 1; /* which aug step */
184: for (ii=0; ii<aug_dim; ii++) {
185: if (lgmres->aug_order[ii] == order) {
186: spot = ii;
187: break; /* must have this because there will be duplicates before aug_ct = aug_dim */
188: }
189: }
191: VecCopy(A_AUGVEC(spot), VEC_VV(1+loc_it));
192: /*note: an alternate implementation choice would be to only save the AUGVECS and
193: not A_AUGVEC and then apply the PC here to the augvec */
194: }
196: /* update hessenberg matrix and do Gram-Schmidt - new direction is in
197: VEC_VV(1+loc_it)*/
198: (*lgmres->orthog)(ksp,loc_it);
200: /* new entry in hessenburg is the 2-norm of our new direction */
201: VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);
203: *HH(loc_it+1,loc_it) = tt;
204: *HES(loc_it+1,loc_it) = tt;
207: /* check for the happy breakdown */
208: hapbnd = PetscAbsScalar(tt / *GRS(loc_it)); /* GRS(loc_it) contains the res_norm from the last iteration */
209: if (hapbnd > lgmres->haptol) hapbnd = lgmres->haptol;
210: if (tt > hapbnd) {
211: tmp = 1.0/tt;
212: VecScale(VEC_VV(loc_it+1),tmp); /* scale new direction by its norm */
213: } else {
214: PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %g tt = %g\n",(double)hapbnd,(double)tt);
215: hapend = PETSC_TRUE;
216: }
218: /* Now apply rotations to new col of hessenberg (and right side of system),
219: calculate new rotation, and get new residual norm at the same time*/
220: KSPLGMRESUpdateHessenberg(ksp,loc_it,hapend,&res);
221: if (ksp->reason) break;
223: loc_it++;
224: lgmres->it = (loc_it-1); /* Add this here in case it has converged */
226: PetscObjectSAWsTakeAccess((PetscObject)ksp);
227: ksp->its++;
228: ksp->rnorm = res;
229: PetscObjectSAWsGrantAccess((PetscObject)ksp);
231: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
233: /* Catch error in happy breakdown and signal convergence and break from loop */
234: if (hapend) {
235: if (!ksp->reason) {
236: if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %g",(double)res);
237: else {
238: ksp->reason = KSP_DIVERGED_BREAKDOWN;
239: break;
240: }
241: }
242: }
243: }
244: /* END OF ITERATION LOOP */
245: KSPLogResidualHistory(ksp,res);
247: /* Monitor if we know that we will not return for a restart */
248: if (ksp->reason || ksp->its >= max_it) {
249: KSPMonitor(ksp, ksp->its, res);
250: }
252: if (itcount) *itcount = loc_it;
254: /*
255: Down here we have to solve for the "best" coefficients of the Krylov
256: columns, add the solution values together, and possibly unwind the
257: preconditioning from the solution
258: */
260: /* Form the solution (or the solution so far) */
261: /* Note: must pass in (loc_it-1) for iteration count so that KSPLGMRESBuildSoln
262: properly navigates */
264: KSPLGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,loc_it-1);
267: /* LGMRES_MOD collect aug vector and A*augvector for future restarts -
268: only if we will be restarting (i.e. this cycle performed it_total
269: iterations) */
270: if (!ksp->reason && ksp->its < max_it && aug_dim > 0) {
272: /*AUG_TEMP contains the new augmentation vector (assigned in KSPLGMRESBuildSoln) */
273: if (!lgmres->aug_ct) {
274: spot = 0;
275: lgmres->aug_ct++;
276: } else if (lgmres->aug_ct < aug_dim) {
277: spot = lgmres->aug_ct;
278: lgmres->aug_ct++;
279: } else { /* truncate */
280: for (ii=0; ii<aug_dim; ii++) {
281: if (lgmres->aug_order[ii] == aug_dim) spot = ii;
282: }
283: }
287: VecCopy(AUG_TEMP, AUGVEC(spot));
288: /*need to normalize */
289: VecNorm(AUGVEC(spot), NORM_2, &tmp_norm);
291: inv_tmp_norm = 1.0/tmp_norm;
293: VecScale(AUGVEC(spot),inv_tmp_norm);
295: /*set new aug vector to order 1 - move all others back one */
296: for (ii=0; ii < aug_dim; ii++) AUG_ORDER(ii)++;
297: AUG_ORDER(spot) = 1;
299: /*now add the A*aug vector to A_AUGVEC(spot) - this is independ. of preconditioning type*/
300: /* want V*H*y - y is in GRS, V is in VEC_VV and H is in HES */
303: /* first do H+*y */
304: avec = lgmres->hwork;
305: PetscMemzero(avec,(it_total+1)*sizeof(*avec));
306: for (ii=0; ii < it_total + 1; ii++) {
307: for (jj=0; jj <= ii+1 && jj < it_total+1; jj++) {
308: avec[jj] += *HES(jj ,ii) * *GRS(ii);
309: }
310: }
312: /*now multiply result by V+ */
313: VecSet(VEC_TEMP,0.0);
314: VecMAXPY(VEC_TEMP, it_total+1, avec, &VEC_VV(0)); /*answer is in VEC_TEMP*/
316: /*copy answer to aug location and scale*/
317: VecCopy(VEC_TEMP, A_AUGVEC(spot));
318: VecScale(A_AUGVEC(spot),inv_tmp_norm);
319: }
320: return(0);
321: }
323: /*
324: KSPSolve_LGMRES - This routine applies the LGMRES method.
327: Input Parameter:
328: . ksp - the Krylov space object that was set to use lgmres
330: Output Parameter:
331: . outits - number of iterations used
333: */
335: PetscErrorCode KSPSolve_LGMRES(KSP ksp)336: {
338: PetscInt cycle_its; /* iterations done in a call to KSPLGMRESCycle */
339: PetscInt itcount; /* running total of iterations, incl. those in restarts */
340: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
341: PetscBool guess_zero = ksp->guess_zero;
342: PetscInt ii; /*LGMRES_MOD variable */
345: if (ksp->calc_sings && !lgmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
347: PetscObjectSAWsTakeAccess((PetscObject)ksp);
349: ksp->its = 0;
350: lgmres->aug_ct = 0;
351: lgmres->matvecs = 0;
353: PetscObjectSAWsGrantAccess((PetscObject)ksp);
355: /* initialize */
356: itcount = 0;
357: ksp->reason = KSP_CONVERGED_ITERATING;
358: /*LGMRES_MOD*/
359: for (ii=0; ii<lgmres->aug_dim; ii++) lgmres->aug_order[ii] = 0;
361: while (!ksp->reason) {
362: /* calc residual - puts in VEC_VV(0) */
363: KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
364: KSPLGMRESCycle(&cycle_its,ksp);
365: itcount += cycle_its;
366: if (itcount >= ksp->max_it) {
367: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
368: break;
369: }
370: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
371: }
372: ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
373: return(0);
374: }
376: /*
378: KSPDestroy_LGMRES - Frees all memory space used by the Krylov method.
380: */
381: PetscErrorCode KSPDestroy_LGMRES(KSP ksp)382: {
383: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
387: PetscFree(lgmres->augvecs);
388: if (lgmres->augwork_alloc) {
389: VecDestroyVecs(lgmres->augwork_alloc,&lgmres->augvecs_user_work[0]);
390: }
391: PetscFree(lgmres->augvecs_user_work);
392: PetscFree(lgmres->aug_order);
393: PetscFree(lgmres->hwork);
394: KSPDestroy_GMRES(ksp);
395: return(0);
396: }
398: /*
399: KSPLGMRESBuildSoln - create the solution from the starting vector and the
400: current iterates.
402: Input parameters:
403: nrs - work area of size it + 1.
404: vguess - index of initial guess
405: vdest - index of result. Note that vguess may == vdest (replace
406: guess with the solution).
407: it - HH upper triangular part is a block of size (it+1) x (it+1)
409: This is an internal routine that knows about the LGMRES internals.
410: */
411: static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar *nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)412: {
413: PetscScalar tt;
415: PetscInt ii,k,j;
416: KSP_LGMRES *lgmres = (KSP_LGMRES*)(ksp->data);
417: /*LGMRES_MOD */
418: PetscInt it_arnoldi, it_aug;
419: PetscInt jj, spot = 0;
422: /* Solve for solution vector that minimizes the residual */
424: /* If it is < 0, no lgmres steps have been performed */
425: if (it < 0) {
426: VecCopy(vguess,vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
427: return(0);
428: }
430: /* so (it+1) lgmres steps HAVE been performed */
432: /* LGMRES_MOD - determine if we need to use augvecs for the soln - do not assume that
433: this is called after the total its allowed for an approx space */
434: if (lgmres->approx_constant) {
435: it_arnoldi = lgmres->max_k - lgmres->aug_ct;
436: } else {
437: it_arnoldi = lgmres->max_k - lgmres->aug_dim;
438: }
439: if (it_arnoldi >= it +1) {
440: it_aug = 0;
441: it_arnoldi = it+1;
442: } else {
443: it_aug = (it + 1) - it_arnoldi;
444: }
446: /* now it_arnoldi indicates the number of matvecs that took place */
447: lgmres->matvecs += it_arnoldi;
450: /* solve the upper triangular system - GRS is the right side and HH is
451: the upper triangular matrix - put soln in nrs */
452: if (*HH(it,it) == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"HH(it,it) is identically zero; it = %D GRS(it) = %g",it,(double)PetscAbsScalar(*GRS(it)));
453: if (*HH(it,it) != 0.0) {
454: nrs[it] = *GRS(it) / *HH(it,it);
455: } else {
456: nrs[it] = 0.0;
457: }
459: for (ii=1; ii<=it; ii++) {
460: k = it - ii;
461: tt = *GRS(k);
462: for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
463: nrs[k] = tt / *HH(k,k);
464: }
466: /* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */
467: VecSet(VEC_TEMP,0.0); /* set VEC_TEMP components to 0 */
469: /*LGMRES_MOD - if augmenting has happened we need to form the solution
470: using the augvecs */
471: if (!it_aug) { /* all its are from arnoldi */
472: VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
473: } else { /*use aug vecs */
474: /*first do regular krylov directions */
475: VecMAXPY(VEC_TEMP,it_arnoldi,nrs,&VEC_VV(0));
476: /*now add augmented portions - add contribution of aug vectors one at a time*/
479: for (ii=0; ii<it_aug; ii++) {
480: for (jj=0; jj<lgmres->aug_dim; jj++) {
481: if (lgmres->aug_order[jj] == (ii+1)) {
482: spot = jj;
483: break; /* must have this because there will be duplicates before aug_ct = aug_dim */
484: }
485: }
486: VecAXPY(VEC_TEMP,nrs[it_arnoldi+ii],AUGVEC(spot));
487: }
488: }
489: /* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
490: preconditioner is "unwound" from right-precondtioning*/
491: VecCopy(VEC_TEMP, AUG_TEMP);
493: KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
495: /* add solution to previous solution */
496: /* put updated solution into vdest.*/
497: VecCopy(vguess,vdest);
498: VecAXPY(vdest,1.0,VEC_TEMP);
499: return(0);
500: }
502: /*
504: KSPLGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
505: Return new residual.
507: input parameters:
509: . ksp - Krylov space object
510: . it - plane rotations are applied to the (it+1)th column of the
511: modified hessenberg (i.e. HH(:,it))
512: . hapend - PETSC_FALSE not happy breakdown ending.
514: output parameters:
515: . res - the new residual
517: */
518: static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)519: {
520: PetscScalar *hh,*cc,*ss,tt;
521: PetscInt j;
522: KSP_LGMRES *lgmres = (KSP_LGMRES*)(ksp->data);
525: hh = HH(0,it); /* pointer to beginning of column to update - so
526: incrementing hh "steps down" the (it+1)th col of HH*/
527: cc = CC(0); /* beginning of cosine rotations */
528: ss = SS(0); /* beginning of sine rotations */
530: /* Apply all the previously computed plane rotations to the new column
531: of the Hessenberg matrix */
532: /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta) */
534: for (j=1; j<=it; j++) {
535: tt = *hh;
536: *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
537: hh++;
538: *hh = *cc++ * *hh - (*ss++ * tt);
539: /* hh, cc, and ss have all been incremented one by end of loop */
540: }
542: /*
543: compute the new plane rotation, and apply it to:
544: 1) the right-hand-side of the Hessenberg system (GRS)
545: note: it affects GRS(it) and GRS(it+1)
546: 2) the new column of the Hessenberg matrix
547: note: it affects HH(it,it) which is currently pointed to
548: by hh and HH(it+1, it) (*(hh+1))
549: thus obtaining the updated value of the residual...
550: */
552: /* compute new plane rotation */
554: if (!hapend) {
555: tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
556: if (tt == 0.0) {
557: ksp->reason = KSP_DIVERGED_NULL;
558: return(0);
559: }
560: *cc = *hh / tt; /* new cosine value */
561: *ss = *(hh+1) / tt; /* new sine value */
563: /* apply to 1) and 2) */
564: *GRS(it+1) = -(*ss * *GRS(it));
565: *GRS(it) = PetscConj(*cc) * *GRS(it);
566: *hh = PetscConj(*cc) * *hh + *ss * *(hh+1);
568: /* residual is the last element (it+1) of right-hand side! */
569: *res = PetscAbsScalar(*GRS(it+1));
571: } else { /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
572: another rotation matrix (so RH doesn't change). The new residual is
573: always the new sine term times the residual from last time (GRS(it)),
574: but now the new sine rotation would be zero...so the residual should
575: be zero...so we will multiply "zero" by the last residual. This might
576: not be exactly what we want to do here -could just return "zero". */
578: *res = 0.0;
579: }
580: return(0);
581: }
583: /*
585: KSPLGMRESGetNewVectors - This routine allocates more work vectors, starting from
586: VEC_VV(it)
588: */
589: static PetscErrorCode KSPLGMRESGetNewVectors(KSP ksp,PetscInt it)590: {
591: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
592: PetscInt nwork = lgmres->nwork_alloc; /* number of work vector chunks allocated */
593: PetscInt nalloc; /* number to allocate */
595: PetscInt k;
598: nalloc = lgmres->delta_allocate; /* number of vectors to allocate
599: in a single chunk */
601: /* Adjust the number to allocate to make sure that we don't exceed the
602: number of available slots (lgmres->vecs_allocated)*/
603: if (it + VEC_OFFSET + nalloc >= lgmres->vecs_allocated) {
604: nalloc = lgmres->vecs_allocated - it - VEC_OFFSET;
605: }
606: if (!nalloc) return(0);
608: lgmres->vv_allocated += nalloc; /* vv_allocated is the number of vectors allocated */
610: /* work vectors */
611: KSPCreateVecs(ksp,nalloc,&lgmres->user_work[nwork],0,NULL);
612: PetscLogObjectParents(ksp,nalloc,lgmres->user_work[nwork]);
613: /* specify size of chunk allocated */
614: lgmres->mwork_alloc[nwork] = nalloc;
616: for (k=0; k < nalloc; k++) {
617: lgmres->vecs[it+VEC_OFFSET+k] = lgmres->user_work[nwork][k];
618: }
621: /* LGMRES_MOD - for now we are preallocating the augmentation vectors */
624: /* increment the number of work vector chunks */
625: lgmres->nwork_alloc++;
626: return(0);
627: }
629: /*
631: KSPBuildSolution_LGMRES
633: Input Parameter:
634: . ksp - the Krylov space object
635: . ptr-
637: Output Parameter:
638: . result - the solution
640: Note: this calls KSPLGMRESBuildSoln - the same function that KSPLGMRESCycle
641: calls directly.
643: */
644: PetscErrorCode KSPBuildSolution_LGMRES(KSP ksp,Vec ptr,Vec *result)645: {
646: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
650: if (!ptr) {
651: if (!lgmres->sol_temp) {
652: VecDuplicate(ksp->vec_sol,&lgmres->sol_temp);
653: PetscLogObjectParent((PetscObject)ksp,(PetscObject)lgmres->sol_temp);
654: }
655: ptr = lgmres->sol_temp;
656: }
657: if (!lgmres->nrs) {
658: /* allocate the work area */
659: PetscMalloc1(lgmres->max_k,&lgmres->nrs);
660: PetscLogObjectMemory((PetscObject)ksp,lgmres->max_k*sizeof(PetscScalar));
661: }
663: KSPLGMRESBuildSoln(lgmres->nrs,ksp->vec_sol,ptr,ksp,lgmres->it);
664: if (result) *result = ptr;
665: return(0);
666: }
668: PetscErrorCode KSPView_LGMRES(KSP ksp,PetscViewer viewer)669: {
670: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
672: PetscBool iascii;
675: KSPView_GMRES(ksp,viewer);
676: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
677: if (iascii) {
678: /*LGMRES_MOD */
679: PetscViewerASCIIPrintf(viewer," aug. dimension=%D\n",lgmres->aug_dim);
680: if (lgmres->approx_constant) {
681: PetscViewerASCIIPrintf(viewer," approx. space size was kept constant.\n");
682: }
683: PetscViewerASCIIPrintf(viewer," number of matvecs=%D\n",lgmres->matvecs);
684: }
685: return(0);
686: }
688: PetscErrorCode KSPSetFromOptions_LGMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)689: {
691: PetscInt aug;
692: KSP_LGMRES *lgmres = (KSP_LGMRES*) ksp->data;
693: PetscBool flg = PETSC_FALSE;
696: KSPSetFromOptions_GMRES(PetscOptionsObject,ksp);
697: PetscOptionsHead(PetscOptionsObject,"KSP LGMRES Options");
698: PetscOptionsBool("-ksp_lgmres_constant","Use constant approx. space size","KSPGMRESSetConstant",lgmres->approx_constant,&lgmres->approx_constant,NULL);
699: PetscOptionsInt("-ksp_lgmres_augment","Number of error approximations to augment the Krylov space with","KSPLGMRESSetAugDim",lgmres->aug_dim,&aug,&flg);
700: if (flg) { KSPLGMRESSetAugDim(ksp,aug); }
701: PetscOptionsTail();
702: return(0);
703: }
705: /*functions for extra lgmres options here*/
706: static PetscErrorCode KSPLGMRESSetConstant_LGMRES(KSP ksp)707: {
708: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
711: lgmres->approx_constant = PETSC_TRUE;
712: return(0);
713: }
715: static PetscErrorCode KSPLGMRESSetAugDim_LGMRES(KSP ksp,PetscInt aug_dim)716: {
717: KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;
720: if (aug_dim < 0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Augmentation dimension must be positive");
721: if (aug_dim > (lgmres->max_k -1)) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Augmentation dimension must be <= (restart size-1)");
722: lgmres->aug_dim = aug_dim;
723: return(0);
724: }
726: /* end new lgmres functions */
728: /*MC
729: KSPLGMRES - Augments the standard GMRES approximation space with approximations to
730: the error from previous restart cycles.
732: Options Database Keys:
733: + -ksp_gmres_restart <restart> - total approximation space size (Krylov directions + error approximations)
734: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
735: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
736: vectors are allocated as needed)
737: . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
738: . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
739: . -ksp_gmres_cgs_refinement_type <never,ifneeded,always> - determine if iterative refinement is used to increase the
740: stability of the classical Gram-Schmidt orthogonalization.
741: . -ksp_gmres_krylov_monitor - plot the Krylov space generated
742: . -ksp_lgmres_augment <k> - number of error approximations to augment the Krylov space with
743: - -ksp_lgmres_constant - use a constant approx. space size (only affects restart cycles < num. error approx.(k), i.e. the first k restarts)
745: To run LGMRES(m, k) as described in the above paper, use:
746: -ksp_gmres_restart <m+k>
747: -ksp_lgmres_augment <k>
749: Level: beginner
751: Notes:
752: Supports both left and right preconditioning, but not symmetric.
754: References:
755: . 1. - A. H. Baker, E.R. Jessup, and T.A. Manteuffel. A technique for accelerating the convergence of restarted GMRES. SIAM Journal on Matrix Analysis and Applications, 26 (2005).
757: Developer Notes:
758: This object is subclassed off of KSPGMRES760: Contributed by: Allison Baker
762: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPGMRES,
763: KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
764: KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
765: KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPLGMRESSetAugDim(),
766: KSPGMRESSetConstant()
768: M*/
770: PETSC_EXTERN PetscErrorCode KSPCreate_LGMRES(KSP ksp)771: {
772: KSP_LGMRES *lgmres;
776: PetscNewLog(ksp,&lgmres);
778: ksp->data = (void*)lgmres;
779: ksp->ops->buildsolution = KSPBuildSolution_LGMRES;
781: ksp->ops->setup = KSPSetUp_LGMRES;
782: ksp->ops->solve = KSPSolve_LGMRES;
783: ksp->ops->destroy = KSPDestroy_LGMRES;
784: ksp->ops->view = KSPView_LGMRES;
785: ksp->ops->setfromoptions = KSPSetFromOptions_LGMRES;
786: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
787: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
789: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
790: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
792: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
793: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
794: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
795: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
796: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
797: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
798: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
799: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);
801: /*LGMRES_MOD add extra functions here - like the one to set num of aug vectors */
802: PetscObjectComposeFunction((PetscObject)ksp,"KSPLGMRESSetConstant_C",KSPLGMRESSetConstant_LGMRES);
803: PetscObjectComposeFunction((PetscObject)ksp,"KSPLGMRESSetAugDim_C",KSPLGMRESSetAugDim_LGMRES);
806: /*defaults */
807: lgmres->haptol = 1.0e-30;
808: lgmres->q_preallocate = 0;
809: lgmres->delta_allocate = LGMRES_DELTA_DIRECTIONS;
810: lgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
811: lgmres->nrs = 0;
812: lgmres->sol_temp = 0;
813: lgmres->max_k = LGMRES_DEFAULT_MAXK;
814: lgmres->Rsvd = 0;
815: lgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
816: lgmres->orthogwork = 0;
818: /*LGMRES_MOD - new defaults */
819: lgmres->aug_dim = LGMRES_DEFAULT_AUGDIM;
820: lgmres->aug_ct = 0; /* start with no aug vectors */
821: lgmres->approx_constant = PETSC_FALSE;
822: lgmres->matvecs = 0;
823: return(0);
824: }