static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
{
PetscScalar *hh,*cc,*ss,tt;
PetscInt j;
KSP_LGMRES *lgmres = (KSP_LGMRES *)(ksp->data);
PetscFunctionBegin;
hh = HH(0,it); /* pointer to beginning of column to update - so
incrementing hh "steps down" the (it+1)th col of HH*/
cc = CC(0); /* beginning of cosine rotations */
ss = SS(0); /* beginning of sine rotations */
/* Apply all the previously computed plane rotations to the new column
of the Hessenberg matrix */
/* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta) */
for (j=1; j<=it; j++) {
tt = *hh;
*hh = PetscConj(*cc) * tt + *ss * *(hh+1);
hh++;
*hh = *cc++ * *hh - (*ss++ * tt);
/* hh, cc, and ss have all been incremented one by end of loop */
}
/*
compute the new plane rotation, and apply it to:
1) the right-hand-side of the Hessenberg system (GRS)
note: it affects GRS(it) and GRS(it+1)
2) the new column of the Hessenberg matrix
note: it affects HH(it,it) which is currently pointed to
by hh and HH(it+1, it) (*(hh+1))
thus obtaining the updated value of the residual...
*/
/* compute new plane rotation */
if (!hapend) {
tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
if (tt == 0.0) {
ksp->reason = KSP_DIVERGED_NULL;
PetscFunctionReturn(0);
}
*cc = *hh / tt; /* new cosine value */
*ss = *(hh+1) / tt; /* new sine value */
/* apply to 1) and 2) */
*GRS(it+1) = - (*ss * *GRS(it));
*GRS(it) = PetscConj(*cc) * *GRS(it);
*hh = PetscConj(*cc) * *hh + *ss * *(hh+1);
/* residual is the last element (it+1) of right-hand side! */
*res = PetscAbsScalar(*GRS(it+1));
} else {
/* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
another rotation matrix (so RH doesn't change). The new residual is
always the new sine term times the residual from last time (GRS(it)),
but now the new sine rotation would be zero...so the residual should
be zero...so we will multiply "zero" by the last residual. This might
not be exactly what we want to do here -could just return "zero". */
*res = 0.0;
}
PetscFunctionReturn(0);
}
/*
. it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this
*/
static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool *hapend,PetscReal *res)
{
PetscScalar *hh,*cc,*ss,*rs;
PetscInt j;
PetscReal hapbnd;
KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data);
PetscErrorCode ierr;
PetscFunctionBegin;
hh = HH(0,it); /* pointer to beginning of column to update */
cc = CC(0); /* beginning of cosine rotations */
ss = SS(0); /* beginning of sine rotations */
rs = RS(0); /* right hand side of least squares system */
/* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
for (j=0; j<=it+1; j++) *HES(j,it) = hh[j];
/* check for the happy breakdown */
hapbnd = PetscMin(PetscAbsScalar(hh[it+1] / rs[it]),pgmres->haptol);
if (PetscAbsScalar(hh[it+1]) < hapbnd) {
ierr = PetscInfo4(ksp,"Detected happy breakdown, current hapbnd = %14.12e H(%D,%D) = %14.12e\n",(double)hapbnd,it+1,it,(double)PetscAbsScalar(*HH(it+1,it)));CHKERRQ(ierr);
*hapend = PETSC_TRUE;
}
/* Apply all the previously computed plane rotations to the new column
of the Hessenberg matrix */
/* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta),
and some refs have [c s ; -conj(s) c] (don't be confused!) */
for (j=0; j<it; j++) {
PetscScalar hhj = hh[j];
hh[j] = PetscConj(cc[j])*hhj + ss[j]*hh[j+1];
hh[j+1] = -ss[j] *hhj + cc[j]*hh[j+1];
}
/*
compute the new plane rotation, and apply it to:
1) the right-hand-side of the Hessenberg system (RS)
note: it affects RS(it) and RS(it+1)
2) the new column of the Hessenberg matrix
note: it affects HH(it,it) which is currently pointed to
by hh and HH(it+1, it) (*(hh+1))
thus obtaining the updated value of the residual...
*/
/* compute new plane rotation */
if (!*hapend) {
PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it+1])));
if (delta == 0.0) {
ksp->reason = KSP_DIVERGED_NULL;
PetscFunctionReturn(0);
}
cc[it] = hh[it] / delta; /* new cosine value */
ss[it] = hh[it+1] / delta; /* new sine value */
hh[it] = PetscConj(cc[it])*hh[it] + ss[it]*hh[it+1];
rs[it+1] = -ss[it]*rs[it];
rs[it] = PetscConj(cc[it])*rs[it];
*res = PetscAbsScalar(rs[it+1]);
} else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
another rotation matrix (so RH doesn't change). The new residual is
always the new sine term times the residual from last time (RS(it)),
but now the new sine rotation would be zero...so the residual should
be zero...so we will multiply "zero" by the last residual. This might
not be exactly what we want to do here -could just return "zero". */
*res = 0.0;
}
PetscFunctionReturn(0);
}
//.........这里部分代码省略.........
/*LGMRES_MOD: decide whether this is an arnoldi step or an aug step */
if (loc_it < it_arnoldi) { /* Arnoldi */
ierr = KSP_PCApplyBAorAB(ksp,VEC_VV(loc_it),VEC_VV(1+loc_it),VEC_TEMP_MATOP);
CHKERRQ(ierr);
} else { /*aug step */
order = loc_it - it_arnoldi + 1; /* which aug step */
for (ii=0; ii<aug_dim; ii++) {
if (lgmres->aug_order[ii] == order) {
spot = ii;
break; /* must have this because there will be duplicates before aug_ct = aug_dim */
}
}
ierr = VecCopy(A_AUGVEC(spot), VEC_VV(1+loc_it));
CHKERRQ(ierr);
/*note: an alternate implementation choice would be to only save the AUGVECS and
not A_AUGVEC and then apply the PC here to the augvec */
}
/* update hessenberg matrix and do Gram-Schmidt - new direction is in
VEC_VV(1+loc_it)*/
ierr = (*lgmres->orthog)(ksp,loc_it);
CHKERRQ(ierr);
/* new entry in hessenburg is the 2-norm of our new direction */
ierr = VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);
CHKERRQ(ierr);
*HH(loc_it+1,loc_it) = tt;
*HES(loc_it+1,loc_it) = tt;
/* check for the happy breakdown */
hapbnd = PetscAbsScalar(tt / *GRS(loc_it));/* GRS(loc_it) contains the res_norm from the last iteration */
if (hapbnd > lgmres->haptol) hapbnd = lgmres->haptol;
if (tt > hapbnd) {
tmp = 1.0/tt;
ierr = VecScale(VEC_VV(loc_it+1),tmp);
CHKERRQ(ierr); /* scale new direction by its norm */
} else {
ierr = PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %G tt = %G\n",hapbnd,tt);
CHKERRQ(ierr);
hapend = PETSC_TRUE;
}
/* Now apply rotations to new col of hessenberg (and right side of system),
calculate new rotation, and get new residual norm at the same time*/
ierr = KSPLGMRESUpdateHessenberg(ksp,loc_it,hapend,&res);
CHKERRQ(ierr);
if (ksp->reason) break;
loc_it++;
lgmres->it = (loc_it-1); /* Add this here in case it has converged */
ierr = PetscObjectTakeAccess(ksp);
CHKERRQ(ierr);
ksp->its++;
ksp->rnorm = res;
ierr = PetscObjectGrantAccess(ksp);
CHKERRQ(ierr);
ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
CHKERRQ(ierr);
/* Catch error in happy breakdown and signal convergence and break from loop */
if (hapend) {
PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
{
PetscErrorCode ierr;
PC_IS *pcis = (PC_IS*)(pc->data);
PetscFunctionBegin;
/*
Neumann solvers.
Applying the inverse of the local Schur complement, i.e, solving a Neumann
Problem with zero at the interior nodes of the RHS and extracting the interface
part of the solution. inverse Schur complement is applied to b and the result
is stored in x.
*/
/* Setting the RHS vec1_N */
ierr = VecSet(vec1_N,0.0);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
/* Checking for consistency of the RHS */
{
PetscBool flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(NULL,"-pc_is_check_consistency",&flg,NULL);CHKERRQ(ierr);
if (flg) {
PetscScalar average;
PetscViewer viewer;
ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);CHKERRQ(ierr);
ierr = VecSum(vec1_N,&average);CHKERRQ(ierr);
average = average / ((PetscReal)pcis->n);
ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
if (pcis->pure_neumann) {
ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
}
ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
}
}
/* Solving the system for vec2_N */
ierr = KSPSolve(pcis->ksp_N,vec1_N,vec2_N);CHKERRQ(ierr);
/* Extracting the local interface vector out of the solution */
ierr = VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar* nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
{
PetscScalar tt;
PetscErrorCode ierr;
PetscInt ii,k,j;
KSP_LGMRES *lgmres = (KSP_LGMRES *)(ksp->data);
/*LGMRES_MOD */
PetscInt it_arnoldi, it_aug;
PetscInt jj, spot = 0;
PetscFunctionBegin;
/* Solve for solution vector that minimizes the residual */
/* If it is < 0, no lgmres steps have been performed */
if (it < 0) {
ierr = VecCopy(vguess,vdest);
CHKERRQ(ierr); /* VecCopy() is smart, exists immediately if vguess == vdest */
PetscFunctionReturn(0);
}
/* so (it+1) lgmres steps HAVE been performed */
/* LGMRES_MOD - determine if we need to use augvecs for the soln - do not assume that
this is called after the total its allowed for an approx space */
if (lgmres->approx_constant) {
it_arnoldi = lgmres->max_k - lgmres->aug_ct;
} else {
it_arnoldi = lgmres->max_k - lgmres->aug_dim;
}
if (it_arnoldi >= it +1) {
it_aug = 0;
it_arnoldi = it+1;
} else {
it_aug = (it + 1) - it_arnoldi;
}
/* now it_arnoldi indicates the number of matvecs that took place */
lgmres->matvecs += it_arnoldi;
/* solve the upper triangular system - GRS is the right side and HH is
the upper triangular matrix - put soln in nrs */
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,PetscAbsScalar(*GRS(it)));
if (*HH(it,it) != 0.0) {
nrs[it] = *GRS(it) / *HH(it,it);
} else {
nrs[it] = 0.0;
}
for (ii=1; ii<=it; ii++) {
k = it - ii;
tt = *GRS(k);
for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
nrs[k] = tt / *HH(k,k);
}
/* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */
ierr = VecSet(VEC_TEMP,0.0);
CHKERRQ(ierr); /* set VEC_TEMP components to 0 */
/*LGMRES_MOD - if augmenting has happened we need to form the solution
using the augvecs */
if (!it_aug) { /* all its are from arnoldi */
ierr = VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
CHKERRQ(ierr);
} else { /*use aug vecs */
/*first do regular krylov directions */
ierr = VecMAXPY(VEC_TEMP,it_arnoldi,nrs,&VEC_VV(0));
CHKERRQ(ierr);
/*now add augmented portions - add contribution of aug vectors one at a time*/
for (ii=0; ii<it_aug; ii++) {
for (jj=0; jj<lgmres->aug_dim; jj++) {
if (lgmres->aug_order[jj] == (ii+1)) {
spot = jj;
break; /* must have this because there will be duplicates before aug_ct = aug_dim */
}
}
ierr = VecAXPY(VEC_TEMP,nrs[it_arnoldi+ii],AUGVEC(spot));
CHKERRQ(ierr);
}
}
/* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
preconditioner is "unwound" from right-precondtioning*/
ierr = VecCopy(VEC_TEMP, AUG_TEMP);
CHKERRQ(ierr);
ierr = KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
CHKERRQ(ierr);
/* add solution to previous solution */
/* put updated solution into vdest.*/
ierr = VecCopy(vguess,vdest);
CHKERRQ(ierr);
ierr = VecAXPY(vdest,1.0,VEC_TEMP);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}
客服电话
APP下载
官方微信
请发表评论