[petsc-dev] Contribution proposal about the computation of (harmonic) Ritz pairs within GMRES
Sylvain Mercier
sylvainmercier85 at gmail.com
Tue Nov 3 10:09:47 CST 2015
Hi everyone,
During my phd thesis, I have worked on solving sequences of linear
systems with slowly varying matrices using GMRES(restart). In
particular, I have developed a preconditioning technique to improve
the action of an existing "first-level preconditioner". This new
method is defined using Ritz or harmonic Ritz vectors obtained at the
end of the solution of a linear system. Then I have developed two
routines in PETSc (that I've called KSPSetComputeRitz and
KSPComputeRitz) in order to compute the (harmonic) Ritz pairs
associated to the smallest or largest (harmonic) Ritz values in
modulus computed from the Hessenberg matrix of the last complete cycle
in GMRES.
Beyond the application to the preconditioning technique that I have
developed, this routine can be used to recover approximated eigenpairs
(and not only eigenvalues as already available in PETSc) of the
preconditioned matrix at the end of a solution with GMRES. That is why
I propose this contribution.
Two new routines has been developed, similarly to the existing
routines KSPSetComputeEigenvalues and KSPComputeEigenvalues.
The first one is called KSPSetComputeRitz and sets a flag so that the
last complete Hessenberg matrix computed with GMRES(restart) will be
stored. Here is the synopsis of the current version (located in
src/ksp/ksp/interface/itfunc.c)
PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
Input Parameters
ksp - iterative context obtained from KSPCreate()
flg - PETSC_TRUE or PETSC_FALSE
The second routine aims at computing the Ritz or harmonic Ritz
associated to the smallest or largest values in modulus. Here is the
synopsis of the current version (located in
src/ksp/ksp/impls/gmres/gmreig.c)
PetscErrorCode KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool
small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
Input Parameter
ksp - iterative context obtained from KSPCreate()
ritz - PETSC_TRUE or PETSC_FALSE to compute Ritz pairs
and harmonic Ritz pairs, respectively
small - PETSC_TRUE or PETSC_FALSE to compute pairs
associated to smallest or largest values in modulus, respectively
Input/Output parameter
n - The number of required and recovered pairs
Output Parameters
S[] - a multidimensional PETSc vector to store
the (harmonic) Ritz vectors, provided by user with a dimension of at
least n
tetar - real part of computed (harmonic) Ritz values, provided
by user with a dimension of at least n
tetai - imaginary part of computed (harmonic) Ritz values,
provided by user with a dimension of at least n
This routine has been developed whithin GMRES for real-valued linear
systems. Then the (harmonic) Ritz pairs are possibly complex-valued
and conjugated. In this case, two successive columns of S are equal to
the real and the imaginary parts of the vectors. Finally, the routine
KSPSolve_GMRES (located in src/ksp/ksp/impls/gmres/gmreig.c) has been
modified in order to store the Hessenberg matrix and the basis vectors
of the Krylov subspace as soon as a complete cycle has been performed.
To conclude, I propose to contribute to PETSc with these developments
(I have followed the conventions detailed in the developer guide).
Please find attached a patch and the modified files.
Regards,
Sylvain
-------------- next part --------------
From 94f3839b9290c48d39946edb8d9f1e510ed23d56 Mon Sep 17 00:00:00 2001
From: Sylvain Mercier <sylvainmercier85 at gmail.com>
Date: Sun, 1 Nov 2015 17:50:22 +0100
Subject: [PATCH] [PATCH] Computation of Ritz pairs within GMRES
Signed-off-by: Sylvain Mercier <sylvainmercier85 at gmail.com>
---
include/petsc/private/kspimpl.h | 2 +
include/petscksp.h | 2 +
src/ksp/ksp/impls/gmres/gmreig.c | 145 +++++++++++++++++++++++++++++++++++
src/ksp/ksp/impls/gmres/gmres.c | 28 ++++++-
src/ksp/ksp/impls/gmres/gmresimpl.h | 4 +
src/ksp/ksp/interface/itfunc.c | 84 ++++++++++++++++++++
6 files changed, 263 insertions(+), 2 deletions(-)
diff --git a/include/petsc/private/kspimpl.h b/include/petsc/private/kspimpl.h
index baf03fa..f4796e8 100644
--- a/include/petsc/private/kspimpl.h
+++ b/include/petsc/private/kspimpl.h
@@ -24,6 +24,7 @@ struct _KSPOps {
PetscErrorCode (*publishoptions)(KSP);
PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
+ PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
PetscErrorCode (*destroy)(KSP);
PetscErrorCode (*view)(KSP,PetscViewer);
PetscErrorCode (*reset)(KSP);
@@ -51,6 +52,7 @@ struct _p_KSP {
KSPFischerGuess guess;
PetscBool guess_zero, /* flag for whether initial guess is 0 */
calc_sings, /* calculate extreme Singular Values */
+ calc_ritz, /* calculate (harmonic) Ritz pairs */
guess_knoll; /* use initial guess of PCApply(ksp->B,b */
PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
diff --git a/include/petscksp.h b/include/petscksp.h
index e186c42..da70700 100644
--- a/include/petscksp.h
+++ b/include/petscksp.h
@@ -95,6 +95,7 @@ PETSC_EXTERN PetscErrorCode KSPGetInitialGuessKnoll(KSP,PetscBool *);
PETSC_EXTERN PetscErrorCode KSPSetErrorIfNotConverged(KSP,PetscBool );
PETSC_EXTERN PetscErrorCode KSPGetErrorIfNotConverged(KSP,PetscBool *);
PETSC_EXTERN PetscErrorCode KSPSetComputeEigenvalues(KSP,PetscBool );
+PETSC_EXTERN PetscErrorCode KSPSetComputeRitz(KSP,PetscBool );
PETSC_EXTERN PetscErrorCode KSPGetComputeEigenvalues(KSP,PetscBool *);
PETSC_EXTERN PetscErrorCode KSPSetComputeSingularValues(KSP,PetscBool );
PETSC_EXTERN PetscErrorCode KSPGetComputeSingularValues(KSP,PetscBool *);
@@ -148,6 +149,7 @@ PETSC_EXTERN PetscErrorCode KSPChebyshevEstEigGetKSP(KSP,KSP*);
PETSC_EXTERN PetscErrorCode KSPComputeExtremeSingularValues(KSP,PetscReal*,PetscReal*);
PETSC_EXTERN PetscErrorCode KSPComputeEigenvalues(KSP,PetscInt,PetscReal[],PetscReal[],PetscInt *);
PETSC_EXTERN PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP,PetscInt,PetscReal[],PetscReal[]);
+PETSC_EXTERN PetscErrorCode KSPComputeRitz(KSP,PetscBool,PetscBool,PetscInt *,Vec[],PetscReal[],PetscReal[]);
/*E
diff --git a/src/ksp/ksp/impls/gmres/gmreig.c b/src/ksp/ksp/impls/gmres/gmreig.c
index 44c3724..4d63bac 100644
--- a/src/ksp/ksp/impls/gmres/gmreig.c
+++ b/src/ksp/ksp/impls/gmres/gmreig.c
@@ -193,6 +193,151 @@ PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp,PetscInt nmax,PetscReal *r,Pe
PetscFunctionReturn(0);
}
+#undef __FUNCT__
+#define __FUNCT__ "KSPComputeRitz_GMRES"
+PetscErrorCode KSPComputeRitz_GMRES(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal *tetar,PetscReal *tetai)
+{
+ KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
+ PetscErrorCode ierr;
+ PetscInt n = gmres->it + 1,N = gmres->max_k + 1,NbrRitz,nb=0;
+ PetscInt i,j,*perm;
+ PetscReal *H,*Q,*Ht; /* H Hessenberg Matrix and Q matrix of eigenvectors of H*/
+ PetscReal *wr,*wi,*modul; /* Real and imaginary part and modul of the Ritz values*/
+ PetscReal *SR,*work;
+ PetscBLASInt bn,bN,lwork,idummy;
+ PetscScalar *t,sdummy;
+
+#if defined(PETSC_USE_COMPLEX)
+ SETERRQ(PetscObjectComm((PetscObject)ksp),-1,"NO SUPPORT FOR COMPLEX VALUES AT THIS TIME");
+#endif
+
+ PetscFunctionBegin;
+ /* n: size of the Hessenberg matrix */
+ if (gmres->fullcycle) n = N-1;
+ /* NbrRitz: number of (harmonic) Ritz pairs to extract */
+ NbrRitz = PetscMin(*nrit,n);
+
+ /* Definition of PetscBLASInt for lapack routines*/
+ ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr);
+ ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr);
+ ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr);
+ ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr);
+ /* Memory allocation */
+ ierr = PetscMalloc1(bN*bN,&H);CHKERRQ(ierr);
+ ierr = PetscMalloc1(bn*bn,&Q);CHKERRQ(ierr);
+ ierr = PetscMalloc1(lwork,&work);CHKERRQ(ierr);
+ ierr = PetscMalloc1(n,&wr);CHKERRQ(ierr);
+ ierr = PetscMalloc1(n,&wi);CHKERRQ(ierr);
+
+ /* copy H matrix to work space */
+ if (gmres->fullcycle) {
+ ierr = PetscMemcpy(H,gmres->hes_ritz,bN*bN*sizeof(PetscReal));CHKERRQ(ierr);
+ } else {
+ ierr = PetscMemcpy(H,gmres->hes_origin,bN*bN*sizeof(PetscReal));CHKERRQ(ierr);
+ }
+
+ /* Modify H to compute Harmonic Ritz pairs H = H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
+ if (!ritz) {
+ /* Transpose the Hessenberg matrix => Ht */
+ ierr = PetscMalloc1(bn*bn,&Ht);CHKERRQ(ierr);
+ for (i=0; i<bn; i++) {
+ for (j=0; j<bn; j++) {
+ Ht[i*bn+j] = H[j*bN+i];
+ }
+ }
+ /* Solve the system H^T*t = h^2_{m+1,m}e_m */
+ ierr = PetscCalloc1(bn,&t);CHKERRQ(ierr);
+ /* t = h^2_{m+1,m}e_m */
+ if (gmres->fullcycle) {
+ t[bn-1] = PetscSqr(gmres->hes_ritz[(bn-1)*bN+bn]);
+ } else {
+ t[bn-1] = PetscSqr(gmres->hes_origin[(bn-1)*bN+bn]);
+ }
+ /* Call the LAPACK routine dgesv to compute t = H^{-T}*t */
+#if defined(PETSC_MISSING_LAPACK_GESV)
+ SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GESV - Lapack routine is unavailable.");
+#else
+ {
+ PetscBLASInt info;
+ PetscBLASInt nrhs = 1;
+ PetscBLASInt *ipiv;
+ ierr = PetscMalloc1(bn,&ipiv);CHKERRQ(ierr);
+ PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&bn,&nrhs,Ht,&bn,ipiv,t,&bn,&info));
+ if (info) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Error while calling the Lapack routine DGESV");
+ ierr = PetscFree(ipiv);CHKERRQ(ierr);
+ ierr = PetscFree(Ht);CHKERRQ(ierr);
+ }
+#endif
+ /* Now form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
+ for (i=0; i<bn; i++) H[(bn-1)*bn+i] += t[i];
+ ierr = PetscFree(t);CHKERRQ(ierr);
+ }
+
+ /* Compute (harmonic) Ritz pairs */
+#if defined(PETSC_MISSING_LAPACK_HSEQR)
+ SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
+#else
+ {
+ PetscBLASInt info;
+ ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
+ PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&bn,H,&bN,wr,wi,&sdummy,&idummy,Q,&bn,work,&lwork,&info));
+ if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine");
+ }
+#endif
+ /* sort the (harmonic) Ritz values */
+ ierr = PetscMalloc1(n,&modul);CHKERRQ(ierr);
+ ierr = PetscMalloc1(n,&perm);CHKERRQ(ierr);
+ for (i=0; i<n; i++) modul[i] = PetscSqrtReal(wr[i]*wr[i]+wi[i]*wi[i]);
+ for (i=0; i<n; i++) perm[i] = i;
+ ierr = PetscSortRealWithPermutation(n,modul,perm);CHKERRQ(ierr);
+ /* count the number of extracted Ritz or Harmonic Ritz pairs (with complex conjugates) */
+ if (small) {
+ while (nb < NbrRitz) {
+ if (!wi[perm[nb]]) nb += 1;
+ else nb += 2;
+ }
+ ierr = PetscMalloc(nb*n*sizeof(PetscReal),&SR);
+ for (i=0; i<nb; i++) {
+ tetar[i] = wr[perm[i]];
+ tetai[i] = wi[perm[i]];
+ ierr = PetscMemcpy(&SR[i*n],&(Q[perm[i]*bn]),n*sizeof(PetscReal));CHKERRQ(ierr);
+ }
+ } else {
+ while (nb < NbrRitz) {
+ if (wi[perm[n-nb-1]] == 0) nb += 1;
+ else nb += 2;
+ }
+ ierr = PetscMalloc(nb*n*sizeof(PetscReal),&SR);
+ for (i=0; i<nb; i++) {
+ tetar[i] = wr[perm[n-nb+i]];
+ tetai[i] = wi[perm[n-nb+i]];
+ ierr = PetscMemcpy(&SR[i*n], &(Q[perm[n-nb+i]*bn]), n*sizeof(PetscReal));CHKERRQ(ierr);
+ }
+ }
+ ierr = PetscFree(modul);CHKERRQ(ierr);
+ ierr = PetscFree(perm);CHKERRQ(ierr);
+
+ /* Form the Ritz or Harmonic Ritz vectors S=VV*Sr,
+ where the columns of VV correspond to the basis of the Krylov subspace */
+ if (gmres->fullcycle) {
+ for (j=0; j<nb; j++) {
+ ierr = VecZeroEntries(S[j]);CHKERRQ(ierr);
+ ierr = VecMAXPY(S[j],n,&SR[j*n],gmres->vecb);CHKERRQ(ierr);
+ }
+ } else {
+ for (j=0; j<nb; j++) {
+ ierr = VecZeroEntries(S[j]);CHKERRQ(ierr);
+ ierr = VecMAXPY(S[j],n,&SR[j*n],&VEC_VV(0));CHKERRQ(ierr);
+ }
+ }
+ *nrit = nb;
+ ierr = PetscFree(H);CHKERRQ(ierr);
+ ierr = PetscFree(Q);CHKERRQ(ierr);
+ ierr = PetscFree(SR);CHKERRQ(ierr);
+ ierr = PetscFree(wr);CHKERRQ(ierr);
+ ierr = PetscFree(wi);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
diff --git a/src/ksp/ksp/impls/gmres/gmres.c b/src/ksp/ksp/impls/gmres/gmres.c
index 1e9c2aa..6da6574 100644
--- a/src/ksp/ksp/impls/gmres/gmres.c
+++ b/src/ksp/ksp/impls/gmres/gmres.c
@@ -219,9 +219,11 @@ PetscErrorCode KSPGMRESCycle(PetscInt *itcount,KSP ksp)
PetscErrorCode KSPSolve_GMRES(KSP ksp)
{
PetscErrorCode ierr;
- PetscInt its,itcount;
+ PetscInt its,itcount,i;
KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
PetscBool guess_zero = ksp->guess_zero;
+ PetscInt N = gmres->max_k + 1;
+ PetscBLASInt bN;
PetscFunctionBegin;
if (ksp->calc_sings && !gmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
@@ -231,10 +233,28 @@ PetscErrorCode KSPSolve_GMRES(KSP ksp)
ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
itcount = 0;
+ gmres->fullcycle = 0;
ksp->reason = KSP_CONVERGED_ITERATING;
while (!ksp->reason) {
ierr = KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);CHKERRQ(ierr);
ierr = KSPGMRESCycle(&its,ksp);CHKERRQ(ierr);
+ /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace
+ if the cycle is complete for the computation of the Ritz pairs */
+ if (its == gmres->max_k) {
+ gmres->fullcycle++;
+ if (ksp->calc_ritz) {
+ if (!gmres->hes_ritz) {
+ ierr = PetscMalloc1(N*N,&gmres->hes_ritz);CHKERRQ(ierr);
+ ierr = PetscLogObjectMemory((PetscObject)ksp,N*N*sizeof(PetscScalar));CHKERRQ(ierr);
+ ierr = VecDuplicateVecs(VEC_VV(0),N,&gmres->vecb);CHKERRQ(ierr);
+ }
+ ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr);
+ ierr = PetscMemcpy(gmres->hes_ritz,gmres->hes_origin,bN*bN*sizeof(PetscReal));CHKERRQ(ierr);
+ for (i=0; i<gmres->max_k+1; i++) {
+ ierr = VecCopy(VEC_VV(i),gmres->vecb[i]);CHKERRQ(ierr);
+ }
+ }
+ }
itcount += its;
if (itcount >= ksp->max_it) {
if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
@@ -256,7 +276,7 @@ PetscErrorCode KSPReset_GMRES(KSP ksp)
PetscFunctionBegin;
/* Free the Hessenberg matrices */
- ierr = PetscFree5(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin);CHKERRQ(ierr);
+ ierr = PetscFree6(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin,gmres->hes_ritz);CHKERRQ(ierr);
/* free work vectors */
ierr = PetscFree(gmres->vecs);CHKERRQ(ierr);
@@ -264,6 +284,9 @@ PetscErrorCode KSPReset_GMRES(KSP ksp)
ierr = VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);CHKERRQ(ierr);
}
gmres->nwork_alloc = 0;
+ if (gmres->vecb) {
+ ierr = VecDestroyVecs(gmres->max_k+1,&gmres->vecb);CHKERRQ(ierr);
+ }
ierr = PetscFree(gmres->user_work);CHKERRQ(ierr);
ierr = PetscFree(gmres->mwork_alloc);CHKERRQ(ierr);
@@ -913,6 +936,7 @@ PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
ksp->ops->setfromoptions = KSPSetFromOptions_GMRES;
ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
+ ksp->ops->computeritz = KSPComputeRitz_GMRES;
ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);CHKERRQ(ierr);
ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);CHKERRQ(ierr);
diff --git a/src/ksp/ksp/impls/gmres/gmresimpl.h b/src/ksp/ksp/impls/gmres/gmresimpl.h
index b52c098..378cc74 100644
--- a/src/ksp/ksp/impls/gmres/gmresimpl.h
+++ b/src/ksp/ksp/impls/gmres/gmresimpl.h
@@ -12,6 +12,7 @@
/* Hessenberg matrix and orthogonalization information. */ \
PetscScalar *hh_origin; /* holds hessenburg matrix that has been multiplied by plane rotations (upper tri) */ \
PetscScalar *hes_origin; /* holds the original (unmodified) hessenberg matrix which may be used to estimate the Singular Values of the matrix */ \
+ PetscScalar *hes_ritz; /* holds the last full Hessenberg matrix to compute (harmonic) Ritz pairs */ \
PetscScalar *cc_origin; /* holds cosines for rotation matrices */ \
PetscScalar *ss_origin; /* holds sines for rotation matrices */ \
PetscScalar *rs_origin; /* holds the right-hand-side of the Hessenberg system */ \
@@ -31,6 +32,7 @@
KSPGMRESCGSRefinementType cgstype; \
\
Vec *vecs; /* the work vectors */ \
+ Vec *vecb; /* holds the last full basis vectors of the Krylov subspace to compute (harmonic) Ritz pairs */ \
PetscInt q_preallocate; /* 0=don't preallocate space for work vectors */ \
PetscInt delta_allocate; /* number of vectors to preallocaate in each block if not preallocated */ \
PetscInt vv_allocated; /* number of allocated gmres direction vectors */ \
@@ -42,6 +44,7 @@
\
/* Information for building solution */ \
PetscInt it; /* Current iteration: inside restart */ \
+ PetscInt fullcycle; /* Current number of complete cycle */ \
PetscScalar *nrs; /* temp that holds the coefficients of the Krylov vectors that form the minimum residual solution */ \
Vec sol_temp; /* used to hold temporary solution */
@@ -54,6 +57,7 @@ PETSC_INTERN PetscErrorCode KSPSetUp_GMRES(KSP);
PETSC_INTERN PetscErrorCode KSPSetFromOptions_GMRES(PetscOptions *PetscOptionsObject,KSP);
PETSC_INTERN PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP,PetscReal*,PetscReal*);
PETSC_INTERN PetscErrorCode KSPComputeEigenvalues_GMRES(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt*);
+PETSC_INTERN PetscErrorCode KSPComputeRitz_GMRES(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
PETSC_INTERN PetscErrorCode KSPReset_GMRES(KSP);
PETSC_INTERN PetscErrorCode KSPDestroy_GMRES(KSP);
PETSC_INTERN PetscErrorCode KSPGMRESGetNewVectors(KSP,PetscInt);
diff --git a/src/ksp/ksp/interface/itfunc.c b/src/ksp/ksp/interface/itfunc.c
index 72eeece..af28dea 100644
--- a/src/ksp/ksp/interface/itfunc.c
+++ b/src/ksp/ksp/interface/itfunc.c
@@ -131,6 +131,59 @@ PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal
}
#undef __FUNCT__
+#define __FUNCT__ "KSPComputeRitz"
+/*@
+ KSPComputeRitz - Computes the Ritz or harmonic Ritz pairs associated to the
+ smallest or largest in modulus, for the preconditioned operator.
+ Called after KSPSolve().
+
+ Not Collective
+
+ Input Parameter:
++ ksp - iterative context obtained from KSPCreate()
+. ritz - PETSC_TRUE or PETSC_FALSE for ritz pairs or harmonic Ritz pairs, respectively
+. small - PETSC_TRUE or PETSC_FALSE for smallest or largest (harmonic) Ritz values, respectively
+. nrit - number of (harmonic) Ritz pairs to compute
+
+ Output Parameters:
++ nrit - actual number of computed (harmonic) Ritz pairs
+. S - multidimensional vector with Ritz vectors
+. tetar - real part of the Ritz values
+. tetai - imaginary part of the Ritz values
+
+ Notes:
+ -For GMRES, the (harmonic) Ritz pairs are computed from the Hessenberg matrix obtained during
+ the last complete cycle, or obtained at the end of the solution if the method is stopped before
+ a restart. Then, the number of actual (harmonic) Ritz pairs computed is less or equal to the restart
+ parameter for GMRES if a complete cycle has been performed or less or equal to the number of GMRES
+ iterations.
+ -Moreover, for real matrices, the (harmonic) Ritz pairs are possibly complex-valued. In such a case,
+ the routine selects the complex (harmonic) Ritz value and its conjugate, and two successive columns of S
+ are equal to the real and the imaginary parts of the associated vectors.
+ -the (harmonic) Ritz pairs are given in order of increasing (harmonic) Ritz values in modulus
+
+
+ One must call KSPSetComputeRitz() before calling KSPSetUp()
+ in order for this routine to work correctly.
+
+ Level: advanced
+
+.keywords: KSP, compute, ritz, values
+
+.seealso: KSPSetComputeRitz()
+@*/
+PetscErrorCode KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
+{
+ PetscErrorCode ierr;
+
+ PetscFunctionBegin;
+ PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
+ if (!ksp->calc_ritz) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Ritz pairs not requested before KSPSetUp()");
+
+ if (ksp->ops->computeritz) {ierr = (*ksp->ops->computeritz)(ksp,ritz,small,nrit,S,tetar,tetai);CHKERRQ(ierr);}
+ PetscFunctionReturn(0);
+}
+#undef __FUNCT__
#define __FUNCT__ "KSPSetUpOnBlocks"
/*@
KSPSetUpOnBlocks - Sets up the preconditioner for each block in
@@ -1434,6 +1487,37 @@ PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
}
#undef __FUNCT__
+#define __FUNCT__ "KSPSetComputeRitz"
+/*@
+ KSPSetComputeRitz - Sets a flag so that the ritz or harmonic ritz pairs
+ will be calculated via a Lanczos or Arnoldi process as the linear
+ system is solved.
+
+ Logically Collective on KSP
+
+ Input Parameters:
++ ksp - iterative context obtained from KSPCreate()
+- flg - PETSC_TRUE or PETSC_FALSE
+
+ Notes:
+ Currently this option is only valid for the GMRES method.
+
+ Level: advanced
+
+.keywords: KSP, set, compute, ritz
+
+.seealso: KSPComputeRitz()
+@*/
+PetscErrorCode KSPSetComputeRitz(KSP ksp, PetscBool flg)
+{
+ PetscFunctionBegin;
+ PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
+ PetscValidLogicalCollectiveBool(ksp,flg,2);
+ ksp->calc_ritz = flg;
+ PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
#define __FUNCT__ "KSPGetRhs"
/*@
KSPGetRhs - Gets the right-hand-side vector for the linear system to
--
1.7.9.5
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Ritz_routines.tar.gz
Type: application/x-gzip
Size: 39447 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20151103/0a7d51da/attachment.gz>
More information about the petsc-dev
mailing list