Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2806 → Rev 2807

/trunk/include/slepcps.h
89,7 → 89,7
 
Level: advanced
 
.seealso: PSAllocate(), PSGetArray(), PSGetArrayReal(), PSComputeVector()
.seealso: PSAllocate(), PSGetArray(), PSGetArrayReal(), PSVectors()
E*/
typedef enum { PS_MAT_A,
PS_MAT_B,
128,7 → 128,7
extern PetscErrorCode PSRestoreArray(PS,PSMatType,PetscScalar *a[]);
extern PetscErrorCode PSGetArrayReal(PS,PSMatType,PetscReal *a[]);
extern PetscErrorCode PSRestoreArrayReal(PS,PSMatType,PetscReal *a[]);
extern PetscErrorCode PSComputeVector(PS,PetscInt,PSMatType,PetscBool*);
extern PetscErrorCode PSVectors(PS,PSMatType,PetscInt*,PetscReal*);
extern PetscErrorCode PSSolve(PS,PetscScalar*,PetscScalar*);
extern PetscErrorCode PSSort(PS,PetscScalar*,PetscScalar*,PetscErrorCode (*)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void*);
extern PetscErrorCode PSCond(PS,PetscReal*);
/trunk/include/slepc-private/psimpl.h
24,7 → 24,7
 
#include <slepcps.h>
 
extern PetscLogEvent PS_Solve,PS_Sort,PS_Other;
extern PetscLogEvent PS_Solve,PS_Sort,PS_Vectors,PS_Other;
extern const char *PSMatName[];
 
typedef struct _PSOps *PSOps;
32,7 → 32,7
struct _PSOps {
PetscErrorCode (*allocate)(PS,PetscInt);
PetscErrorCode (*view)(PS,PetscViewer);
PetscErrorCode (*computevector)(PS,PetscInt,PSMatType,PetscBool*);
PetscErrorCode (*vectors)(PS,PSMatType,PetscInt*,PetscReal*);
PetscErrorCode (*solve)(PS,PetscScalar*,PetscScalar*);
PetscErrorCode (*sort)(PS,PetscScalar*,PetscScalar*,PetscErrorCode(*)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void*);
PetscErrorCode (*cond)(PS,PetscReal*);
/trunk/src/ps/pshep.c
92,6 → 92,36
}
 
#undef __FUNCT__
#define __FUNCT__ "PSVectors_HEP"
PetscErrorCode PSVectors_HEP(PS ps,PSMatType mat,PetscInt *k,PetscReal *rnorm)
{
PetscScalar *Q = ps->mat[PS_MAT_Q];
PetscInt ld = ps->ld;
PetscErrorCode ierr;
 
PetscFunctionBegin;
if (ps->state<PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");
switch (mat) {
case PS_MAT_X:
case PS_MAT_Y:
if (k) {
ierr = PetscMemcpy(ps->mat[mat]+(*k)*ld,Q+(*k)*ld,ld*sizeof(PetscScalar));CHKERRQ(ierr);
} else {
ierr = PetscMemcpy(ps->mat[mat],Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
}
if (rnorm) *rnorm = PetscAbsScalar(Q[ps->n-1+(*k)*ld]);
break;
case PS_MAT_U:
case PS_MAT_VT:
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented yet");
break;
default:
SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PSSolve_HEP"
PetscErrorCode PSSolve_HEP(PS ps,PetscScalar *wr,PetscScalar *wi)
{
284,7 → 314,7
ps->nmeth = 1;
ps->ops->allocate = PSAllocate_HEP;
ps->ops->view = PSView_HEP;
//ps->ops->computevector = PSComputeVector_HEP;
ps->ops->vectors = PSVectors_HEP;
ps->ops->solve = PSSolve_HEP;
ps->ops->sort = PSSort_HEP;
ps->ops->cond = PSCond_HEP;
/trunk/src/ps/psbasic.c
26,7 → 26,7
PetscFList PSList = 0;
PetscBool PSRegisterAllCalled = PETSC_FALSE;
PetscClassId PS_CLASSID = 0;
PetscLogEvent PS_Solve = 0,PS_Sort = 0,PS_Other = 0;
PetscLogEvent PS_Solve = 0,PS_Sort = 0,PS_Vectors = 0,PS_Other = 0;
static PetscBool PSPackageInitialized = PETSC_FALSE;
const char *PSMatName[PS_NUM_MAT] = {"A","B","C","T","Q","X","Y","U","VT","W"};
 
80,6 → 80,7
/* Register Events */
ierr = PetscLogEventRegister("PSSolve",PS_CLASSID,&PS_Solve);CHKERRQ(ierr);
ierr = PetscLogEventRegister("PSSort",PS_CLASSID,&PS_Sort);CHKERRQ(ierr);
ierr = PetscLogEventRegister("PSVectors",PS_CLASSID,&PS_Vectors);CHKERRQ(ierr);
ierr = PetscLogEventRegister("PSOther",PS_CLASSID,&PS_Other);CHKERRQ(ierr);
/* Process info exclusions */
ierr = PetscOptionsGetString(PETSC_NULL,"-info_exclude",logList,256,&opt);CHKERRQ(ierr);
521,7 → 522,7
/*@
PSAllocate - Allocates memory for internal storage or matrices in PS.
 
Collective on PS
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
/trunk/src/ps/psops.c
54,7 → 54,7
/*@
PSSetState - Change the state of the PS object.
 
Collective on PS
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
125,7 → 125,7
/*@
PSSetDimensions - Resize the matrices in the PS object.
 
Collective on PS
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
333,7 → 333,7
/*@
PSSolve - Solves the problem.
 
Not Collective
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
370,71 → 370,126
}
 
#undef __FUNCT__
#define __FUNCT__ "PSCond"
#define __FUNCT__ "PSSort"
/*@C
PSCond - Compute the inf-norm condition number of the first matrix
as cond(A) = norm(A)*norm(inv(A)).
PSSort - Reorders the condensed form computed by PSSolve() according to
a given sorting criterion.
 
Not Collective
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
- cond - the computed condition number
. eigr - array to store the sorted eigenvalues (real part)
- eigi - array to store the sorted eigenvalues (imaginary part)
 
Level: advanced
 
.seealso: PSSolve()
@*/
PetscErrorCode PSCond(PS ps,PetscReal *cond)
PetscErrorCode PSSort(PS ps,PetscScalar *eigr,PetscScalar *eigi,PetscErrorCode (*comp_func)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void *comp_ctx)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(ps,PS_CLASSID,1);
PetscValidPointer(cond,2);
if (!ps->ops->cond) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_SUP,"PS type %s",((PetscObject)ps)->type_name);
ierr = PetscLogEventBegin(PS_Other,ps,0,0,0);CHKERRQ(ierr);
PetscValidPointer(eigr,2);
if (ps->state!=PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");
if (!ps->ops->sort) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_SUP,"PS type %s",((PetscObject)ps)->type_name);
if (ps->state>=PS_STATE_SORTED) PetscFunctionReturn(0);
ierr = PetscLogEventBegin(PS_Sort,ps,0,0,0);CHKERRQ(ierr);
ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
ierr = (*ps->ops->cond)(ps,cond);CHKERRQ(ierr);
ierr = (*ps->ops->sort)(ps,eigr,eigi,comp_func,comp_ctx);CHKERRQ(ierr);
ierr = PetscFPTrapPop();CHKERRQ(ierr);
ierr = PetscLogEventEnd(PS_Other,ps,0,0,0);CHKERRQ(ierr);
ierr = PetscLogEventEnd(PS_Sort,ps,0,0,0);CHKERRQ(ierr);
ps->state = PS_STATE_SORTED;
ierr = PetscObjectStateIncrease((PetscObject)ps);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PSSort"
#define __FUNCT__ "PSVectors"
/*@
PSVectors - Compute vectors associated to the projected system such
as eigenvectors.
 
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
- mat - the matrix, used to indicate which vectors are required
 
Input/Output Parameter:
- k - (optional) index of vector to be computed
 
Output Parameter:
. rnorm - (optional) computed residual norm
 
Notes:
Allowed values for mat are PS_MAT_X, PS_MAT_Y, PS_MAT_U and PS_MAT_VT, to
compute right or left eigenvectors, or left or right singular vectors,
respectively.
 
If PETSC_NULL is passed in argument k then all vectors are computed,
otherwise k indicates which vector must be computed. In real non-symmetric
problems, on exit the index k will be incremented when a complex conjugate
pair is found.
 
This function can be invoked after the projected problem has been solved,
to get the residual norm estimate of the associated Ritz pair. In that
case, the relevant information is returned in rnorm.
 
Level: advanced
 
.seealso: PSSolve()
@*/
PetscErrorCode PSVectors(PS ps,PSMatType mat,PetscInt *k,PetscReal *rnorm)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(ps,PS_CLASSID,1);
PetscValidLogicalCollectiveEnum(ps,mat,2);
if (!ps->ld) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSAllocate() first");
if (!ps->ops->vectors) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_SUP,"PS type %s",((PetscObject)ps)->type_name);
ierr = PSAllocateMat_Private(ps,mat);CHKERRQ(ierr);
ierr = PetscLogEventBegin(PS_Vectors,ps,0,0,0);CHKERRQ(ierr);
ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
ierr = (*ps->ops->vectors)(ps,mat,k,rnorm);CHKERRQ(ierr);
ierr = PetscFPTrapPop();CHKERRQ(ierr);
ierr = PetscLogEventEnd(PS_Vectors,ps,0,0,0);CHKERRQ(ierr);
ierr = PetscObjectStateIncrease((PetscObject)ps);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PSCond"
/*@C
PSSort - Reorders the condensed form computed by PSSolve() according to
a given sorting criterion.
PSCond - Compute the inf-norm condition number of the first matrix
as cond(A) = norm(A)*norm(inv(A)).
 
Not Collective
 
Input Parameters:
+ ps - the projected system context
. eigr - array to store the sorted eigenvalues (real part)
- eigi - array to store the sorted eigenvalues (imaginary part)
- cond - the computed condition number
 
Level: advanced
 
.seealso: PSSolve()
@*/
PetscErrorCode PSSort(PS ps,PetscScalar *eigr,PetscScalar *eigi,PetscErrorCode (*comp_func)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void *comp_ctx)
PetscErrorCode PSCond(PS ps,PetscReal *cond)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(ps,PS_CLASSID,1);
PetscValidPointer(eigr,2);
if (ps->state!=PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");
if (!ps->ops->sort) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_SUP,"PS type %s",((PetscObject)ps)->type_name);
if (ps->state>=PS_STATE_SORTED) PetscFunctionReturn(0);
ierr = PetscLogEventBegin(PS_Sort,ps,0,0,0);CHKERRQ(ierr);
PetscValidPointer(cond,2);
if (!ps->ops->cond) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_SUP,"PS type %s",((PetscObject)ps)->type_name);
ierr = PetscLogEventBegin(PS_Other,ps,0,0,0);CHKERRQ(ierr);
ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
ierr = (*ps->ops->sort)(ps,eigr,eigi,comp_func,comp_ctx);CHKERRQ(ierr);
ierr = (*ps->ops->cond)(ps,cond);CHKERRQ(ierr);
ierr = PetscFPTrapPop();CHKERRQ(ierr);
ierr = PetscLogEventEnd(PS_Sort,ps,0,0,0);CHKERRQ(ierr);
ps->state = PS_STATE_SORTED;
ierr = PetscObjectStateIncrease((PetscObject)ps);CHKERRQ(ierr);
ierr = PetscLogEventEnd(PS_Other,ps,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
443,7 → 498,7
/*@C
PSTranslateHarmonic - Computes a translation of the projected matrix.
 
Not Collective
Logically Collective on PS
 
Input Parameters:
+ ps - the projected system context
/trunk/src/ps/examples/tests/makefile
27,7 → 27,7
EXAMPLESC = test1.c test2.c test3.c
EXAMPLESF =
MANSEC = PS
TESTS = test1 test2 test2
TESTS = test1 test2 test3
 
TESTEXAMPLES_C = test1.PETSc runtest1_1 test1.rm \
= test2.PETSc runtest2_1 test2.rm \
/trunk/src/ps/psnhep.c
49,6 → 49,155
}
 
#undef __FUNCT__
#define __FUNCT__ "PSVectors_NHEP_Eigen_Some"
PetscErrorCode PSVectors_NHEP_Eigen_Some(PS ps,PetscInt *k,PetscReal *rnorm,PetscBool left)
{
#if defined(SLEPC_MISSING_LAPACK_TREVC)
PetscFunctionBegin;
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");
#else
PetscErrorCode ierr;
PetscInt i;
PetscBLASInt mm,mout,info,ld,n,inc = 1;
PetscScalar tmp,done=1.0,zero=0.0;
PetscReal norm;
PetscBool iscomplex = PETSC_FALSE;
PetscBLASInt *select;
PetscScalar *A = ps->mat[PS_MAT_A];
PetscScalar *Q = ps->mat[PS_MAT_Q];
PetscScalar *X = ps->mat[PS_MAT_X];
PetscScalar *Y;
#if defined(PETSC_USE_COMPLEX)
PetscReal *rwork;
#endif
 
PetscFunctionBegin;
if (ps->state<PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");
if (left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for left eigenvectors");
n = PetscBLASIntCast(ps->n);
ld = PetscBLASIntCast(ps->ld);
if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
ierr = PSAllocateWork_Private(ps,0,0,ld);CHKERRQ(ierr);
select = ps->iwork;
for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;
 
/* Compute k-th eigenvector Y of A */
Y = X+(*k)*ld;
mm = iscomplex? 2: 1;
select[*k] = (PetscBLASInt)PETSC_TRUE;
#if !defined(PETSC_USE_COMPLEX)
if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;
ierr = PSAllocateWork_Private(ps,3*ld,0,0);CHKERRQ(ierr);
LAPACKtrevc_("R","S",select,&n,A,&ld,PETSC_NULL,&ld,Y,&ld,&mm,&mout,ps->work,&info);
#else
ierr = PSAllocateWork_Private(ps,2*ld,ld,0);CHKERRQ(ierr);
LAPACKtrevc_("R","S",select,&n,A,&ld,PETSC_NULL,&ld,Y,&ld,&mm,&mout,ps->work,ps->rwork,&info);
#endif
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTREVC %d",info);
if (mout != mm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Inconsistent arguments");
ierr = PetscMemcpy(ps->work,Y,mout*ld*sizeof(PetscScalar));CHKERRQ(ierr);
 
/* accumulate and normalize eigenvectors */
BLASgemv_("N",&n,&n,&done,Q,&ld,ps->work,&inc,&zero,Y,&inc);
#if !defined(PETSC_USE_COMPLEX)
if (iscomplex) BLASgemv_("N",&n,&n,&done,Q,&ld,ps->work+ld,&inc,&zero,Y+ld,&inc);
#endif
norm = BLASnrm2_(&n,Y,&inc);
#if !defined(PETSC_USE_COMPLEX)
tmp = BLASnrm2_(&n,Y+ld,&inc);
norm = SlepcAbsEigenvalue(norm,tmp);
#endif
tmp = 1.0 / norm;
BLASscal_(&n,&tmp,Y,&inc);
#if !defined(PETSC_USE_COMPLEX)
BLASscal_(&n,&tmp,Y+ld,&inc);
#endif
 
/* set output arguments */
if (iscomplex) (*k)++;
if (rnorm) {
if (iscomplex) *rnorm = SlepcAbsEigenvalue(Y[n-1],Y[n-1+ld]);
else *rnorm = PetscAbsScalar(Y[n-1]);
}
PetscFunctionReturn(0);
#endif
}
 
#undef __FUNCT__
#define __FUNCT__ "PSVectors_NHEP_Eigen_All"
PetscErrorCode PSVectors_NHEP_Eigen_All(PS ps,PetscBool left)
{
#if defined(SLEPC_MISSING_LAPACK_TREVC)
PetscFunctionBegin;
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");
#else
PetscErrorCode ierr;
PetscBLASInt n,ld,mout,info;
PetscScalar *A = ps->mat[PS_MAT_A];
PetscScalar *Q = ps->mat[PS_MAT_Q];
PetscScalar *X,*Y;
const char *side;
 
PetscFunctionBegin;
if (ps->state<PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");
n = PetscBLASIntCast(ps->n);
ld = PetscBLASIntCast(ps->ld);
if (left) {
X = PETSC_NULL;
Y = ps->mat[PS_MAT_Y];
side = "L";
ierr = PetscMemcpy(Y,Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
} else {
X = ps->mat[PS_MAT_X];
Y = PETSC_NULL;
side = "R";
ierr = PetscMemcpy(X,Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
}
#if !defined(PETSC_USE_COMPLEX)
ierr = PSAllocateWork_Private(ps,3*ld,0,0);CHKERRQ(ierr);
LAPACKtrevc_(side,"B",PETSC_NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ps->work,&info);
#else
ierr = PSAllocateWork_Private(ps,2*ld,ld,0);CHKERRQ(ierr);
LAPACKtrevc_(side,"B",PETSC_NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ps->work,ps->rwork,&info);
#endif
if (info) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
PetscFunctionReturn(0);
#endif
}
 
#undef __FUNCT__
#define __FUNCT__ "PSVectors_NHEP"
PetscErrorCode PSVectors_NHEP(PS ps,PSMatType mat,PetscInt *k,PetscReal *rnorm)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
switch (mat) {
case PS_MAT_X:
if (k) {
ierr = PSVectors_NHEP_Eigen_Some(ps,k,rnorm,PETSC_FALSE);CHKERRQ(ierr);
} else {
ierr = PSVectors_NHEP_Eigen_All(ps,PETSC_FALSE);CHKERRQ(ierr);
}
break;
case PS_MAT_Y:
if (k) {
ierr = PSVectors_NHEP_Eigen_Some(ps,k,rnorm,PETSC_TRUE);CHKERRQ(ierr);
} else {
ierr = PSVectors_NHEP_Eigen_All(ps,PETSC_TRUE);CHKERRQ(ierr);
}
break;
case PS_MAT_U:
case PS_MAT_VT:
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented yet");
break;
default:
SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PSSolve_NHEP"
PetscErrorCode PSSolve_NHEP(PS ps,PetscScalar *wr,PetscScalar *wi)
{
338,7 → 487,7
ps->nmeth = 1;
ps->ops->allocate = PSAllocate_NHEP;
ps->ops->view = PSView_NHEP;
//ps->ops->computevector = PSComputeVector_NHEP;
ps->ops->vectors = PSVectors_NHEP;
ps->ops->solve = PSSolve_NHEP;
ps->ops->sort = PSSort_NHEP;
ps->ops->cond = PSCond_NHEP;
/trunk/src/eps/interface/opts.c
838,12 → 838,12
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
PetscValidLogicalCollectiveEnum(eps,conv,2);
switch(conv) {
case EPS_CONV_EIG: eps->conv_func = EPSConvergedEigRelative; break;
case EPS_CONV_NORM: eps->conv_func = EPSConvergedNormRelative; break;
case EPS_CONV_ABS: eps->conv_func = EPSConvergedAbsolute; break;
case EPS_CONV_USER: break;
default:
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'conv' value");
case EPS_CONV_EIG: eps->conv_func = EPSConvergedEigRelative; break;
case EPS_CONV_NORM: eps->conv_func = EPSConvergedNormRelative; break;
case EPS_CONV_ABS: eps->conv_func = EPSConvergedAbsolute; break;
case EPS_CONV_USER: break;
default:
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'conv' value");
}
eps->conv = conv;
PetscFunctionReturn(0);