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