/*
The ST (spectral transformation) interface routines, callable by users.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
This file is part of SLEPc.
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
the Free Software Foundation.
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
#include <private/stimpl.h> /*I "slepcst.h" I*/
#undef __FUNCT__
#define __FUNCT__ "STApply"
/*@
STApply - Applies the spectral transformation operator to a vector, for
instance (A - sB)^-1 B in the case of the shift-and-invert tranformation
and generalized eigenproblem.
Collective on ST and Vec
Input Parameters:
+ st - the spectral transformation context
- x - input vector
Output Parameter:
. y - output vector
Level: developer
.seealso: STApplyTranspose()
@*/
PetscErrorCode STApply(ST st,Vec x,Vec y)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidHeaderSpecific(x,VEC_CLASSID,2);
PetscValidHeaderSpecific(y,VEC_CLASSID,3);
if (x == y) SETERRQ(((PetscObject)st)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
if (!st->setupcalled) { ierr = STSetUp(st);CHKERRQ(ierr); }
if (!st->ops->apply) SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST does not have apply");
ierr = PetscLogEventBegin(ST_Apply,st,x,y,0);CHKERRQ(ierr);
st->applys++;
if (st->D) { /* with balancing */
ierr = VecPointwiseDivide(st->wb,x,st->D);CHKERRQ(ierr);
ierr = (*st->ops->apply)(st,st->wb,y);CHKERRQ(ierr);
ierr = VecPointwiseMult(y,y,st->D);CHKERRQ(ierr);
}
else {
ierr = (*st->ops->apply)(st,x,y);CHKERRQ(ierr);
}
ierr = PetscLogEventEnd(ST_Apply,st,x,y,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetBilinearForm"
/*@
STGetBilinearForm - Returns the matrix used in the bilinear form with a
generalized problem with semi-definite B.
Not collective, though a parallel Mat may be returned
Input Parameters:
. st - the spectral transformation context
Output Parameter:
. B - output matrix
Notes:
The output matrix B must be destroyed after use. It will be PETSC_NULL in
case of standard eigenproblems.
Level: developer
@*/
PetscErrorCode STGetBilinearForm(ST st,Mat *B)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidPointer(B,2);
ierr = (*st->ops->getbilinearform)(st,B);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetBilinearForm_Default"
PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
{
PetscErrorCode ierr;
PetscFunctionBegin;
*B = st->B;
if (*B) {
ierr = PetscObjectReference((PetscObject)*B);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STApplyTranspose"
/*@
STApplyTranspose - Applies the transpose of the operator to a vector, for
instance B^T(A - sB)^-T in the case of the shift-and-invert tranformation
and generalized eigenproblem.
Collective on ST and Vec
Input Parameters:
+ st - the spectral transformation context
- x - input vector
Output Parameter:
. y - output vector
Level: developer
.seealso: STApply()
@*/
PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidHeaderSpecific(x,VEC_CLASSID,2);
PetscValidHeaderSpecific(y,VEC_CLASSID,3);
if (x == y) SETERRQ(((PetscObject)st)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
if (!st->setupcalled) { ierr = STSetUp(st);CHKERRQ(ierr); }
if (!st->ops->applytrans) SETERRQ(((PetscObject)st)->comm,PETSC_ERR_SUP,"ST does not have applytrans");
ierr = PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);CHKERRQ(ierr);
st->applys++;
if (st->D) { /* with balancing */
ierr = VecPointwiseMult(st->wb,x,st->D);CHKERRQ(ierr);
ierr = (*st->ops->applytrans)(st,st->wb,y);CHKERRQ(ierr);
ierr = VecPointwiseDivide(y,y,st->D);CHKERRQ(ierr);
}
else {
ierr = (*st->ops->applytrans)(st,x,y);CHKERRQ(ierr);
}
ierr = PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STComputeExplicitOperator"
/*@
STComputeExplicitOperator - Computes the explicit operator associated
to the eigenvalue problem with the specified spectral transformation.
Collective on ST
Input Parameter:
. st - the spectral transform context
Output Parameter:
. mat - the explicit operator
Notes:
This routine builds a matrix containing the explicit operator. For
example, in generalized problems with shift-and-invert spectral
transformation the result would be matrix (A - s B)^-1 B.
This computation is done by applying the operator to columns of the
identity matrix. Note that the result is a dense matrix.
Level: advanced
.seealso: STApply()
@*/
PetscErrorCode STComputeExplicitOperator(ST st,Mat *mat)
{
PetscErrorCode ierr;
Vec in,out;
PetscInt i,M,m,*rows,start,end;
const PetscScalar *array;
PetscScalar one = 1.0;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidPointer(mat,2);
ierr = MatGetVecs(st->A,&in,&out);CHKERRQ(ierr);
ierr = VecGetSize(out,&M);CHKERRQ(ierr);
ierr = VecGetLocalSize(out,&m);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(out,&start,&end);CHKERRQ(ierr);
ierr = PetscMalloc(m*sizeof(PetscInt),&rows);CHKERRQ(ierr);
for (i=0; i<m; i++) rows[i] = start + i;
ierr = MatCreateMPIDense(((PetscObject)st)->comm,m,m,M,M,PETSC_NULL,mat);CHKERRQ(ierr);
for (i=0; i<M; i++) {
ierr = VecSet(in,0.0);CHKERRQ(ierr);
ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
ierr = VecAssemblyEnd(in);CHKERRQ(ierr);
ierr = STApply(st,in,out);CHKERRQ(ierr);
ierr = VecGetArrayRead(out,&array);CHKERRQ(ierr);
ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(out,&array);CHKERRQ(ierr);
}
ierr = PetscFree(rows);CHKERRQ(ierr);
ierr = VecDestroy(&in);CHKERRQ(ierr);
ierr = VecDestroy(&out);CHKERRQ(ierr);
ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetUp"
/*@
STSetUp - Prepares for the use of a spectral transformation.
Collective on ST
Input Parameter:
. st - the spectral transformation context
Level: advanced
.seealso: STCreate(), STApply(), STDestroy()
@*/
PetscErrorCode STSetUp(ST st)
{
PetscInt n,k;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (!st->A) {SETERRQ(((PetscObject)st)->comm,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}
if (st->setupcalled) PetscFunctionReturn(0);
PetscInfo(st,"Setting up new ST\n");
ierr = PetscLogEventBegin(ST_SetUp,st,0,0,0);CHKERRQ(ierr);
if (!((PetscObject)st)->type_name) {
ierr = STSetType(st,STSHIFT);CHKERRQ(ierr);
}
ierr = VecDestroy(&st->w);CHKERRQ(ierr);
ierr = MatGetVecs(st->A,&st->w,PETSC_NULL);CHKERRQ(ierr);
if (st->D) {
ierr = MatGetLocalSize(st->A,PETSC_NULL,&n);CHKERRQ(ierr);
ierr = VecGetLocalSize(st->D,&k);CHKERRQ(ierr);
if (n != k) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %D (should be %D)",k,n);
ierr = VecDestroy(&st->wb);CHKERRQ(ierr);
ierr = VecDuplicate(st->D,&st->wb);CHKERRQ(ierr);
}
if (st->ops->setup) { ierr = (*st->ops->setup)(st);CHKERRQ(ierr); }
st->setupcalled = 1;
ierr = PetscLogEventEnd(ST_SetUp,st,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STPostSolve"
/*@
STPostSolve - Optional post-solve phase, intended for any actions that must
be performed on the ST object after the eigensolver has finished.
Collective on ST
Input Parameters:
. st - the spectral transformation context
Level: developer
.seealso: EPSSolve()
@*/
PetscErrorCode STPostSolve(ST st)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (st->ops->postsolve) {
ierr = (*st->ops->postsolve)(st);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STBackTransform"
/*@
STBackTransform - Back-transformation phase, intended for
spectral transformations which require to transform the computed
eigenvalues back to the original eigenvalue problem.
Not Collective
Input Parameters:
st - the spectral transformation context
eigr - real part of a computed eigenvalue
eigi - imaginary part of a computed eigenvalue
Level: developer
@*/
PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (st->ops->backtr) {
ierr = (*st->ops->backtr)(st,n,eigr,eigi);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}