/*
|
/*
|
The ST (spectral transformation) interface routines related to the
|
The ST (spectral transformation) interface routines related to the
|
KSP object associated to it.
|
KSP object associated to it.
|
|
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
|
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
|
|
|
This file is part of SLEPc.
|
This file is part of SLEPc.
|
|
|
SLEPc is free software: you can redistribute it and/or modify it under the
|
SLEPc is free software: you can redistribute it and/or modify it under the
|
terms of version 3 of the GNU Lesser General Public License as published by
|
terms of version 3 of the GNU Lesser General Public License as published by
|
the Free Software Foundation.
|
the Free Software Foundation.
|
|
|
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
|
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
more details.
|
more details.
|
|
|
You should have received a copy of the GNU Lesser General Public License
|
You should have received a copy of the GNU Lesser General Public License
|
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
|
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
*/
|
*/
|
|
|
#include "private/stimpl.h" /*I "slepcst.h" I*/
|
#include <private/stimpl.h> /*I "slepcst.h" I*/
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STAssociatedKSPSolve"
|
#define __FUNCT__ "STAssociatedKSPSolve"
|
/*
|
/*
|
STAssociatedKSPSolve - Solves the linear system of equations associated
|
STAssociatedKSPSolve - Solves the linear system of equations associated
|
to the spectral transformation.
|
to the spectral transformation.
|
|
|
Input Parameters:
|
Input Parameters:
|
. st - the spectral transformation context
|
. st - the spectral transformation context
|
. b - right hand side vector
|
. b - right hand side vector
|
|
|
Output Parameter:
|
Output Parameter:
|
. x - computed solution
|
. x - computed solution
|
*/
|
*/
|
PetscErrorCode STAssociatedKSPSolve(ST st,Vec b,Vec x)
|
PetscErrorCode STAssociatedKSPSolve(ST st,Vec b,Vec x)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt its;
|
PetscInt its;
|
KSPConvergedReason reason;
|
KSPConvergedReason reason;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(b,VEC_CLASSID,2);
|
PetscValidHeaderSpecific(b,VEC_CLASSID,2);
|
PetscValidHeaderSpecific(x,VEC_CLASSID,3);
|
PetscValidHeaderSpecific(x,VEC_CLASSID,3);
|
if (!st->ksp) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST has no associated KSP"); }
|
if (!st->ksp) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST has no associated KSP"); }
|
ierr = KSPSolve(st->ksp,b,x);CHKERRQ(ierr);
|
ierr = KSPSolve(st->ksp,b,x);CHKERRQ(ierr);
|
ierr = KSPGetConvergedReason(st->ksp,&reason);CHKERRQ(ierr);
|
ierr = KSPGetConvergedReason(st->ksp,&reason);CHKERRQ(ierr);
|
if (reason<0) { SETERRQ1(((PetscObject)st)->comm,0,"Warning: KSP did not converge (%d)",reason); }
|
if (reason<0) { SETERRQ1(((PetscObject)st)->comm,0,"Warning: KSP did not converge (%d)",reason); }
|
ierr = KSPGetIterationNumber(st->ksp,&its);CHKERRQ(ierr);
|
ierr = KSPGetIterationNumber(st->ksp,&its);CHKERRQ(ierr);
|
st->lineariterations += its;
|
st->lineariterations += its;
|
PetscInfo1(st,"Linear solve iterations=%d\n",its);
|
PetscInfo1(st,"Linear solve iterations=%d\n",its);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STAssociatedKSPSolveTranspose"
|
#define __FUNCT__ "STAssociatedKSPSolveTranspose"
|
/*
|
/*
|
STAssociatedKSPSolveTranspose - Solves the transpose of the linear
|
STAssociatedKSPSolveTranspose - Solves the transpose of the linear
|
system of equations associated to the spectral transformation.
|
system of equations associated to the spectral transformation.
|
|
|
Input Parameters:
|
Input Parameters:
|
. st - the spectral transformation context
|
. st - the spectral transformation context
|
. b - right hand side vector
|
. b - right hand side vector
|
|
|
Output Parameter:
|
Output Parameter:
|
. x - computed solution
|
. x - computed solution
|
*/
|
*/
|
PetscErrorCode STAssociatedKSPSolveTranspose(ST st,Vec b,Vec x)
|
PetscErrorCode STAssociatedKSPSolveTranspose(ST st,Vec b,Vec x)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt its;
|
PetscInt its;
|
KSPConvergedReason reason;
|
KSPConvergedReason reason;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(b,VEC_CLASSID,2);
|
PetscValidHeaderSpecific(b,VEC_CLASSID,2);
|
PetscValidHeaderSpecific(x,VEC_CLASSID,3);
|
PetscValidHeaderSpecific(x,VEC_CLASSID,3);
|
if (!st->ksp) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST has no associated KSP"); }
|
if (!st->ksp) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST has no associated KSP"); }
|
ierr = KSPSolveTranspose(st->ksp,b,x);CHKERRQ(ierr);
|
ierr = KSPSolveTranspose(st->ksp,b,x);CHKERRQ(ierr);
|
ierr = KSPGetConvergedReason(st->ksp,&reason);CHKERRQ(ierr);
|
ierr = KSPGetConvergedReason(st->ksp,&reason);CHKERRQ(ierr);
|
if (reason<0) { SETERRQ1(((PetscObject)st)->comm,0,"Warning: KSP did not converge (%d)",reason); }
|
if (reason<0) { SETERRQ1(((PetscObject)st)->comm,0,"Warning: KSP did not converge (%d)",reason); }
|
ierr = KSPGetIterationNumber(st->ksp,&its);CHKERRQ(ierr);
|
ierr = KSPGetIterationNumber(st->ksp,&its);CHKERRQ(ierr);
|
st->lineariterations += its;
|
st->lineariterations += its;
|
PetscInfo1(st,"Linear solve iterations=%d\n",its);
|
PetscInfo1(st,"Linear solve iterations=%d\n",its);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STSetKSP"
|
#define __FUNCT__ "STSetKSP"
|
/*@
|
/*@
|
STSetKSP - Sets the KSP object associated with the spectral
|
STSetKSP - Sets the KSP object associated with the spectral
|
transformation.
|
transformation.
|
|
|
Not collective
|
Not collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ st - the spectral transformation context
|
+ st - the spectral transformation context
|
- ksp - the linear system context
|
- ksp - the linear system context
|
|
|
Level: advanced
|
Level: advanced
|
|
|
@*/
|
@*/
|
PetscErrorCode STSetKSP(ST st,KSP ksp)
|
PetscErrorCode STSetKSP(ST st,KSP ksp)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
|
PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
|
PetscCheckSameComm(st,1,ksp,2);
|
PetscCheckSameComm(st,1,ksp,2);
|
ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
|
ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
|
if (st->ksp) {
|
if (st->ksp) {
|
ierr = KSPDestroy(st->ksp);CHKERRQ(ierr);
|
ierr = KSPDestroy(st->ksp);CHKERRQ(ierr);
|
}
|
}
|
st->ksp = ksp;
|
st->ksp = ksp;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STGetKSP"
|
#define __FUNCT__ "STGetKSP"
|
/*@
|
/*@
|
STGetKSP - Gets the KSP object associated with the spectral
|
STGetKSP - Gets the KSP object associated with the spectral
|
transformation.
|
transformation.
|
|
|
Not collective
|
Not collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. st - the spectral transformation context
|
. st - the spectral transformation context
|
|
|
Output Parameter:
|
Output Parameter:
|
. ksp - the linear system context
|
. ksp - the linear system context
|
|
|
Notes:
|
Notes:
|
On output, the value of ksp can be PETSC_NULL if the combination of
|
On output, the value of ksp can be PETSC_NULL if the combination of
|
eigenproblem type and selected transformation does not require to
|
eigenproblem type and selected transformation does not require to
|
solve a linear system of equations.
|
solve a linear system of equations.
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
@*/
|
@*/
|
PetscErrorCode STGetKSP(ST st,KSP* ksp)
|
PetscErrorCode STGetKSP(ST st,KSP* ksp)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
if (!((PetscObject)st)->type_name) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call STSetType first"); }
|
if (!((PetscObject)st)->type_name) { SETERRQ(((PetscObject)st)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call STSetType first"); }
|
if (ksp) *ksp = st->ksp;
|
if (ksp) *ksp = st->ksp;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STGetOperationCounters"
|
#define __FUNCT__ "STGetOperationCounters"
|
/*@
|
/*@
|
STGetOperationCounters - Gets the total number of operator applications
|
STGetOperationCounters - Gets the total number of operator applications
|
and linear solver iterations used by the ST object.
|
and linear solver iterations used by the ST object.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameter:
|
Input Parameter:
|
. st - the spectral transformation context
|
. st - the spectral transformation context
|
|
|
Output Parameter:
|
Output Parameter:
|
+ ops - number of operator applications
|
+ ops - number of operator applications
|
- lits - number of linear solver iterations
|
- lits - number of linear solver iterations
|
|
|
Notes:
|
Notes:
|
Any output parameter may be PETSC_NULL on input if not needed.
|
Any output parameter may be PETSC_NULL on input if not needed.
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
.seealso: STResetOperationCounters()
|
.seealso: STResetOperationCounters()
|
@*/
|
@*/
|
PetscErrorCode STGetOperationCounters(ST st,PetscInt* ops,PetscInt* lits)
|
PetscErrorCode STGetOperationCounters(ST st,PetscInt* ops,PetscInt* lits)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
if (ops) *ops = st->applys;
|
if (ops) *ops = st->applys;
|
if (lits) *lits = st->lineariterations;
|
if (lits) *lits = st->lineariterations;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STResetOperationCounters"
|
#define __FUNCT__ "STResetOperationCounters"
|
/*@
|
/*@
|
STResetOperationCounters - Resets the counters for operator applications,
|
STResetOperationCounters - Resets the counters for operator applications,
|
inner product operations and total number of linear iterations used by
|
inner product operations and total number of linear iterations used by
|
the ST object.
|
the ST object.
|
|
|
Collective on ST
|
Collective on ST
|
|
|
Input Parameter:
|
Input Parameter:
|
. st - the spectral transformation context
|
. st - the spectral transformation context
|
|
|
Level: intermediate
|
Level: intermediate
|
|
|
.seealso: STGetOperationCounters()
|
.seealso: STGetOperationCounters()
|
@*/
|
@*/
|
PetscErrorCode STResetOperationCounters(ST st)
|
PetscErrorCode STResetOperationCounters(ST st)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
PetscValidHeaderSpecific(st,ST_CLASSID,1);
|
st->lineariterations = 0;
|
st->lineariterations = 0;
|
st->applys = 0;
|
st->applys = 0;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STCheckNullSpace_Default"
|
#define __FUNCT__ "STCheckNullSpace_Default"
|
PetscErrorCode STCheckNullSpace_Default(ST st,PetscInt n,const Vec V[])
|
PetscErrorCode STCheckNullSpace_Default(ST st,PetscInt n,const Vec V[])
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i,c;
|
PetscInt i,c;
|
PetscReal norm;
|
PetscReal norm;
|
Vec *T,w;
|
Vec *T,w;
|
Mat A;
|
Mat A;
|
PC pc;
|
PC pc;
|
MatNullSpace nullsp;
|
MatNullSpace nullsp;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscMalloc(n*sizeof(Vec),&T);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(Vec),&T);CHKERRQ(ierr);
|
ierr = KSPGetPC(st->ksp,&pc);CHKERRQ(ierr);
|
ierr = KSPGetPC(st->ksp,&pc);CHKERRQ(ierr);
|
ierr = PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
ierr = PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
ierr = MatGetVecs(A,PETSC_NULL,&w);CHKERRQ(ierr);
|
ierr = MatGetVecs(A,PETSC_NULL,&w);CHKERRQ(ierr);
|
c = 0;
|
c = 0;
|
for (i=0;i<n;i++) {
|
for (i=0;i<n;i++) {
|
ierr = MatMult(A,V[i],w);CHKERRQ(ierr);
|
ierr = MatMult(A,V[i],w);CHKERRQ(ierr);
|
ierr = VecNorm(w,NORM_2,&norm);CHKERRQ(ierr);
|
ierr = VecNorm(w,NORM_2,&norm);CHKERRQ(ierr);
|
if (norm < 1e-8) {
|
if (norm < 1e-8) {
|
PetscInfo2(st,"Vector %i norm=%g\n",i,norm);
|
PetscInfo2(st,"Vector %i norm=%g\n",i,norm);
|
T[c] = V[i];
|
T[c] = V[i];
|
c++;
|
c++;
|
}
|
}
|
}
|
}
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
ierr = VecDestroy(w);CHKERRQ(ierr);
|
if (c>0) {
|
if (c>0) {
|
ierr = MatNullSpaceCreate(((PetscObject)st)->comm,PETSC_FALSE,c,T,&nullsp);CHKERRQ(ierr);
|
ierr = MatNullSpaceCreate(((PetscObject)st)->comm,PETSC_FALSE,c,T,&nullsp);CHKERRQ(ierr);
|
ierr = KSPSetNullSpace(st->ksp,nullsp);CHKERRQ(ierr);
|
ierr = KSPSetNullSpace(st->ksp,nullsp);CHKERRQ(ierr);
|
ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr);
|
ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr);
|
}
|
}
|
ierr = PetscFree(T);CHKERRQ(ierr);
|
ierr = PetscFree(T);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "STCheckNullSpace"
|
#define __FUNCT__ "STCheckNullSpace"
|
/*@
|
/*@
|
STCheckNullSpace - Given a set of vectors, this function tests each of
|
STCheckNullSpace - Given a set of vectors, this function tests each of
|
them to be a nullspace vector of the coefficient matrix of the associated
|
them to be a nullspace vector of the coefficient matrix of the associated
|
KSP object. All these nullspace vectors are passed to the KSP object.
|
KSP object. All these nullspace vectors are passed to the KSP object.
|
|
|
Collective on ST
|
Collective on ST
|
|
|
Input Parameters:
|
Input Parameters:
|
+ st - the spectral transformation context
|
+ st - the spectral transformation context
|
. n - number of vectors
|
. n - number of vectors
|
- V - vectors to be checked
|
- V - vectors to be checked
|
|
|
Note:
|
Note:
|
This function allows to handle singular pencils and to solve some problems
|
This function allows to handle singular pencils and to solve some problems
|
in which the nullspace is important (see the users guide for details).
|
in which the nullspace is important (see the users guide for details).
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSSetDeflationSpace()
|
.seealso: EPSSetDeflationSpace()
|
@*/
|
@*/
|
PetscErrorCode STCheckNullSpace(ST st,PetscInt n,const Vec V[])
|
PetscErrorCode STCheckNullSpace(ST st,PetscInt n,const Vec V[])
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
if (n>0 && st->checknullspace) {
|
if (n>0 && st->checknullspace) {
|
ierr = (*st->checknullspace)(st,n,V);CHKERRQ(ierr);
|
ierr = (*st->checknullspace)(st,n,V);CHKERRQ(ierr);
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
|
|
|
|