<div dir="ltr">Thanks Jed and Barry,<div><br></div><div>So, just to confirm, </div><div><br></div><div>-- From the KSP_GMRES structure, if I call <font face="courier new, monospace">*HH(a,b)</font>, that will return the row a, column b entry of the Hessenberg matrix (while the back end array <font face="courier new, monospace">*hh_origin</font> array is ordering using the Fortran convention)<br></div><div><br></div><div>-- Matrices are passed into and returned from <font face="courier new, monospace">PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_()</font> using Fortran indexing, and need to be transposed to get back to C ordering</div><div><br></div><div>Are both of these statements correct?</div><div><br></div><div>Cheers, Dave.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, May 20, 2019 at 4:34 PM Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
The little work arrays in GMRES tend to be stored in Fortran ordering; there is no C style p[][] indexing into such arrays. Thus the arrays can safely be sent to LAPACK. The only trick is knowing the two dimensions and as Jed say the "leading dimension parameter. He gave you a place to look<br>
<br>
> On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> <br>
> Dave Lee via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> writes:<br>
> <br>
>> Hi Petsc,<br>
>> <br>
>> I'm attempting to implement a "hookstep" for the SNES trust region solver.<br>
>> Essentially what I'm trying to do is replace the solution of the least<br>
>> squares problem at the end of each GMRES solve with a modified solution<br>
>> with a norm that is constrained to be within the size of the trust region.<br>
>> <br>
>> In order to do this I need to perform an SVD on the Hessenberg matrix,<br>
>> which copying the function KSPComputeExtremeSingularValues(), I'm trying to<br>
>> do by accessing the LAPACK function dgesvd() via the PetscStackCallBLAS()<br>
>> machinery. One thing I'm confused about however is the ordering of the 2D<br>
>> arrays into and out of this function, given that that C and FORTRAN arrays<br>
>> use reverse indexing, ie: C[j+1][i+1] = F[i,j].<br>
>> <br>
>> Given that the Hessenberg matrix has k+1 rows and k columns, should I be<br>
>> still be initializing this as H[row][col] and passing this into<br>
>> PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))<br>
>> or should I be transposing this before passing it in?<br>
> <br>
> LAPACK terminology is with respect to Fortran ordering. There is a<br>
> "leading dimension" parameter so that you can operate on non-contiguous<br>
> blocks. See KSPComputeExtremeSingularValues_GMRES for an example.<br>
> <br>
>> Also for the left and right singular vector matrices that are returned by<br>
>> this function, should I be transposing these before I interpret them as C<br>
>> arrays?<br>
>> <br>
>> I've attached my modified version of gmres.c in case this is helpful. If<br>
>> you grep for DRL (my initials) then you'll see my changes to the code.<br>
>> <br>
>> Cheers, Dave.<br>
>> <br>
>> /*<br>
>> This file implements GMRES (a Generalized Minimal Residual) method.<br>
>> Reference: Saad and Schultz, 1986.<br>
>> <br>
>> <br>
>> Some comments on left vs. right preconditioning, and restarts.<br>
>> Left and right preconditioning.<br>
>> If right preconditioning is chosen, then the problem being solved<br>
>> by gmres is actually<br>
>> My = AB^-1 y = f<br>
>> so the initial residual is<br>
>> r = f - Mx<br>
>> Note that B^-1 y = x or y = B x, and if x is non-zero, the initial<br>
>> residual is<br>
>> r = f - A x<br>
>> The final solution is then<br>
>> x = B^-1 y<br>
>> <br>
>> If left preconditioning is chosen, then the problem being solved is<br>
>> My = B^-1 A x = B^-1 f,<br>
>> and the initial residual is<br>
>> r = B^-1(f - Ax)<br>
>> <br>
>> Restarts: Restarts are basically solves with x0 not equal to zero.<br>
>> Note that we can eliminate an extra application of B^-1 between<br>
>> restarts as long as we don't require that the solution at the end<br>
>> of an unsuccessful gmres iteration always be the solution x.<br>
>> */<br>
>> <br>
>> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I "petscksp.h" I*/<br>
>> #include <petscblaslapack.h> // DRL<br>
>> #define GMRES_DELTA_DIRECTIONS 10<br>
>> #define GMRES_DEFAULT_MAXK 30<br>
>> static PetscErrorCode KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);<br>
>> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);<br>
>> <br>
>> PetscErrorCode KSPSetUp_GMRES(KSP ksp)<br>
>> {<br>
>> PetscInt hh,hes,rs,cc;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt max_k,k;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> max_k = gmres->max_k; /* restart size */<br>
>> hh = (max_k + 2) * (max_k + 1);<br>
>> hes = (max_k + 1) * (max_k + 1);<br>
>> rs = (max_k + 2);<br>
>> cc = (max_k + 1);<br>
>> <br>
>> ierr = PetscCalloc5(hh,&gmres->hh_origin,hes,&gmres->hes_origin,rs,&gmres->rs_origin,cc,&gmres->cc_origin,cc,&gmres->ss_origin);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);<br>
>> <br>
>> if (ksp->calc_sings) {<br>
>> /* Allocate workspace to hold Hessenberg matrix needed by lapack */<br>
>> ierr = PetscMalloc1((max_k + 3)*(max_k + 9),&gmres->Rsvd);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 9)*sizeof(PetscScalar));CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(6*(max_k+2),&gmres->Dsvd);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);<br>
>> }<br>
>> <br>
>> /* Allocate array to hold pointers to user vectors. Note that we need<br>
>> 4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */<br>
>> gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;<br>
>> <br>
>> ierr = PetscMalloc1(gmres->vecs_allocated,&gmres->vecs);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->user_work);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->mwork_alloc);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt)) + gmres->vecs_allocated*sizeof(Vec));CHKERRQ(ierr);<br>
>> <br>
>> if (gmres->q_preallocate) {<br>
>> gmres->vv_allocated = VEC_OFFSET + 2 + max_k;<br>
>> <br>
>> ierr = KSPCreateVecs(ksp,gmres->vv_allocated,&gmres->user_work[0],0,NULL);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);CHKERRQ(ierr);<br>
>> <br>
>> gmres->mwork_alloc[0] = gmres->vv_allocated;<br>
>> gmres->nwork_alloc = 1;<br>
>> for (k=0; k<gmres->vv_allocated; k++) {<br>
>> gmres->vecs[k] = gmres->user_work[0][k];<br>
>> }<br>
>> } else {<br>
>> gmres->vv_allocated = 5;<br>
>> <br>
>> ierr = KSPCreateVecs(ksp,5,&gmres->user_work[0],0,NULL);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectParents(ksp,5,gmres->user_work[0]);CHKERRQ(ierr);<br>
>> <br>
>> gmres->mwork_alloc[0] = 5;<br>
>> gmres->nwork_alloc = 1;<br>
>> for (k=0; k<gmres->vv_allocated; k++) {<br>
>> gmres->vecs[k] = gmres->user_work[0][k];<br>
>> }<br>
>> }<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*<br>
>> Run gmres, possibly with restart. Return residual history if requested.<br>
>> input parameters:<br>
>> <br>
>> . gmres - structure containing parameters and work areas<br>
>> <br>
>> output parameters:<br>
>> . nres - residuals (from preconditioned system) at each step.<br>
>> If restarting, consider passing nres+it. If null,<br>
>> ignored<br>
>> . itcount - number of iterations used. nres[0] to nres[itcount]<br>
>> are defined. If null, ignored.<br>
>> <br>
>> Notes:<br>
>> On entry, the value in vector VEC_VV(0) should be the initial residual<br>
>> (this allows shortcuts where the initial preconditioned residual is 0).<br>
>> */<br>
>> PetscErrorCode KSPGMRESCycle(PetscInt *itcount,KSP ksp)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);<br>
>> PetscReal res_norm,res,hapbnd,tt;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt it = 0, max_k = gmres->max_k;<br>
>> PetscBool hapend = PETSC_FALSE;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> if (itcount) *itcount = 0;<br>
>> ierr = VecNormalize(VEC_VV(0),&res_norm);CHKERRQ(ierr);<br>
>> KSPCheckNorm(ksp,res_norm);<br>
>> res = res_norm;<br>
>> *GRS(0) = res_norm;<br>
>> <br>
>> /* check for the convergence */<br>
>> ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);<br>
>> ksp->rnorm = res;<br>
>> ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);<br>
>> gmres->it = (it - 1);<br>
>> ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);<br>
>> ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);<br>
>> if (!res) {<br>
>> ksp->reason = KSP_CONVERGED_ATOL;<br>
>> ierr = PetscInfo(ksp,"Converged due to zero residual norm on entry\n");CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);<br>
>> while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {<br>
>> if (it) {<br>
>> ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);<br>
>> ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);<br>
>> }<br>
>> gmres->it = (it - 1);<br>
>> if (gmres->vv_allocated <= it + VEC_OFFSET + 1) {<br>
>> ierr = KSPGMRESGetNewVectors(ksp,it+1);CHKERRQ(ierr);<br>
>> }<br>
>> ierr = KSP_PCApplyBAorAB(ksp,VEC_VV(it),VEC_VV(1+it),VEC_TEMP_MATOP);CHKERRQ(ierr);<br>
>> <br>
>> /* update hessenberg matrix and do Gram-Schmidt */<br>
>> ierr = (*gmres->orthog)(ksp,it);CHKERRQ(ierr);<br>
>> if (ksp->reason) break;<br>
>> <br>
>> /* vv(i+1) . vv(i+1) */<br>
>> ierr = VecNormalize(VEC_VV(it+1),&tt);CHKERRQ(ierr);<br>
>> <br>
>> /* save the magnitude */<br>
>> *HH(it+1,it) = tt;<br>
>> *HES(it+1,it) = tt;<br>
>> <br>
>> /* check for the happy breakdown */<br>
>> hapbnd = PetscAbsScalar(tt / *GRS(it));<br>
>> if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;<br>
>> if (tt < hapbnd) {<br>
>> ierr = PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n",(double)hapbnd,(double)tt);CHKERRQ(ierr);<br>
>> hapend = PETSC_TRUE;<br>
>> }<br>
>> ierr = KSPGMRESUpdateHessenberg(ksp,it,hapend,&res);CHKERRQ(ierr);<br>
>> <br>
>> it++;<br>
>> gmres->it = (it-1); /* For converged */<br>
>> ksp->its++;<br>
>> ksp->rnorm = res;<br>
>> if (ksp->reason) break;<br>
>> <br>
>> ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);<br>
>> <br>
>> /* Catch error in happy breakdown and signal convergence and break from loop */<br>
>> if (hapend) {<br>
>> if (!ksp->reason) {<br>
>> 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);<br>
>> else {<br>
>> ksp->reason = KSP_DIVERGED_BREAKDOWN;<br>
>> break;<br>
>> }<br>
>> }<br>
>> }<br>
>> }<br>
>> <br>
>> /* Monitor if we know that we will not return for a restart */<br>
>> if (it && (ksp->reason || ksp->its >= ksp->max_it)) {<br>
>> ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);<br>
>> ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);<br>
>> }<br>
>> <br>
>> if (itcount) *itcount = it;<br>
>> <br>
>> <br>
>> /*<br>
>> Down here we have to solve for the "best" coefficients of the Krylov<br>
>> columns, add the solution values together, and possibly unwind the<br>
>> preconditioning from the solution<br>
>> */<br>
>> /* Form the solution (or the solution so far) */<br>
>> ierr = KSPGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-1);CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPSolve_GMRES(KSP ksp)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> PetscInt its,itcount,i;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscBool guess_zero = ksp->guess_zero;<br>
>> PetscInt N = gmres->max_k + 1;<br>
>> PetscBLASInt bN;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> if (ksp->calc_sings && !gmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");<br>
>> <br>
>> ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);<br>
>> ksp->its = 0;<br>
>> ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);<br>
>> <br>
>> itcount = 0;<br>
>> gmres->fullcycle = 0;<br>
>> ksp->reason = KSP_CONVERGED_ITERATING;<br>
>> while (!ksp->reason) {<br>
>> ierr = KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);CHKERRQ(ierr);<br>
>> ierr = KSPGMRESCycle(&its,ksp);CHKERRQ(ierr);<br>
>> /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace<br>
>> if the cycle is complete for the computation of the Ritz pairs */<br>
>> if (its == gmres->max_k) {<br>
>> gmres->fullcycle++;<br>
>> if (ksp->calc_ritz) {<br>
>> if (!gmres->hes_ritz) {<br>
>> ierr = PetscMalloc1(N*N,&gmres->hes_ritz);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,N*N*sizeof(PetscScalar));CHKERRQ(ierr);<br>
>> ierr = VecDuplicateVecs(VEC_VV(0),N,&gmres->vecb);CHKERRQ(ierr);<br>
>> }<br>
>> ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr);<br>
>> ierr = PetscMemcpy(gmres->hes_ritz,gmres->hes_origin,bN*bN*sizeof(PetscReal));CHKERRQ(ierr);<br>
>> for (i=0; i<gmres->max_k+1; i++) {<br>
>> ierr = VecCopy(VEC_VV(i),gmres->vecb[i]);CHKERRQ(ierr);<br>
>> }<br>
>> }<br>
>> }<br>
>> itcount += its;<br>
>> if (itcount >= ksp->max_it) {<br>
>> if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;<br>
>> break;<br>
>> }<br>
>> ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */<br>
>> }<br>
>> ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPReset_GMRES(KSP ksp)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt i;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> /* Free the Hessenberg matrices */<br>
>> ierr = PetscFree6(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin,gmres->hes_ritz);CHKERRQ(ierr);<br>
>> <br>
>> /* free work vectors */<br>
>> ierr = PetscFree(gmres->vecs);CHKERRQ(ierr);<br>
>> for (i=0; i<gmres->nwork_alloc; i++) {<br>
>> ierr = VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);CHKERRQ(ierr);<br>
>> }<br>
>> gmres->nwork_alloc = 0;<br>
>> if (gmres->vecb) {<br>
>> ierr = VecDestroyVecs(gmres->max_k+1,&gmres->vecb);CHKERRQ(ierr);<br>
>> }<br>
>> <br>
>> ierr = PetscFree(gmres->user_work);CHKERRQ(ierr);<br>
>> ierr = PetscFree(gmres->mwork_alloc);CHKERRQ(ierr);<br>
>> ierr = PetscFree(gmres->nrs);CHKERRQ(ierr);<br>
>> ierr = VecDestroy(&gmres->sol_temp);CHKERRQ(ierr);<br>
>> ierr = PetscFree(gmres->Rsvd);CHKERRQ(ierr);<br>
>> ierr = PetscFree(gmres->Dsvd);CHKERRQ(ierr);<br>
>> ierr = PetscFree(gmres->orthogwork);CHKERRQ(ierr);<br>
>> <br>
>> gmres->sol_temp = 0;<br>
>> gmres->vv_allocated = 0;<br>
>> gmres->vecs_allocated = 0;<br>
>> gmres->sol_temp = 0;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPDestroy_GMRES(KSP ksp)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = KSPReset_GMRES(ksp);CHKERRQ(ierr);<br>
>> ierr = PetscFree(ksp->data);CHKERRQ(ierr);<br>
>> /* clear composed functions */<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",NULL);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",NULL);CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> /*<br>
>> KSPGMRESBuildSoln - create the solution from the starting vector and the<br>
>> current iterates.<br>
>> <br>
>> Input parameters:<br>
>> nrs - work area of size it + 1.<br>
>> vs - index of initial guess<br>
>> vdest - index of result. Note that vs may == vdest (replace<br>
>> guess with the solution).<br>
>> <br>
>> This is an internal routine that knows about the GMRES internals.<br>
>> */<br>
>> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *nrs,Vec vs,Vec vdest,KSP ksp,PetscInt it)<br>
>> {<br>
>> PetscScalar tt;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt ii,k,j;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);<br>
>> <br>
>> PetscFunctionBegin;<br>
>> /* Solve for solution vector that minimizes the residual */<br>
>> <br>
>> /* If it is < 0, no gmres steps have been performed */<br>
>> if (it < 0) {<br>
>> ierr = VecCopy(vs,vdest);CHKERRQ(ierr); /* VecCopy() is smart, exists immediately if vguess == vdest */<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> if (*HH(it,it) != 0.0) {<br>
>> nrs[it] = *GRS(it) / *HH(it,it);<br>
>> } else {<br>
>> ksp->reason = KSP_DIVERGED_BREAKDOWN;<br>
>> <br>
>> ierr = PetscInfo2(ksp,"Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %D GRS(it) = %g\n",it,(double)PetscAbsScalar(*GRS(it)));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> for (ii=1; ii<=it; ii++) {<br>
>> k = it - ii;<br>
>> tt = *GRS(k);<br>
>> for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];<br>
>> if (*HH(k,k) == 0.0) {<br>
>> ksp->reason = KSP_DIVERGED_BREAKDOWN;<br>
>> <br>
>> ierr = PetscInfo1(ksp,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> nrs[k] = tt / *HH(k,k);<br>
>> }<br>
>> <br>
>> /* Perform the hookstep correction - DRL */<br>
>> if(gmres->delta > 0.0 && gmres->it > 0) { // Apply the hookstep to correct the GMRES solution (if required)<br>
>> printf("\t\tapplying hookstep: initial delta: %lf", gmres->delta);<br>
>> PetscInt N = gmres->max_k+2, ii, jj, j0;<br>
>> PetscBLASInt nRows, nCols, lwork, lierr;<br>
>> PetscScalar *R, *work;<br>
>> PetscReal* S;<br>
>> PetscScalar *U, *VT, *p, *q, *y;<br>
>> PetscScalar bnorm, mu, qMag, qMag2, delta2;<br>
>> <br>
>> ierr = PetscMalloc1((gmres->max_k + 3)*(gmres->max_k + 9),&R);CHKERRQ(ierr);<br>
>> work = R + N*N;<br>
>> ierr = PetscMalloc1(6*(gmres->max_k+2),&S);CHKERRQ(ierr);<br>
>> <br>
>> ierr = PetscBLASIntCast(gmres->it+1,&nRows);CHKERRQ(ierr);<br>
>> ierr = PetscBLASIntCast(gmres->it+0,&nCols);CHKERRQ(ierr);<br>
>> ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr);<br>
>> //ierr = PetscMemcpy(R,gmres->hes_origin,(gmres->max_k+2)*(gmres->max_k+1)*sizeof(PetscScalar));CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(nRows*nCols,&R);CHKERRQ(ierr);<br>
>> for (ii = 0; ii < nRows; ii++) {<br>
>> for (jj = 0; jj < nCols; jj++) {<br>
>> R[ii*nCols+jj] = *HH(ii,jj);<br>
>> // Ensure Hessenberg structure<br>
>> //if (ii > jj+1) R[ii*nCols+jj] = 0.0;<br>
>> }<br>
>> }<br>
>> <br>
>> ierr = PetscMalloc1(nRows*nRows,&U);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(nCols*nCols,&VT);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(nRows,&p);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(nCols,&q);CHKERRQ(ierr);<br>
>> ierr = PetscMalloc1(nRows,&y);CHKERRQ(ierr);<br>
>> <br>
>> printf("\n\n");for(ii=0;ii<nRows;ii++){for(jj=0;jj<nCols;jj++){printf("\t%g",R[ii*nCols+jj]);}printf("\n");}printf("\n");<br>
>> <br>
>> // Perform an SVD on the Hessenberg matrix<br>
>> ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);<br>
>> PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nRows,&nCols,R,&nRows,S,U,&nRows,VT,&nCols,work,&lwork,&lierr));<br>
>> if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr);<br>
>> ierr = PetscFPTrapPop();CHKERRQ(ierr);<br>
>> <br>
>> // Compute p = ||b|| U^T e_1<br>
>> ierr = VecNorm(ksp->vec_rhs,NORM_2,&bnorm);CHKERRQ(ierr);<br>
>> for (ii=0; ii<nRows; ii++) {<br>
>> p[ii] = bnorm*U[ii*nRows];<br>
>> }<br>
>> <br>
>> // Solve the root finding problem for \mu such that ||q|| < \delta (where \delta is the radius of the trust region)<br>
>> // This step is largely copied from Ashley Willis' openpipeflow: <a href="http://doi.org/10.1016/j.softx.2017.05.003" rel="noreferrer" target="_blank">doi.org/10.1016/j.softx.2017.05.003</a><br>
>> mu = S[nCols-1]*S[nCols-1]*1.0e-6;<br>
>> if (mu < 1.0e-99) mu = 1.0e-99;<br>
>> qMag = 1.0e+99;<br>
>> <br>
>> while (qMag > gmres->delta) {<br>
>> mu *= 1.1;<br>
>> qMag2 = 0.0;<br>
>> for (ii=0; ii<nCols; ii++) {<br>
>> q[ii] = p[ii]*S[ii]/(mu + S[ii]*S[ii]);<br>
>> qMag2 += q[ii]*q[ii];<br>
>> }<br>
>> qMag = PetscSqrtScalar(qMag2);<br>
>> }<br>
>> <br>
>> // Expand y in terms of the right singular vectors as y = V q<br>
>> for (ii=0; ii<nCols; ii++) {<br>
>> y[ii] = 0.0;<br>
>> for (jj=0; jj<nCols; jj++) {<br>
>> y[ii] += VT[jj*nCols+ii]*q[jj]; // transpose of the transpose<br>
>> }<br>
>> }<br>
>> <br>
>> // Recompute the size of the trust region, \delta<br>
>> delta2 = 0.0;<br>
>> for (ii=0; ii<nRows; ii++) {<br>
>> j0 = (ii < 2) ? 0 : ii - 1;<br>
>> p[ii] = 0.0;<br>
>> for (jj=j0; jj<nCols; jj++) {<br>
>> p[ii] -= R[ii*nCols+jj]*y[jj];<br>
>> }<br>
>> if (ii == 0) {<br>
>> p[ii] += bnorm;<br>
>> }<br>
>> delta2 += p[ii]*p[ii];<br>
>> }<br>
>> gmres->delta = PetscSqrtScalar(delta2);<br>
>> printf("\t\t...final delta: %lf.\n", gmres->delta);<br>
>> <br>
>> // Pass the orthnomalized Krylov vector weights back out<br>
>> for (ii=0; ii<nCols; ii++) {<br>
>> nrs[ii] = y[ii];<br>
>> }<br>
>> <br>
>> ierr = PetscFree(R);CHKERRQ(ierr);<br>
>> ierr = PetscFree(S);CHKERRQ(ierr);<br>
>> ierr = PetscFree(U);CHKERRQ(ierr);<br>
>> ierr = PetscFree(VT);CHKERRQ(ierr);<br>
>> ierr = PetscFree(p);CHKERRQ(ierr);<br>
>> ierr = PetscFree(q);CHKERRQ(ierr);<br>
>> ierr = PetscFree(y);CHKERRQ(ierr);<br>
>> }<br>
>> /*** DRL ***/<br>
>> <br>
>> /* Accumulate the correction to the solution of the preconditioned problem in TEMP */<br>
>> ierr = VecSet(VEC_TEMP,0.0);CHKERRQ(ierr);<br>
>> if (gmres->delta > 0.0) {<br>
>> ierr = VecMAXPY(VEC_TEMP,it,nrs,&VEC_VV(0));CHKERRQ(ierr); // DRL<br>
>> } else {<br>
>> ierr = VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));CHKERRQ(ierr); <br>
>> }<br>
>> <br>
>> ierr = KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);CHKERRQ(ierr);<br>
>> /* add solution to previous solution */<br>
>> if (vdest != vs) {<br>
>> ierr = VecCopy(vs,vdest);CHKERRQ(ierr);<br>
>> }<br>
>> ierr = VecAXPY(vdest,1.0,VEC_TEMP);CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> /*<br>
>> Do the scalar work for the orthogonalization. Return new residual norm.<br>
>> */<br>
>> static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)<br>
>> {<br>
>> PetscScalar *hh,*cc,*ss,tt;<br>
>> PetscInt j;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);<br>
>> <br>
>> PetscFunctionBegin;<br>
>> hh = HH(0,it);<br>
>> cc = CC(0);<br>
>> ss = SS(0);<br>
>> <br>
>> /* Apply all the previously computed plane rotations to the new column<br>
>> of the Hessenberg matrix */<br>
>> for (j=1; j<=it; j++) {<br>
>> tt = *hh;<br>
>> *hh = PetscConj(*cc) * tt + *ss * *(hh+1);<br>
>> hh++;<br>
>> *hh = *cc++ * *hh - (*ss++ * tt);<br>
>> }<br>
>> <br>
>> /*<br>
>> compute the new plane rotation, and apply it to:<br>
>> 1) the right-hand-side of the Hessenberg system<br>
>> 2) the new column of the Hessenberg matrix<br>
>> thus obtaining the updated value of the residual<br>
>> */<br>
>> if (!hapend) {<br>
>> tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));<br>
>> if (tt == 0.0) {<br>
>> ksp->reason = KSP_DIVERGED_NULL;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> *cc = *hh / tt;<br>
>> *ss = *(hh+1) / tt;<br>
>> *GRS(it+1) = -(*ss * *GRS(it));<br>
>> *GRS(it) = PetscConj(*cc) * *GRS(it);<br>
>> *hh = PetscConj(*cc) * *hh + *ss * *(hh+1);<br>
>> *res = PetscAbsScalar(*GRS(it+1));<br>
>> } else {<br>
>> /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply<br>
>> another rotation matrix (so RH doesn't change). The new residual is<br>
>> always the new sine term times the residual from last time (GRS(it)),<br>
>> but now the new sine rotation would be zero...so the residual should<br>
>> be zero...so we will multiply "zero" by the last residual. This might<br>
>> not be exactly what we want to do here -could just return "zero". */<br>
>> <br>
>> *res = 0.0;<br>
>> }<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> /*<br>
>> This routine allocates more work vectors, starting from VEC_VV(it).<br>
>> */<br>
>> PetscErrorCode KSPGMRESGetNewVectors(KSP ksp,PetscInt it)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt nwork = gmres->nwork_alloc,k,nalloc;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> nalloc = PetscMin(ksp->max_it,gmres->delta_allocate);<br>
>> /* Adjust the number to allocate to make sure that we don't exceed the<br>
>> number of available slots */<br>
>> if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) {<br>
>> nalloc = gmres->vecs_allocated - it - VEC_OFFSET;<br>
>> }<br>
>> if (!nalloc) PetscFunctionReturn(0);<br>
>> <br>
>> gmres->vv_allocated += nalloc;<br>
>> <br>
>> ierr = KSPCreateVecs(ksp,nalloc,&gmres->user_work[nwork],0,NULL);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectParents(ksp,nalloc,gmres->user_work[nwork]);CHKERRQ(ierr);<br>
>> <br>
>> gmres->mwork_alloc[nwork] = nalloc;<br>
>> for (k=0; k<nalloc; k++) {<br>
>> gmres->vecs[it+VEC_OFFSET+k] = gmres->user_work[nwork][k];<br>
>> }<br>
>> gmres->nwork_alloc++;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec *result)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> if (!ptr) {<br>
>> if (!gmres->sol_temp) {<br>
>> ierr = VecDuplicate(ksp->vec_sol,&gmres->sol_temp);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)gmres->sol_temp);CHKERRQ(ierr);<br>
>> }<br>
>> ptr = gmres->sol_temp;<br>
>> }<br>
>> if (!gmres->nrs) {<br>
>> /* allocate the work area */<br>
>> ierr = PetscMalloc1(gmres->max_k,&gmres->nrs);CHKERRQ(ierr);<br>
>> ierr = PetscLogObjectMemory((PetscObject)ksp,gmres->max_k*sizeof(PetscScalar));CHKERRQ(ierr);<br>
>> }<br>
>> <br>
>> ierr = KSPGMRESBuildSoln(gmres->nrs,ksp->vec_sol,ptr,ksp,gmres->it);CHKERRQ(ierr);<br>
>> if (result) *result = ptr;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPView_GMRES(KSP ksp,PetscViewer viewer)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> const char *cstr;<br>
>> PetscErrorCode ierr;<br>
>> PetscBool iascii,isstring;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);<br>
>> ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);<br>
>> if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {<br>
>> switch (gmres->cgstype) {<br>
>> case (KSP_GMRES_CGS_REFINE_NEVER):<br>
>> cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";<br>
>> break;<br>
>> case (KSP_GMRES_CGS_REFINE_ALWAYS):<br>
>> cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";<br>
>> break;<br>
>> case (KSP_GMRES_CGS_REFINE_IFNEEDED):<br>
>> cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";<br>
>> break;<br>
>> default:<br>
>> SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Unknown orthogonalization");<br>
>> }<br>
>> } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {<br>
>> cstr = "Modified Gram-Schmidt Orthogonalization";<br>
>> } else {<br>
>> cstr = "unknown orthogonalization";<br>
>> }<br>
>> if (iascii) {<br>
>> ierr = PetscViewerASCIIPrintf(viewer," restart=%D, using %s\n",gmres->max_k,cstr);CHKERRQ(ierr);<br>
>> ierr = PetscViewerASCIIPrintf(viewer," happy breakdown tolerance %g\n",(double)gmres->haptol);CHKERRQ(ierr);<br>
>> } else if (isstring) {<br>
>> ierr = PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,gmres->max_k);CHKERRQ(ierr);<br>
>> }<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*@C<br>
>> KSPGMRESMonitorKrylov - Calls VecView() for each new direction in the GMRES accumulated Krylov space.<br>
>> <br>
>> Collective on KSP<br>
>> <br>
>> Input Parameters:<br>
>> + ksp - the KSP context<br>
>> . its - iteration number<br>
>> . fgnorm - 2-norm of residual (or gradient)<br>
>> - dummy - an collection of viewers created with KSPViewerCreate()<br>
>> <br>
>> Options Database Keys:<br>
>> . -ksp_gmres_kyrlov_monitor<br>
>> <br>
>> Notes: A new PETSCVIEWERDRAW is created for each Krylov vector so they can all be simultaneously viewed<br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, nonlinear, vector, monitor, view, Krylov space<br>
>> <br>
>> .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView(), KSPViewersCreate(), KSPViewersDestroy()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)<br>
>> {<br>
>> PetscViewers viewers = (PetscViewers)dummy;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscErrorCode ierr;<br>
>> Vec x;<br>
>> PetscViewer viewer;<br>
>> PetscBool flg;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = PetscViewersGetViewer(viewers,gmres->it+1,&viewer);CHKERRQ(ierr);<br>
>> ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&flg);CHKERRQ(ierr);<br>
>> if (!flg) {<br>
>> ierr = PetscViewerSetType(viewer,PETSCVIEWERDRAW);CHKERRQ(ierr);<br>
>> ierr = PetscViewerDrawSetInfo(viewer,NULL,"Krylov GMRES Monitor",PETSC_DECIDE,PETSC_DECIDE,300,300);CHKERRQ(ierr);<br>
>> }<br>
>> x = VEC_VV(gmres->it+1);<br>
>> ierr = VecView(x,viewer);CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPSetFromOptions_GMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> PetscInt restart;<br>
>> PetscReal haptol;<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscBool flg;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = PetscOptionsHead(PetscOptionsObject,"KSP GMRES Options");CHKERRQ(ierr);<br>
>> ierr = PetscOptionsInt("-ksp_gmres_restart","Number of Krylov search directions","KSPGMRESSetRestart",gmres->max_k,&restart,&flg);CHKERRQ(ierr);<br>
>> if (flg) { ierr = KSPGMRESSetRestart(ksp,restart);CHKERRQ(ierr); }<br>
>> ierr = PetscOptionsReal("-ksp_gmres_haptol","Tolerance for exact convergence (happy ending)","KSPGMRESSetHapTol",gmres->haptol,&haptol,&flg);CHKERRQ(ierr);<br>
>> if (flg) { ierr = KSPGMRESSetHapTol(ksp,haptol);CHKERRQ(ierr); }<br>
>> flg = PETSC_FALSE;<br>
>> ierr = PetscOptionsBool("-ksp_gmres_preallocate","Preallocate Krylov vectors","KSPGMRESSetPreAllocateVectors",flg,&flg,NULL);CHKERRQ(ierr);<br>
>> if (flg) {ierr = KSPGMRESSetPreAllocateVectors(ksp);CHKERRQ(ierr);}<br>
>> ierr = PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt","Classical (unmodified) Gram-Schmidt (fast)","KSPGMRESSetOrthogonalization",&flg);CHKERRQ(ierr);<br>
>> if (flg) {ierr = KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization);CHKERRQ(ierr);}<br>
>> ierr = PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt","Modified Gram-Schmidt (slow,more stable)","KSPGMRESSetOrthogonalization",&flg);CHKERRQ(ierr);<br>
>> if (flg) {ierr = KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization);CHKERRQ(ierr);}<br>
>> ierr = PetscOptionsEnum("-ksp_gmres_cgs_refinement_type","Type of iterative refinement for classical (unmodified) Gram-Schmidt","KSPGMRESSetCGSRefinementType",<br>
>> KSPGMRESCGSRefinementTypes,(PetscEnum)gmres->cgstype,(PetscEnum*)&gmres->cgstype,&flg);CHKERRQ(ierr);<br>
>> flg = PETSC_FALSE;<br>
>> ierr = PetscOptionsBool("-ksp_gmres_krylov_monitor","Plot the Krylov directions","KSPMonitorSet",flg,&flg,NULL);CHKERRQ(ierr);<br>
>> if (flg) {<br>
>> PetscViewers viewers;<br>
>> ierr = PetscViewersCreate(PetscObjectComm((PetscObject)ksp),&viewers);CHKERRQ(ierr);<br>
>> ierr = KSPMonitorSet(ksp,KSPGMRESMonitorKrylov,viewers,(PetscErrorCode (*)(void**))PetscViewersDestroy);CHKERRQ(ierr);<br>
>> }<br>
>> ierr = PetscOptionsTail();CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> if (tol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");<br>
>> gmres->haptol = tol;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt *max_k)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> *max_k = gmres->max_k;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> if (max_k < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");<br>
>> if (!ksp->setupstage) {<br>
>> gmres->max_k = max_k;<br>
>> } else if (gmres->max_k != max_k) {<br>
>> gmres->max_k = max_k;<br>
>> ksp->setupstage = KSP_SETUP_NEW;<br>
>> /* free the data structures, then create them again */<br>
>> ierr = KSPReset_GMRES(ksp);CHKERRQ(ierr);<br>
>> }<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)<br>
>> {<br>
>> PetscFunctionBegin;<br>
>> ((KSP_GMRES*)ksp->data)->orthog = fcn;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN *fcn)<br>
>> {<br>
>> PetscFunctionBegin;<br>
>> *fcn = ((KSP_GMRES*)ksp->data)->orthog;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)<br>
>> {<br>
>> KSP_GMRES *gmres;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> gmres = (KSP_GMRES*)ksp->data;<br>
>> gmres->q_preallocate = 1;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> gmres->cgstype = type;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType *type)<br>
>> {<br>
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> *type = gmres->cgstype;<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*@<br>
>> KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use<br>
>> in the classical Gram Schmidt orthogonalization.<br>
>> <br>
>> Logically Collective on KSP<br>
>> <br>
>> Input Parameters:<br>
>> + ksp - the Krylov space context<br>
>> - type - the type of refinement<br>
>> <br>
>> Options Database:<br>
>> . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always><br>
>> <br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, GMRES, iterative refinement<br>
>> <br>
>> .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESGetCGSRefinementType(),<br>
>> KSPGMRESGetOrthogonalization()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);<br>
>> PetscValidLogicalCollectiveEnum(ksp,type,2);<br>
>> ierr = PetscTryMethod(ksp,"KSPGMRESSetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType),(ksp,type));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*@<br>
>> KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use<br>
>> in the classical Gram Schmidt orthogonalization.<br>
>> <br>
>> Not Collective<br>
>> <br>
>> Input Parameter:<br>
>> . ksp - the Krylov space context<br>
>> <br>
>> Output Parameter:<br>
>> . type - the type of refinement<br>
>> <br>
>> Options Database:<br>
>> . -ksp_gmres_cgs_refinement_type <never,ifneeded,always><br>
>> <br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, GMRES, iterative refinement<br>
>> <br>
>> .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),<br>
>> KSPGMRESGetOrthogonalization()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType *type)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);<br>
>> ierr = PetscUseMethod(ksp,"KSPGMRESGetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType*),(ksp,type));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> <br>
>> /*@<br>
>> KSPGMRESSetRestart - Sets number of iterations at which GMRES, FGMRES and LGMRES restarts.<br>
>> <br>
>> Logically Collective on KSP<br>
>> <br>
>> Input Parameters:<br>
>> + ksp - the Krylov space context<br>
>> - restart - integer restart value<br>
>> <br>
>> Options Database:<br>
>> . -ksp_gmres_restart <positive integer><br>
>> <br>
>> Note: The default value is 30.<br>
>> <br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, GMRES, restart, iterations<br>
>> <br>
>> .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESGetRestart()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> PetscValidLogicalCollectiveInt(ksp,restart,2);<br>
>> <br>
>> ierr = PetscTryMethod(ksp,"KSPGMRESSetRestart_C",(KSP,PetscInt),(ksp,restart));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*@<br>
>> KSPGMRESGetRestart - Gets number of iterations at which GMRES, FGMRES and LGMRES restarts.<br>
>> <br>
>> Not Collective<br>
>> <br>
>> Input Parameter:<br>
>> . ksp - the Krylov space context<br>
>> <br>
>> Output Parameter:<br>
>> . restart - integer restart value<br>
>> <br>
>> Note: The default value is 30.<br>
>> <br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, GMRES, restart, iterations<br>
>> <br>
>> .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetRestart()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = PetscUseMethod(ksp,"KSPGMRESGetRestart_C",(KSP,PetscInt*),(ksp,restart));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*@<br>
>> KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in GMRES, FGMRES and LGMRES.<br>
>> <br>
>> Logically Collective on KSP<br>
>> <br>
>> Input Parameters:<br>
>> + ksp - the Krylov space context<br>
>> - tol - the tolerance<br>
>> <br>
>> Options Database:<br>
>> . -ksp_gmres_haptol <positive real value><br>
>> <br>
>> Note: Happy breakdown is the rare case in GMRES where an 'exact' solution is obtained after<br>
>> a certain number of iterations. If you attempt more iterations after this point unstable<br>
>> things can happen hence very occasionally you may need to set this value to detect this condition<br>
>> <br>
>> Level: intermediate<br>
>> <br>
>> .keywords: KSP, GMRES, tolerance<br>
>> <br>
>> .seealso: KSPSetTolerances()<br>
>> @*/<br>
>> PetscErrorCode KSPGMRESSetHapTol(KSP ksp,PetscReal tol)<br>
>> {<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> PetscValidLogicalCollectiveReal(ksp,tol,2);<br>
>> ierr = PetscTryMethod((ksp),"KSPGMRESSetHapTol_C",(KSP,PetscReal),((ksp),(tol)));CHKERRQ(ierr);<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>> <br>
>> /*MC<br>
>> KSPGMRES - Implements the Generalized Minimal Residual method.<br>
>> (Saad and Schultz, 1986) with restart<br>
>> <br>
>> <br>
>> Options Database Keys:<br>
>> + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against<br>
>> . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)<br>
>> . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of<br>
>> vectors are allocated as needed)<br>
>> . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)<br>
>> . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)<br>
>> . -ksp_gmres_cgs_refinement_type <never,ifneeded,always> - determine if iterative refinement is used to increase the<br>
>> stability of the classical Gram-Schmidt orthogonalization.<br>
>> - -ksp_gmres_krylov_monitor - plot the Krylov space generated<br>
>> <br>
>> Level: beginner<br>
>> <br>
>> Notes: Left and right preconditioning are supported, but not symmetric preconditioning.<br>
>> <br>
>> References:<br>
>> . 1. - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.<br>
>> SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.<br>
>> <br>
>> .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,<br>
>> KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),<br>
>> KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),<br>
>> KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()<br>
>> <br>
>> M*/<br>
>> <br>
>> PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)<br>
>> {<br>
>> KSP_GMRES *gmres;<br>
>> PetscErrorCode ierr;<br>
>> <br>
>> PetscFunctionBegin;<br>
>> ierr = PetscNewLog(ksp,&gmres);CHKERRQ(ierr);<br>
>> ksp->data = (void*)gmres;<br>
>> <br>
>> ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,4);CHKERRQ(ierr);<br>
>> ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);CHKERRQ(ierr);<br>
>> ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,2);CHKERRQ(ierr);<br>
>> ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);CHKERRQ(ierr);<br>
>> ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);<br>
>> <br>
>> ksp->ops->buildsolution = KSPBuildSolution_GMRES;<br>
>> ksp->ops->setup = KSPSetUp_GMRES;<br>
>> ksp->ops->solve = KSPSolve_GMRES;<br>
>> ksp->ops->reset = KSPReset_GMRES;<br>
>> ksp->ops->destroy = KSPDestroy_GMRES;<br>
>> ksp->ops->view = KSPView_GMRES;<br>
>> ksp->ops->setfromoptions = KSPSetFromOptions_GMRES;<br>
>> ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;<br>
>> ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;<br>
>> #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_HAVE_ESSL)<br>
>> ksp->ops->computeritz = KSPComputeRitz_GMRES;<br>
>> #endif<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);CHKERRQ(ierr);<br>
>> ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);CHKERRQ(ierr);<br>
>> <br>
>> gmres->haptol = 1.0e-30;<br>
>> gmres->q_preallocate = 0;<br>
>> gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;<br>
>> gmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;<br>
>> gmres->nrs = 0;<br>
>> gmres->sol_temp = 0;<br>
>> gmres->max_k = GMRES_DEFAULT_MAXK;<br>
>> gmres->Rsvd = 0;<br>
>> gmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;<br>
>> gmres->orthogwork = 0;<br>
>> gmres->delta = -1.0; // DRL<br>
>> PetscFunctionReturn(0);<br>
>> }<br>
<br>
</blockquote></div>