Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2 → Rev 6

/trunk/src/st/impls/sinvert/sinvert.h
0,0 → 1,16
/*
Implements the shift-and-invert technique for eigenvalue problems.
*/
 
#if !defined(__SINVERT_H)
#define __SINVERT_H
 
typedef struct {
Mat A, B;
Vec w;
PetscScalar sigma;
} CTX_SINV;
 
extern int MatCreateMatSinvert(ST,Mat*);
 
#endif
/trunk/src/st/impls/sinvert/shellmat.c
0,0 → 1,101
/*
This file contains the subroutines which implement various operations
of the matrix associated to the shift-and-invert technique for eigenvalue
problems, and also a subroutine to create it.
*/
 
#include "src/st/stimpl.h"
#include "sinvert.h"
 
#undef __FUNCT__
#define __FUNCT__ "MatSinvert_Mult"
static int MatSinvert_Mult(Mat A,Vec x,Vec y)
{
int ierr;
CTX_SINV *ctx;
PetscScalar alpha;
 
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
alpha = -ctx->sigma;
 
if (ctx->B) { /* y = (A - sB) x */
ierr = MatMult(ctx->B,x,ctx->w);CHKERRQ(ierr);
ierr = MatMult(ctx->A,x,y);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,ctx->w,y);CHKERRQ(ierr);
}
else { /* y = (A - sI) x */
ierr = MatMult(ctx->A,x,y);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,x,y);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatSinvert_GetDiagonal"
static int MatSinvert_GetDiagonal(Mat A,Vec diag)
{
int ierr;
CTX_SINV *ctx;
PetscScalar alpha;
Vec diagb;
 
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
alpha = -ctx->sigma;
 
ierr = MatGetDiagonal(ctx->A,diag);CHKERRQ(ierr);
if (ctx->B) {
ierr = VecDuplicate(diag,&diagb);CHKERRQ(ierr);
ierr = MatGetDiagonal(ctx->B,diagb);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,diagb,diag);CHKERRQ(ierr);
ierr = VecDestroy(diagb);CHKERRQ(ierr);
}
else {
ierr = VecShift(&alpha,diag);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatSinvert_Destroy"
static int MatSinvert_Destroy(Mat A)
{
CTX_SINV *ctx;
int ierr;
 
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
if (ctx->B) { ierr = VecDestroy(ctx->w);CHKERRQ(ierr); }
ierr = PetscFree(ctx);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "MatCreateMatSinvert"
int MatCreateMatSinvert(ST st,Mat *mat)
{
int n, m, N, M, ierr;
PetscTruth hasA, hasB;
CTX_SINV *ctx;
 
PetscFunctionBegin;
ierr = PetscNew(CTX_SINV,&ctx);CHKERRQ(ierr);
PetscMemzero(ctx,sizeof(CTX_SINV));
PetscLogObjectMemory(st,sizeof(CTX_SINV));
ctx->A = st->A;
ctx->B = st->B;
ctx->sigma = st->sigma;
if (st->B) { ierr = VecDuplicate(st->vec,&ctx->w);CHKERRQ(ierr); }
ierr = MatGetSize(st->A,&M,&N);CHKERRQ(ierr);
ierr = MatGetLocalSize(st->A,&m,&n);CHKERRQ(ierr);
ierr = MatCreateShell(st->comm,m,n,M,N,(void*)ctx,mat);CHKERRQ(ierr);
ierr = MatShellSetOperation(*mat,MATOP_MULT,(void(*)(void))MatSinvert_Mult);CHKERRQ(ierr);
ierr = MatShellSetOperation(*mat,MATOP_DESTROY,(void(*)(void))MatSinvert_Destroy);CHKERRQ(ierr);
 
ierr = MatHasOperation(st->A,MATOP_GET_DIAGONAL,&hasA);CHKERRQ(ierr);
if (st->B) { ierr = MatHasOperation(st->B,MATOP_GET_DIAGONAL,&hasB);CHKERRQ(ierr); }
if ( (hasA && !st->B) || (hasA && hasB) ) {
ierr = MatShellSetOperation(*mat,MATOP_GET_DIAGONAL,(void(*)(void))MatSinvert_GetDiagonal);CHKERRQ(ierr);
}
 
PetscFunctionReturn(0);
}
 
/trunk/src/st/impls/sinvert/sinvert.c
0,0 → 1,354
/*
Implements the shift-and-invert technique for eigenvalue problems.
*/
#include "src/st/stimpl.h" /*I "slepcst.h" I*/
#include "sinvert.h"
 
typedef struct {
PetscTruth shift_matrix; /* shift matrix rather than use shell mat */
MatStructure str; /* whether matrices have the same pattern or not */
Mat mat;
Vec w;
} ST_SINV;
 
#undef __FUNCT__
#define __FUNCT__ "STApply_Sinvert"
static int STApply_Sinvert(ST st,Vec x,Vec y)
{
int ierr;
ST_SINV *ctx = (ST_SINV *) st->data;
 
PetscFunctionBegin;
if (st->B) {
/* generalized eigenproblem: y = (A - sB)^-1 B x */
ierr = MatMult(st->B,x,ctx->w);CHKERRQ(ierr);
ierr = STAssociatedSLESSolve(st,ctx->w,y);CHKERRQ(ierr);
}
else {
/* standard eigenproblem: y = (A - sI)^-1 x */
ierr = STAssociatedSLESSolve(st,x,y);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STApplyNoB_Sinvert"
static int STApplyNoB_Sinvert(ST st,Vec x,Vec y)
{
int ierr;
 
PetscFunctionBegin;
ierr = STAssociatedSLESSolve(st,x,y);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STBackTransform_Sinvert"
int STBackTransform_Sinvert(ST st,PetscScalar *eigr,PetscScalar *eigi)
{
PetscFunctionBegin;
/* Note that this is not correct in the case of the RQI solver */
if (eigr) *eigr = 1.0 / *eigr + st->sigma;
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STPost_Sinvert"
int STPost_Sinvert(ST st)
{
ST_SINV *ctx = (ST_SINV *) st->data;
PetscScalar alpha;
int ierr;
 
PetscFunctionBegin;
if( ctx->shift_matrix ) {
alpha = st->sigma;
if( st->B ) { ierr = MatAXPY(&alpha,st->B,st->A,ctx->str);CHKERRQ(ierr); }
else { ierr = MatShift( &alpha, st->A ); CHKERRQ(ierr); }
st->setupcalled = 0;
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STSetUp_Sinvert"
static int STSetUp_Sinvert(ST st)
{
int ierr;
ST_SINV *ctx = (ST_SINV *) st->data;
PetscScalar alpha;
 
PetscFunctionBegin;
 
if (ctx->shift_matrix) {
alpha = -st->sigma;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,ctx->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
ierr = SLESSetOperators(st->sles,st->A,st->A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
}
else {
ierr = MatCreateMatSinvert(st,&ctx->mat);CHKERRQ(ierr);
ierr = SLESSetOperators(st->sles,ctx->mat,ctx->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
}
if (st->B && !ctx->w) { ierr = VecDuplicate(st->vec,&ctx->w);CHKERRQ(ierr); }
ierr = SLESSetUp(st->sles,st->vec,st->vec);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STSetShift_Sinvert"
static int STSetShift_Sinvert(ST st,PetscScalar newshift)
{
int ierr;
ST_SINV *stctx = (ST_SINV *) st->data;
PetscScalar alpha;
CTX_SINV *ctx;
 
PetscFunctionBegin;
 
/* Nothing to be done if STSetUp has not been called yet */
if (!st->setupcalled) PetscFunctionReturn(0);
 
if (stctx->shift_matrix) {
/* Undo previous operations */
alpha = st->sigma;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,stctx->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
/* Apply new shift */
alpha = -newshift;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,stctx->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
ierr = SLESSetOperators(st->sles,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
}
else {
ierr = MatShellGetContext(stctx->mat,(void**)&ctx);CHKERRQ(ierr);
ctx->sigma = newshift;
ierr = SLESSetOperators(st->sles,stctx->mat,stctx->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
}
ierr = SLESSetUp(st->sles,st->vec,st->vec);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STDestroy_Sinvert"
static int STDestroy_Sinvert(ST st)
{
ST_SINV *ctx = (ST_SINV *) st->data;
int ierr;
 
PetscFunctionBegin;
if (!ctx->shift_matrix) { ierr = MatDestroy(ctx->mat);CHKERRQ(ierr); }
if (st->B) { ierr = VecDestroy(ctx->w);CHKERRQ(ierr); }
ierr = PetscFree(ctx);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STView_Sinvert"
static int STView_Sinvert(ST st,PetscViewer viewer)
{
ST_SINV *ctx = (ST_SINV *) st->data;
int ierr;
PetscTruth isascii;
char *str;
 
PetscFunctionBegin;
ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
if (!isascii) {
SETERRQ1(1,"Viewer type %s not supported for STSINV",((PetscObject)viewer)->type_name);
}
if (ctx->shift_matrix) {
ierr = PetscViewerASCIIPrintf(viewer,"Shifting the matrix and unshifting at exit\n");CHKERRQ(ierr);
if (st->B) {
switch (ctx->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;
}
ierr = PetscViewerASCIIPrintf(viewer,"Matrices A and B have %s\n",str);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "STSetFromOptions_Sinvert"
static int STSetFromOptions_Sinvert(ST st)
{
int ierr;
PetscTruth flg;
PC pc;
 
PetscFunctionBegin;
ierr = PetscOptionsHead("ST Shift-and-invert Options");CHKERRQ(ierr);
ierr = PetscOptionsName("-st_sinvert_shift_mat","Shift matrix explicitly","STSinvertSetShiftMat",&flg);CHKERRQ(ierr);
if (flg) {
ierr = STSinvertSetShiftMat(st);CHKERRQ(ierr);
}
else {
/* if shift_mat is set then the default preconditioner is ILU,
otherwise set Jacobi as the default */
ierr = SLESGetPC(st->sles,&pc); CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
}
ierr = PetscOptionsLogicalGroupBegin("-st_sinvert_same_pattern","same nonzero pattern","STSinvertSetMatStructure",&flg);CHKERRQ(ierr);
if (flg) {ierr = STSinvertSetMatStructure(st,SAME_NONZERO_PATTERN);CHKERRQ(ierr);}
ierr = PetscOptionsLogicalGroup("-st_sinvert_different_pattern","different nonzero pattern","STSinvertSetMatStructure",&flg);CHKERRQ(ierr);
if (flg) {ierr = STSinvertSetMatStructure(st,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);}
ierr = PetscOptionsLogicalGroupEnd("-st_sinvert_subset_pattern","subset nonzero pattern","STSinvertSetMatStructure",&flg);CHKERRQ(ierr);
if (flg) {ierr = STSinvertSetMatStructure(st,SUBSET_NONZERO_PATTERN);CHKERRQ(ierr);}
ierr = PetscOptionsTail();CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
/* -------------------------------------------------------------------------*/
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "STSinvertSetShiftMat_Sinvert"
int STSinvertSetShiftMat_Sinvert(ST st)
{
ST_SINV *ctx = (ST_SINV *) st->data;
 
PetscFunctionBegin;
ctx->shift_matrix = PETSC_TRUE;
PetscFunctionReturn(0);
}
EXTERN_C_END
 
#undef __FUNCT__
#define __FUNCT__ "STSinvertSetShiftMat"
/*@
STSinvertSetShiftMat - Sets a flag to indicate that the matrix is
being shifted at STSetUp() and unshifted at the end of the computations.
 
Collective on ST
 
Input Parameters:
. st - the spectral transformation context
 
Options Database Key:
. -st_sinvert_shift_mat - Activates STSinvertSetShiftMat()
 
Note:
By default, the matrix is not shifted explicitly. Instead, the solver
works with an implicit shell matrix that represents the shifted matrix,
in which case only the Jacobi preconditioning is available for the linear
solves performed in each iteration of the eigensolver.
Level: intermediate
 
.seealso: STSetOperators()
@*/
int STSinvertSetShiftMat(ST st)
{
int ierr, (*f)(ST);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_COOKIE);
ierr = PetscObjectQueryFunction((PetscObject)st,"STSinvertSetShiftMat_C",(void (**)(void))&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(st);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "STSinvertSetMatStructure_Sinvert"
int STSinvertSetMatStructure_Sinvert(ST st,MatStructure str)
{
ST_SINV *ctx = (ST_SINV *) st->data;
 
PetscFunctionBegin;
ctx->str = str;
PetscFunctionReturn(0);
}
EXTERN_C_END
 
#undef __FUNCT__
#define __FUNCT__ "STSinvertSetMatStructure"
/*@
STSinvertSetMatStructure - Sets an internal MatStructure attribute to
indicate which is the relation of the sparsity pattern of the two matrices
A and B constituting the generalized eigenvalue problem. This function
has no effect in the case of standard eigenproblems.
 
Collective on ST
 
Input Parameters:
+ st - the spectral transformation context
- str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or
SUBSET_NONZERO_PATTERN
 
Options Database Key:
+ -st_sinvert_same_pattern - Indicates A and B have the same nonzero pattern
. -st_sinvert_different_pattern - Indicates A and B have different nonzero pattern
- -st_sinvert_subset_pattern - Indicates B's nonzero pattern is a subset of B's
 
Note:
By default, the sparsity patterns are assumed to be different. If the
patterns are equal or a subset then it is recommended to set this attribute
for efficiency reasons (in particular, for internal MatAXPY operations).
Level: advanced
 
.seealso: STSetOperators()
@*/
int STSinvertSetMatStructure(ST st,MatStructure str)
{
int ierr, (*f)(ST,MatStructure);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_COOKIE);
ierr = PetscObjectQueryFunction((PetscObject)st,"STSinvertSetMatStructure_C",(void (**)(void))&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(st,str);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
/* ---------------------------------------------------------------------------*/
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "STCreate_Sinvert"
int STCreate_Sinvert(ST st)
{
int ierr;
char *prefix;
ST_SINV *ctx;
 
PetscFunctionBegin;
ierr = PetscNew(ST_SINV,&ctx); CHKERRQ(ierr);
PetscMemzero(ctx,sizeof(ST_SINV));
PetscLogObjectMemory(st,sizeof(ST_SINV));
st->numberofshifts = 1;
st->data = (void *) ctx;
 
st->ops->apply = STApply_Sinvert;
st->ops->applynoB = STApplyNoB_Sinvert;
st->ops->postsolve = STPost_Sinvert;
st->ops->backtr = STBackTransform_Sinvert;
st->ops->setup = STSetUp_Sinvert;
st->ops->setshift = STSetShift_Sinvert;
st->ops->destroy = STDestroy_Sinvert;
st->ops->setfromoptions = STSetFromOptions_Sinvert;
st->ops->view = STView_Sinvert;
 
ierr = SLESCreate(st->comm,&st->sles);CHKERRQ(ierr);
ierr = STGetOptionsPrefix(st,&prefix);CHKERRQ(ierr);
ierr = SLESSetOptionsPrefix(st->sles,prefix);CHKERRQ(ierr);
ierr = SLESAppendOptionsPrefix(st->sles,"st_");CHKERRQ(ierr);
ctx->shift_matrix = PETSC_FALSE;
ctx->str = DIFFERENT_NONZERO_PATTERN;
 
ierr = PetscObjectComposeFunctionDynamic((PetscObject)st,"STSinvertSetShiftMat_C","STSinvertSetShiftMat_Sinvert",
STSinvertSetShiftMat_Sinvert);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)st,"STSinvertSetMatStructure_C","STSinvertSetMatStructure_Sinvert",
STSinvertSetMatStructure_Sinvert);CHKERRQ(ierr);
 
PetscFunctionReturn(0);
}
EXTERN_C_END
 
/trunk/src/st/impls/sinvert/makefile
0,0 → 1,17
 
ALL: lib
 
CFLAGS =
FFLAGS =
SOURCEC = sinvert.c shellmat.c
SOURCEF =
SOURCEH = sinvert.h
OBJSC = sinvert.o shellmat.o
LIBBASE = libslepc
DIRS =
MANSEC = ST
LOCDIR = src/st/impls/sinvert/
 
include ${SLEPC_DIR}/bmake/slepc_common