/*
The ST (spectral transformation) interface routines, callable by users.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2010, Universidad 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*/
PetscClassId ST_CLASSID = 0;
PetscLogEvent ST_SetUp = 0, ST_Apply = 0, ST_ApplyTranspose = 0;
static PetscTruth STPackageInitialized = PETSC_FALSE;
#undef __FUNCT__
#define __FUNCT__ "STFinalizePackage"
/*@C
STFinalizePackage - This function destroys everything in the Slepc interface to the ST package. It is
called from SlepcFinalize().
Level: developer
.seealso: SlepcFinalize()
@*/
PetscErrorCode STFinalizePackage(void)
{
PetscFunctionBegin;
STPackageInitialized = PETSC_FALSE;
STList = 0;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STInitializePackage"
/*@C
STInitializePackage - This function initializes everything in the ST package. It is called
from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to STCreate()
when using static libraries.
Input Parameter:
path - The dynamic library path, or PETSC_NULL
Level: developer
.seealso: SlepcInitialize()
@*/
PetscErrorCode STInitializePackage(const char *path) {
char logList[256];
char *className;
PetscTruth opt;
PetscErrorCode ierr;
PetscFunctionBegin;
if (STPackageInitialized) PetscFunctionReturn(0);
STPackageInitialized = PETSC_TRUE;
/* Register Classes */
ierr = PetscClassIdRegister("Spectral Transform",&ST_CLASSID);CHKERRQ(ierr);
/* Register Constructors */
ierr = STRegisterAll(path);CHKERRQ(ierr);
/* Register Events */
ierr = PetscLogEventRegister("STSetUp",ST_CLASSID,&ST_SetUp);CHKERRQ(ierr);
ierr = PetscLogEventRegister("STApply",ST_CLASSID,&ST_Apply);CHKERRQ(ierr);
ierr = PetscLogEventRegister("STApplyTranspose",ST_CLASSID,&ST_ApplyTranspose); CHKERRQ(ierr);
/* Process info exclusions */
ierr = PetscOptionsGetString(PETSC_NULL, "-info_exclude", logList, 256, &opt);CHKERRQ(ierr);
if (opt) {
ierr = PetscStrstr(logList, "st", &className);CHKERRQ(ierr);
if (className) {
ierr = PetscInfoDeactivateClass(ST_CLASSID);CHKERRQ(ierr);
}
}
/* Process summary exclusions */
ierr = PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);CHKERRQ(ierr);
if (opt) {
ierr = PetscStrstr(logList, "st", &className);CHKERRQ(ierr);
if (className) {
ierr = PetscLogEventDeactivateClass(ST_CLASSID);CHKERRQ(ierr);
}
}
ierr = PetscRegisterFinalize(STFinalizePackage);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STDestroy"
/*@
STDestroy - Destroys ST context that was created with STCreate().
Collective on ST
Input Parameter:
. st - the spectral transformation context
Level: beginner
.seealso: STCreate(), STSetUp()
@*/
PetscErrorCode STDestroy(ST st)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (--((PetscObject)st)->refct > 0) PetscFunctionReturn(0);
/* if memory was published with AMS then destroy it */
ierr = PetscObjectDepublish(st);CHKERRQ(ierr);
if (st->ops->destroy) { ierr = (*st->ops->destroy)(st);CHKERRQ(ierr); }
if (st->A) { ierr = MatDestroy(st->A);CHKERRQ(ierr); }
if (st->B) { ierr = MatDestroy(st->B);CHKERRQ(ierr); }
if (st->ksp) { ierr = KSPDestroy(st->ksp);CHKERRQ(ierr); }
if (st->w) { ierr = VecDestroy(st->w);CHKERRQ(ierr); }
if (st->D) { ierr = VecDestroy(st->D);CHKERRQ(ierr); }
if (st->wb){ ierr = VecDestroy(st->wb);CHKERRQ(ierr); }
if (st->shift_matrix != ST_MATMODE_INPLACE && st->mat) {
ierr = MatDestroy(st->mat);CHKERRQ(ierr);
}
ierr = PetscHeaderDestroy(st);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STCreate"
/*@C
STCreate - Creates a spectral transformation context.
Collective on MPI_Comm
Input Parameter:
. comm - MPI communicator
Output Parameter:
. st - location to put the spectral transformation context
Level: beginner
.seealso: STSetUp(), STApply(), STDestroy(), ST
@*/
PetscErrorCode STCreate(MPI_Comm comm,ST *newst)
{
PetscErrorCode ierr;
ST st;
const char *prefix;
PetscFunctionBegin;
PetscValidPointer(newst,2);
*newst = 0;
ierr = PetscHeaderCreate(st,_p_ST,struct _STOps,ST_CLASSID,-1,"ST",comm,STDestroy,STView);CHKERRQ(ierr);
ierr = PetscMemzero(st->ops,sizeof(struct _STOps));CHKERRQ(ierr);
st->A = 0;
st->B = 0;
st->sigma = 0.0;
st->sigma_set = PETSC_FALSE;
st->defsigma = 0.0;
st->data = 0;
st->setupcalled = 0;
st->w = 0;
st->D = 0;
st->wb = 0;
st->shift_matrix = ST_MATMODE_COPY;
st->str = DIFFERENT_NONZERO_PATTERN;
ierr = KSPCreate(((PetscObject)st)->comm,&st->ksp);CHKERRQ(ierr);
ierr = STGetOptionsPrefix(st,&prefix);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(st->ksp,prefix);CHKERRQ(ierr);
ierr = KSPAppendOptionsPrefix(st->ksp,"st_");CHKERRQ(ierr);
ierr = PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1);CHKERRQ(ierr);
*newst = st;
ierr = PetscPublishAll(st);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetOperators"
/*@
STSetOperators - Sets the matrices associated with the eigenvalue problem.
Collective on ST and Mat
Input Parameters:
+ st - the spectral transformation context
. A - the matrix associated with the eigensystem
- B - the second matrix in the case of generalized eigenproblems
Notes:
To specify a standard eigenproblem, use PETSC_NULL for B.
Level: intermediate
.seealso: STGetOperators()
@*/
PetscErrorCode STSetOperators(ST st,Mat A,Mat B)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidHeaderSpecific(A,MAT_CLASSID,2);
if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
PetscCheckSameComm(st,1,A,2);
if (B) PetscCheckSameComm(st,1,B,3);
ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
if (st->A) { ierr = MatDestroy(st->A);CHKERRQ(ierr); }
st->A = A;
if (B) { ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); }
if (st->B) { ierr = MatDestroy(st->B);CHKERRQ(ierr); }
st->B = B;
st->setupcalled = 0;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetOperators"
/*@C
STGetOperators - Gets the matrices associated with the eigensystem.
Not collective, though parallel Mats are returned if the ST is parallel
Input Parameter:
. st - the spectral transformation context
Output Parameters:
. A - the matrix associated with the eigensystem
- B - the second matrix in the case of generalized eigenproblems
Level: intermediate
.seealso: STSetOperators()
@*/
PetscErrorCode STGetOperators(ST st,Mat *A,Mat *B)
{
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (A) *A = st->A;
if (B) *B = st->B;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetShift"
/*@
STSetShift - Sets the shift associated with the spectral transformation.
Collective on ST
Input Parameters:
+ st - the spectral transformation context
- shift - the value of the shift
Note:
In some spectral transformations, changing the shift may have associated
a lot of work, for example recomputing a factorization.
Level: beginner
@*/
PetscErrorCode STSetShift(ST st,PetscScalar shift)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (st->sigma != shift) {
if (st->ops->setshift) {
ierr = (*st->ops->setshift)(st,shift); CHKERRQ(ierr);
}
}
st->sigma = shift;
st->sigma_set = PETSC_TRUE;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetShift"
/*@
STGetShift - Gets the shift associated with the spectral transformation.
Not collective
Input Parameter:
. st - the spectral transformation context
Output Parameter:
. shift - the value of the shift
Level: beginner
@*/
PetscErrorCode STGetShift(ST st,PetscScalar* shift)
{
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (shift) *shift = st->sigma;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetDefaultShift"
/*@
STSetDefaultShift - Sets the value of the shift that should be employed if
the user did not specify one.
Collective on ST
Input Parameters:
+ st - the spectral transformation context
- defaultshift - the default value of the shift
Level: developer
@*/
PetscErrorCode STSetDefaultShift(ST st,PetscScalar defaultshift)
{
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
st->defsigma = defaultshift;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetBalanceMatrix"
/*@
STSetBalanceMatrix - Sets the diagonal matrix to be used for balancing.
Collective on ST and Vec
Input Parameters:
+ st - the spectral transformation context
- D - the diagonal matrix (represented as a vector)
Notes:
If this matrix is set, STApply will effectively apply D*OP*D^{-1}.
Balancing is usually set via EPSSetBalance, but the advanced user may use
this function to bypass the usual balancing methods.
Level: developer
.seealso: EPSSetBalance(), STApply(), STGetBalanceMatrix()
@*/
PetscErrorCode STSetBalanceMatrix(ST st,Vec D)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidHeaderSpecific(D,VEC_CLASSID,2);
PetscCheckSameComm(st,1,D,2);
ierr = PetscObjectReference((PetscObject)D);CHKERRQ(ierr);
if (st->D) {
ierr = VecDestroy(st->D); CHKERRQ(ierr);
}
st->D = D;
if (!st->wb) {
ierr = VecDuplicate(st->D,&st->wb); CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetBalanceMatrix"
/*@
STGetBalanceMatrix - Gets the balance matrix used by the spectral transformation.
Not collective, but vector is shared by all processors that share the ST
Input Parameter:
. st - the spectral transformation context
Output Parameter:
. D - the diagonal matrix (represented as a vector)
Note:
If the matrix was not set, a null pointer will be returned.
Level: developer
.seealso: STSetBalanceMatrix()
@*/
PetscErrorCode STGetBalanceMatrix(ST st,Vec *D)
{
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
PetscValidPointer(D,2);
*D = st->D;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STSetOptionsPrefix"
/*@C
STSetOptionsPrefix - Sets the prefix used for searching for all
ST options in the database.
Collective on ST
Input Parameters:
+ st - the spectral transformation context
- prefix - the prefix string to prepend to all ST option requests
Notes:
A hyphen (-) must NOT be given at the beginning of the prefix name.
The first character of all runtime options is AUTOMATICALLY the
hyphen.
Level: advanced
.seealso: STAppendOptionsPrefix(), STGetOptionsPrefix()
@*/
PetscErrorCode STSetOptionsPrefix(ST st,const char *prefix)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
ierr = PetscObjectSetOptionsPrefix((PetscObject)st,prefix);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(st->ksp,prefix);CHKERRQ(ierr);
ierr = KSPAppendOptionsPrefix(st->ksp,"st_");CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STAppendOptionsPrefix"
/*@C
STAppendOptionsPrefix - Appends to the prefix used for searching for all
ST options in the database.
Collective on ST
Input Parameters:
+ st - the spectral transformation context
- prefix - the prefix string to prepend to all ST option requests
Notes:
A hyphen (-) must NOT be given at the beginning of the prefix name.
The first character of all runtime options is AUTOMATICALLY the
hyphen.
Level: advanced
.seealso: STSetOptionsPrefix(), STGetOptionsPrefix()
@*/
PetscErrorCode STAppendOptionsPrefix(ST st,const char *prefix)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
ierr = PetscObjectAppendOptionsPrefix((PetscObject)st,prefix);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);CHKERRQ(ierr);
ierr = KSPAppendOptionsPrefix(st->ksp,"st_");CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STGetOptionsPrefix"
/*@C
STGetOptionsPrefix - Gets the prefix used for searching for all
ST options in the database.
Not Collective
Input Parameters:
. st - the spectral transformation context
Output Parameters:
. prefix - pointer to the prefix string used, is returned
Notes: On the Fortran side, the user should pass in a string 'prefix' of
sufficient length to hold the prefix.
Level: advanced
.seealso: STSetOptionsPrefix(), STAppendOptionsPrefix()
@*/
PetscErrorCode STGetOptionsPrefix(ST st,const char *prefix[])
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
ierr = PetscObjectGetOptionsPrefix((PetscObject)st, prefix);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STView"
/*@C
STView - Prints the ST data structure.
Collective on ST
Input Parameters:
+ ST - the ST context
- viewer - optional visualization context
Note:
The available visualization contexts include
+ PETSC_VIEWER_STDOUT_SELF - standard output (default)
- PETSC_VIEWER_STDOUT_WORLD - synchronized standard
output where only the first processor opens
the file. All other processors send their
data to the first processor to print.
The user can open an alternative visualization contexts with
PetscViewerASCIIOpen() (output to a specified file).
Level: beginner
.seealso: EPSView(), PetscViewerASCIIOpen()
@*/
PetscErrorCode STView(ST st,PetscViewer viewer)
{
PetscErrorCode ierr;
const STType cstr;
const char* str;
PetscTruth isascii,isstring;
PetscViewerFormat format;
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_CLASSID,1);
if (!viewer) viewer = PETSC_VIEWER_STDOUT_(((PetscObject)st)->comm);
PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
PetscCheckSameComm(st,1,viewer,2);
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
if (isascii) {
ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"ST Object:\n");CHKERRQ(ierr);
ierr = STGetType(st,&cstr);CHKERRQ(ierr);
if (cstr) {
ierr = PetscViewerASCIIPrintf(viewer," type: %s\n",cstr);CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," type: not yet set\n");CHKERRQ(ierr);
}
#if !defined(PETSC_USE_COMPLEX)
ierr = PetscViewerASCIIPrintf(viewer," shift: %g\n",st->sigma);CHKERRQ(ierr);
#else
ierr = PetscViewerASCIIPrintf(viewer," shift: %g+%g i\n",PetscRealPart(st->sigma),PetscImaginaryPart(st->sigma));CHKERRQ(ierr);
#endif
switch (st->shift_matrix) {
case ST_MATMODE_COPY:
break;
case ST_MATMODE_INPLACE:
ierr = PetscViewerASCIIPrintf(viewer,"Shifting the matrix and unshifting at exit\n");CHKERRQ(ierr);
break;
case ST_MATMODE_SHELL:
ierr = PetscViewerASCIIPrintf(viewer,"Using a shell matrix\n");CHKERRQ(ierr);
break;
}
if (st->B && st->shift_matrix != ST_MATMODE_SHELL) {
switch (st->str) {
case SAME_NONZERO_PATTERN: str = "same nonzero pattern";break;
case DIFFERENT_NONZERO_PATTERN: str = "different nonzero pattern";break;
case SUBSET_NONZERO_PATTERN: str = "subset nonzero pattern";break;
default: SETERRQ(1,"Wrong structure flag");
}
ierr = PetscViewerASCIIPrintf(viewer,"Matrices A and B have %s\n",str);CHKERRQ(ierr);
}
if (st->ops->view) {
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = (*st->ops->view)(st,viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
}
} else if (isstring) {
ierr = STGetType(st,&cstr);CHKERRQ(ierr);
ierr = PetscViewerStringSPrintf(viewer," %-7.7s",cstr);CHKERRQ(ierr);
if (st->ops->view) {ierr = (*st->ops->view)(st,viewer);CHKERRQ(ierr);}
} else {
SETERRQ1(1,"Viewer type %s not supported by ST",((PetscObject)viewer)->type_name);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STView_Default"
PetscErrorCode STView_Default(ST st,PetscViewer viewer)
{
PetscErrorCode ierr;
PetscTruth isascii,isstring;
PetscFunctionBegin;
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);CHKERRQ(ierr);
if (isascii) {
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"Associated KSP object\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"------------------------------\n");CHKERRQ(ierr);
ierr = KSPView(st->ksp,viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"------------------------------\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
} else if (isstring) {
ierr = KSPView(st->ksp,viewer);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
/*MC
STRegisterDynamic - Adds a method to the spectral transformation package.
Synopsis:
STRegisterDynamic(char *name_solver,char *path,char *name_create,PetscErrorCode (*routine_create)(ST))
Not collective
Input Parameters:
+ name_solver - name of a new user-defined solver
. path - path (either absolute or relative) the library containing this solver
. name_create - name of routine to create method context
- routine_create - routine to create method context
Notes:
STRegisterDynamic() may be called multiple times to add several user-defined spectral transformations.
If dynamic libraries are used, then the fourth input argument (routine_create)
is ignored.
Sample usage:
.vb
STRegisterDynamic("my_solver","/home/username/my_lib/lib/libO/solaris/mylib.a",
"MySolverCreate",MySolverCreate);
.ve
Then, your solver can be chosen with the procedural interface via
$ STSetType(st,"my_solver")
or at runtime via the option
$ -st_type my_solver
Level: advanced
$PETSC_DIR, $PETSC_ARCH and $PETSC_LIB_DIR occuring in pathname will be replaced with appropriate values.
.seealso: STRegisterDestroy(), STRegisterAll()
M*/
#undef __FUNCT__
#define __FUNCT__ "STRegister"
/*@C
STRegister - See STRegisterDynamic()
Level: advanced
@*/
PetscErrorCode STRegister(const char *sname,const char *path,const char *name,PetscErrorCode (*function)(ST))
{
PetscErrorCode ierr;
char fullname[256];
PetscFunctionBegin;
ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr);
ierr = PetscFListAdd(&STList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "STRegisterDestroy"
/*@
STRegisterDestroy - Frees the list of ST methods that were
registered by STRegisterDynamic().
Not Collective
Level: advanced
.seealso: STRegisterDynamic(), STRegisterAll()
@*/
PetscErrorCode STRegisterDestroy(void)
{
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscFListDestroy(&STList);CHKERRQ(ierr);
ierr = STRegisterAll(PETSC_NULL);CHKERRQ(ierr);
PetscFunctionReturn(0);
}