Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2602 → Rev 2603

/branches/slepc-3_2-proj/include/slepceps.h
347,7 → 347,6
extern PetscErrorCode EPSLanczosGetReorthog(EPS,EPSLanczosReorthogType*);
 
extern PetscErrorCode EPSBlzpackSetBlockSize(EPS,PetscInt);
extern PetscErrorCode EPSBlzpackSetInterval(EPS,PetscReal,PetscReal);
extern PetscErrorCode EPSBlzpackSetNSteps(EPS,PetscInt);
 
/*E
/branches/slepc-3_2-proj/include/private/epsimpl.h
139,7 → 139,7
extern PetscErrorCode EPSBasicArnoldi(EPS,PetscBool,PetscScalar*,PetscInt,Vec*,PetscInt,PetscInt*,Vec,PetscReal*,PetscBool*);
extern PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,Vec*,PetscInt,PetscInt*,Vec,PetscReal*,PetscBool*);
extern PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,Vec*,PetscInt,PetscInt*,Vec,PetscReal*,PetscBool*);
extern PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*,Vec*,PetscInt,PetscReal,PetscReal,PetscInt*,PetscScalar*);
extern PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscBool,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*,Vec*,PetscInt,PetscReal,PetscReal,PetscInt*,PetscScalar*);
extern PetscErrorCode EPSFullLanczos(EPS,PetscReal*,PetscReal*,Vec*,PetscInt,PetscInt*,Vec,PetscBool*);
extern PetscErrorCode EPSTranslateHarmonic(PetscInt,PetscScalar*,PetscInt,PetscScalar,PetscScalar,PetscScalar*,PetscScalar*);
extern PetscErrorCode EPSBuildBalance_Krylov(EPS);
/branches/slepc-3_2-proj/src/eps/impls/external/trlan/trlan.c
37,10 → 37,11
EPS_TRLAN *tr = (EPS_TRLAN *)eps->data;
 
PetscFunctionBegin;
tr->maxlan = PetscBLASIntCast(PetscMax(7,eps->nev+PetscMin(eps->nev,6)));
if (eps->ncv) {
if (eps->ncv<eps->nev) SETERRQ(((PetscObject)eps)->comm,1,"The value of ncv must be at least nev");
}
else eps->ncv = eps->nev;
else eps->ncv = tr->maxlan;
if (eps->mpd) { ierr = PetscInfo(eps,"Warning: parameter mpd ignored\n");CHKERRQ(ierr); }
if (!eps->max_it) eps->max_it = PetscMax(1000,eps->n);
55,7 → 56,6
SETERRQ(((PetscObject)eps)->comm,1,"Wrong value of eps->which");
 
tr->restart = 0;
tr->maxlan = PetscBLASIntCast(eps->nev+PetscMin(eps->nev,6));
if (tr->maxlan+1-eps->ncv<=0) { tr->lwork = PetscBLASIntCast(tr->maxlan*(tr->maxlan+10)); }
else { tr->lwork = PetscBLASIntCast(eps->nloc*(tr->maxlan+1-eps->ncv) + tr->maxlan*(tr->maxlan+10)); }
ierr = PetscMalloc(tr->lwork*sizeof(PetscReal),&tr->work);CHKERRQ(ierr);
/branches/slepc-3_2-proj/src/eps/impls/external/blzpack/blzpack.c
70,9 → 70,7
PetscErrorCode ierr;
PetscInt listor,lrstor,ncuv,k1,k2,k3,k4;
EPS_BLZPACK *blz = (EPS_BLZPACK *)eps->data;
PetscBool flg;
KSP ksp;
PC pc;
PetscBool issinv;
 
PetscFunctionBegin;
if (eps->ncv) {
85,21 → 83,26
 
if (!eps->ishermitian)
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
if (eps->which==EPS_ALL) {
if (eps->inta==0.0 && eps->intb==0.0) SETERRQ(((PetscObject)eps)->comm,1,"Must define a computational interval when using EPS_ALL");
blz->slice = 1;
}
ierr = PetscTypeCompare((PetscObject)eps->OP,STSINVERT,&issinv);CHKERRQ(ierr);
if (blz->slice || eps->isgeneralized) {
ierr = PetscTypeCompare((PetscObject)eps->OP,STSINVERT,&flg);CHKERRQ(ierr);
if (!flg)
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Shift-and-invert ST is needed for generalized problems or spectrum slicing");
ierr = STGetKSP(eps->OP,&ksp);CHKERRQ(ierr);
ierr = PetscTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);CHKERRQ(ierr);
if (!flg)
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Preonly KSP is needed for generalized problems or spectrum slicing");
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PetscTypeCompare((PetscObject)pc,PCCHOLESKY,&flg);CHKERRQ(ierr);
if (!flg)
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Cholesky PC is needed for generalized problems or spectrum slicing");
if (!issinv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Shift-and-invert ST is needed for generalized problems or spectrum slicing");
}
if (!eps->which) eps->which = EPS_SMALLEST_REAL;
if (eps->which!=EPS_SMALLEST_REAL)
if (blz->slice) {
if (eps->intb >= PETSC_MAX_REAL) { /* right-open interval */
if (eps->inta <= PETSC_MIN_REAL) SETERRQ(((PetscObject)eps)->comm,1,"The defined computational interval should have at least one of their sides bounded");
ierr = STSetDefaultShift(eps->OP,eps->inta);CHKERRQ(ierr);
}
else { ierr = STSetDefaultShift(eps->OP,eps->intb);CHKERRQ(ierr); }
}
if (!eps->which) {
if (issinv) eps->which = EPS_TARGET_REAL;
else eps->which = EPS_SMALLEST_REAL;
}
if ((issinv && eps->which!=EPS_TARGET_REAL && eps->which!=EPS_TARGET_MAGNITUDE && eps->which!=EPS_ALL) || (!issinv && eps->which!=EPS_SMALLEST_REAL))
SETERRQ(((PetscObject)eps)->comm,1,"Wrong value of eps->which");
 
k1 = PetscMin(eps->n,180);
114,6 → 117,7
 
if (blz->slice) lrstor = eps->nloc*(k2*4+k1*2+k4)+k3;
else lrstor = eps->nloc*(k2*4+k1)+k3;
lrstor*=10;
ierr = PetscFree(blz->rstor);CHKERRQ(ierr);
ierr = PetscMalloc((4+lrstor)*sizeof(PetscReal),&blz->rstor);CHKERRQ(ierr);
blz->rstor[3] = lrstor;
162,10 → 166,10
blz->rstor[1] = sigma; /* upper limit of eigenvalue interval */
} else {
sigma = 0.0;
blz->rstor[0] = blz->initial; /* lower limit of eigenvalue interval */
blz->rstor[1] = blz->final; /* upper limit of eigenvalue interval */
blz->rstor[0] = eps->inta; /* lower limit of eigenvalue interval */
blz->rstor[1] = eps->intb; /* upper limit of eigenvalue interval */
}
nneig = 0; /* no. of eigs less than sigma */
nneig = 0; /* no. of eigs less than sigma */
 
blz->istor[0] = PetscBLASIntCast(eps->nloc); /* number of rows of U, V, X*/
blz->istor[1] = PetscBLASIntCast(eps->nloc); /* leading dimension of U, V, X */
305,7 → 309,6
PetscFunctionBegin;
ierr = PetscFree(eps->data);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetInterval_C","",PETSC_NULL);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetNSteps_C","",PETSC_NULL);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
324,7 → 327,10
SETERRQ1(((PetscObject)eps)->comm,1,"Viewer type %s not supported for EPSBLZPACK",((PetscObject)viewer)->type_name);
}
ierr = PetscViewerASCIIPrintf(viewer," BLZPACK: block size=%d\n",blz->block_size);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer," BLZPACK: computational interval [%f,%f]\n",blz->initial,blz->final);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer," BLZPACK: maximum number of steps per run=%d\n",blz->nsteps);CHKERRQ(ierr);
if (blz->slice) {
ierr = PetscViewerASCIIPrintf(viewer," BLZPACK: computational interval [%f,%f]\n",eps->inta,eps->intb);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
335,10 → 341,7
PetscErrorCode ierr;
EPS_BLZPACK *blz = (EPS_BLZPACK *)eps->data;
PetscInt bs,n;
PetscReal interval[2];
PetscBool flg;
KSP ksp;
PC pc;
 
PetscFunctionBegin;
ierr = PetscOptionsHead("EPS BLZPACK Options");CHKERRQ(ierr);
351,23 → 354,6
ierr = PetscOptionsInt("-eps_blzpack_nsteps","Number of steps","EPSBlzpackSetNSteps",n,&n,&flg);CHKERRQ(ierr);
if (flg) {ierr = EPSBlzpackSetNSteps(eps,n);CHKERRQ(ierr);}
 
interval[0] = blz->initial;
interval[1] = blz->final;
n = 2;
ierr = PetscOptionsRealArray("-eps_blzpack_interval","Computational interval","EPSBlzpackSetInterval",interval,&n,&flg);CHKERRQ(ierr);
if (flg) {
if (n==1) interval[1]=interval[0];
ierr = EPSBlzpackSetInterval(eps,interval[0],interval[1]);CHKERRQ(ierr);
}
 
if (blz->slice || eps->isgeneralized) {
ierr = STSetType(eps->OP,STSINVERT);CHKERRQ(ierr);
ierr = STGetKSP(eps->OP,&ksp);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
}
 
ierr = PetscOptionsTail();CHKERRQ(ierr);
PetscFunctionReturn(0);
}
403,8 → 389,6
. -eps_blzpack_block_size - Sets the value of the block size
 
Level: advanced
 
.seealso: EPSBlzpackSetInterval()
@*/
PetscErrorCode EPSBlzpackSetBlockSize(EPS eps,PetscInt bs)
{
419,63 → 403,6
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "EPSBlzpackSetInterval_BLZPACK"
PetscErrorCode EPSBlzpackSetInterval_BLZPACK(EPS eps,PetscReal initial,PetscReal final)
{
PetscErrorCode ierr;
EPS_BLZPACK *blz = (EPS_BLZPACK*)eps->data;;
 
PetscFunctionBegin;
blz->initial = initial;
blz->final = final;
blz->slice = 1;
ierr = STSetShift(eps->OP,initial);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_END
 
#undef __FUNCT__
#define __FUNCT__ "EPSBlzpackSetInterval"
/*@
EPSBlzpackSetInterval - Sets the computational interval for the BLZPACK
package.
 
Collective on EPS
 
Input Parameters:
+ eps - the eigenproblem solver context
. initial - lower bound of the interval
- final - upper bound of the interval
 
Options Database Key:
. -eps_blzpack_interval - Sets the bounds of the interval (two values
separated by commas)
 
Note:
The following possibilities are accepted (see Blzpack user's guide for
details).
initial>final: start seeking for eigenpairs in the upper bound
initial<final: start in the lower bound
initial=final: run around a single value (no interval)
Level: advanced
 
.seealso: EPSBlzpackSetBlockSize()
@*/
PetscErrorCode EPSBlzpackSetInterval(EPS eps,PetscReal initial,PetscReal final)
{
PetscErrorCode ierr;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
PetscValidLogicalCollectiveReal(eps,initial,2);
PetscValidLogicalCollectiveReal(eps,final,3);
ierr = PetscTryMethod(eps,"EPSBlzpackSetInterval_C",(EPS,PetscReal,PetscReal),(eps,initial,final));CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "EPSBlzpackSetNSteps_BLZPACK"
PetscErrorCode EPSBlzpackSetNSteps_BLZPACK(EPS eps,PetscInt nsteps)
{
537,13 → 464,10
eps->ops->computevectors = EPSComputeVectors_Default;
 
blzpack->block_size = 3;
blzpack->initial = 0.0;
blzpack->final = 0.0;
blzpack->slice = 0;
blzpack->nsteps = 0;
 
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetBlockSize_C","EPSBlzpackSetBlockSize_BLZPACK",EPSBlzpackSetBlockSize_BLZPACK);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetInterval_C","EPSBlzpackSetInterval_BLZPACK",EPSBlzpackSetInterval_BLZPACK);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSBlzpackSetNSteps_C","EPSBlzpackSetNSteps_BLZPACK",EPSBlzpackSetNSteps_BLZPACK);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/branches/slepc-3_2-proj/src/eps/impls/external/blzpack/blzpackp.h
26,7 → 26,6
 
typedef struct {
PetscBLASInt block_size; /* block size */
PetscReal initial,final; /* computational interval */
PetscBLASInt slice; /* use spectrum slicing */
PetscBLASInt nsteps; /* maximum number of steps per run */
PetscBLASInt *istor;
/branches/slepc-3_2-proj/src/eps/impls/external/blopex/blopex.c
199,24 → 199,41
PetscErrorCode EPSSolve_BLOPEX(EPS eps)
{
EPS_BLOPEX *blopex = (EPS_BLOPEX *)eps->data;
int i,info,its;
int i,j,info,its,nconv;
double *lambdahist=PETSC_NULL,*residhist=PETSC_NULL;
PetscErrorCode ierr;
PetscFunctionBegin;
if (eps->numbermonitors>0) {
ierr = PetscMalloc(eps->ncv*(eps->max_it+1)*sizeof(double),&lambdahist);CHKERRQ(ierr);
ierr = PetscMalloc(eps->ncv*(eps->max_it+1)*sizeof(double),&residhist);CHKERRQ(ierr);
}
 
#if defined(PETSC_USE_COMPLEX)
info = lobpcg_solve_complex(blopex->eigenvectors,eps,OperatorAMultiVector,
eps->isgeneralized?eps:PETSC_NULL,eps->isgeneralized?OperatorBMultiVector:PETSC_NULL,
eps,Precond_FnMultiVector,blopex->Y,
blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
(komplex*)eps->eigr,PETSC_NULL,0,eps->errest,PETSC_NULL,0);
(komplex*)eps->eigr,lambdahist,eps->ncv,eps->errest,residhist,eps->ncv);
#else
info = lobpcg_solve_double(blopex->eigenvectors,eps,OperatorAMultiVector,
eps->isgeneralized?eps:PETSC_NULL,eps->isgeneralized?OperatorBMultiVector:PETSC_NULL,
eps,Precond_FnMultiVector,blopex->Y,
blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
eps->eigr,PETSC_NULL,0,eps->errest,PETSC_NULL,0);
eps->eigr,lambdahist,eps->ncv,eps->errest,residhist,eps->ncv);
#endif
if (info>0) SETERRQ1(((PetscObject)eps)->comm,PETSC_ERR_LIB,"Error in blopex (code=%d)",info);
 
if (eps->numbermonitors>0) {
for (i=0;i<its;i++) {
nconv = 0;
for (j=0;j<eps->ncv;j++) { if (residhist[j+i*eps->ncv]>eps->tol) break; else nconv++; }
ierr = EPSMonitor(eps,i,nconv,lambdahist+i*eps->ncv,eps->eigi,residhist+i*eps->ncv,eps->ncv);CHKERRQ(ierr);
}
ierr = PetscFree(lambdahist);CHKERRQ(ierr);
ierr = PetscFree(residhist);CHKERRQ(ierr);
}
 
eps->its = its;
eps->nconv = eps->ncv;
if (eps->OP->sigma != 0.0) {
/branches/slepc-3_2-proj/src/eps/impls/krylov/krylov.c
75,6 → 75,7
Input Parameters:
eps - the eigensolver; some error estimates are updated in eps->errest
issym - whether the projected problem is symmetric or not
trackall - whether all residuals must be computed
kini - initial value of k (the loop variable)
nits - number of iterations of the loop
S - Schur form of projected matrix (not referenced if issym)
91,7 → 92,7
Workspace:
work is workspace to store 5*nv scalars, nv booleans and nv reals (only if !issym)
*/
PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool issym,PetscInt kini,PetscInt nits,PetscScalar *S,PetscInt lds,PetscScalar *Q,Vec *V,PetscInt nv,PetscReal beta,PetscReal corrf,PetscInt *kout,PetscScalar *work)
PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool issym,PetscBool trackall,PetscInt kini,PetscInt nits,PetscScalar *S,PetscInt lds,PetscScalar *Q,Vec *V,PetscInt nv,PetscReal beta,PetscReal corrf,PetscInt *kout,PetscScalar *work)
{
PetscErrorCode ierr;
PetscInt k,marker;
129,7 → 130,7
ierr = (*eps->conv_func)(eps,re,im,resnorm,&eps->errest[k],eps->conv_ctx);CHKERRQ(ierr);
if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
if (iscomplex) { eps->errest[k+1] = eps->errest[k]; k++; }
if (marker!=-1 && !eps->trackall) break;
if (marker!=-1 && !trackall) break;
}
if (marker!=-1) k = marker;
*kout = k;
/branches/slepc-3_2-proj/src/eps/impls/krylov/arnoldi/arnoldi.c
416,7 → 416,7
ierr = EPSProjectedArnoldi(eps,H,eps->ncv,U,nv);CHKERRQ(ierr);
 
/* Check convergence */
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,H,eps->ncv,U,eps->V,nv,beta,corrf,&k,work);CHKERRQ(ierr);
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->trackall,eps->nconv,nv-eps->nconv,H,eps->ncv,U,eps->V,nv,beta,corrf,&k,work);CHKERRQ(ierr);
 
ierr = EPSUpdateVectors(eps,nv,eps->V,eps->nconv,PetscMin(k+1,nv),U,nv,Hcopy,eps->ncv);CHKERRQ(ierr);
eps->nconv = k;
/branches/slepc-3_2-proj/src/eps/impls/krylov/krylovschur/ks-symm.c
203,7 → 203,7
ierr = EPSProjectedKSSym(eps,nv,l,a,b,eps->eigr+eps->nconv,Q,work,iwork);CHKERRQ(ierr);
 
/* Check convergence */
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv,PETSC_NULL,nv,Q,eps->V+eps->nconv,nv,beta,1.0,&k,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->trackall,eps->nconv,nv,PETSC_NULL,nv,Q,eps->V+eps->nconv,nv,beta,1.0,&k,PETSC_NULL);CHKERRQ(ierr);
if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
if (k >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
/branches/slepc-3_2-proj/src/eps/impls/krylov/krylovschur/ks-harm.c
196,7 → 196,7
ierr = EPSProjectedKSNonsym(eps,l,S,eps->ncv,Q,nv);CHKERRQ(ierr);
 
/* Check convergence */
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,S,eps->ncv,Q,eps->V,nv,beta,PetscSqrtReal(1.0+gnorm),&k,work);CHKERRQ(ierr);
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->trackall,eps->nconv,nv-eps->nconv,S,eps->ncv,Q,eps->V,nv,beta,PetscSqrtReal(1.0+gnorm),&k,work);CHKERRQ(ierr);
if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
if (k >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
/branches/slepc-3_2-proj/src/eps/impls/krylov/krylovschur/ks-slice.c
38,49 → 38,44
 
extern PetscErrorCode EPSProjectedKSSym(EPS,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscScalar*,PetscScalar*,PetscReal*,PetscInt*);
 
/**/
PetscInt allKs,def,deg,db;
 
/* Type of data characterizing a shift (place from where an eps is applied) */
typedef struct _n_shift *shift;
struct _n_shift{
PetscReal value;
PetscInt inertia;
PetscBool comp[2]; //
shift neighb[2];//
PetscInt index; //
PetscInt neigs; //
PetscReal ext[2]; //
PetscInt nsch[2]; //
PetscReal pert; //
PetscBool deg; //not used
PetscBool comp[2]; /* Shows completion of subintervals (left and right) */
shift neighb[2];/* Adjacent shifts */
PetscInt index;/* Index in eig where found values are stored */
PetscInt neigs; /* Number of values found */
PetscReal ext[2]; /* Limits for accepted values */
PetscInt nsch[2]; /* Number of missing values for each subinterval */
PetscInt nconv[2]; /* Converged on each side (accepted or not)*/
};
 
/* Type of data for storing the state of spectrum slicing*/
struct _n_SR{
PetscReal int0,int1; // extremes of the interval
PetscInt dir; // determines the order of values in eig (+1 incr, -1 decr)
PetscBool hasEnd; // tells whether the interval has an end
PetscReal int0,int1; /* Extremes of the interval */
PetscInt dir; /* Determines the order of values in eig (+1 incr, -1 decr) */
PetscBool hasEnd; /* Tells whether the interval has an end */
PetscInt inertia0,inertia1;
Vec *V;
PetscScalar *eig;
PetscInt *perm;// permutation for keeping the eigenvalues in order
PetscInt numEigs; // number of eigenvalues in the interval
PetscScalar *eig,*eigi,*monit;
PetscReal *errest;
PetscInt *perm;/* Permutation for keeping the eigenvalues in order */
PetscInt numEigs; /* Number of eigenvalues in the interval */
PetscInt indexEig;
shift sPres; // present shift
shift *pending;//pending shifts array
PetscInt nPend;// number of pending shifts
PetscInt maxPend;// size of "pending" array
Vec *VDef; // vector for deflation
PetscInt *idxDef;//for deflation
shift sPres; /* Present shift */
shift *pending;/* Pending shifts array */
PetscInt nPend;/* Number of pending shifts */
PetscInt maxPend;/* Size of "pending" array */
Vec *VDef; /* Vector for deflation */
PetscInt *idxDef;/* For deflation */
PetscInt nMAXCompl;
PetscInt iterCompl;
PetscInt itsKs;
PetscInt nShifts;//number of computed shifts
shift s0;// initial shift
PetscReal tolDeg;
PetscInt nDeg;//number of values coinciding with a shift
PetscInt defMin; //minimum amount of values for deflation
PetscInt itsKs; /* Krylovschur restarts */
PetscInt nleap;
shift s0;/* Initial shift */
};
typedef struct _n_SR *SR;
 
109,12 → 104,11
s->comp[1] = PETSC_FALSE;
s->index = -1;
s->neigs = 0;
s->deg = PETSC_FALSE;
s->nconv[0] = s->nconv[1] = 0;
s->nsch[0] = s->nsch[1]=0;
// inserts in the stack of pending shifts
// if needed, the array is resized
/* Inserts in the stack of pending shifts */
/* If needed, the array is resized */
if(sr->nPend >= sr->maxPend){
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"resizing pending shifts array\n");CHKERRQ(ierr);}
sr->maxPend *= 2;
ierr = PetscMalloc((sr->maxPend)*sizeof(shift),&pending2);CHKERRQ(ierr);
for(i=0;i < sr->nPend; i++)pending2[i] = sr->pending[i];
122,7 → 116,6
sr->pending = pending2;
}
sr->pending[sr->nPend++]=s;
sr->nShifts++;
PetscFunctionReturn(0);
}
 
150,15 → 143,15
eps->target = sr->sPres->value;
eps->nconv = 0;
eps->reason = EPS_CONVERGED_ITERATING;
eps->its =0;
eps->its = 0;
}else sr->sPres = PETSC_NULL;
PetscFunctionReturn(0);
}
/*
Symmetric KrylovSchur adapted to spectrum slicing:
allows searching an specific amount of eigenvalues in the subintervals left and right.
returns whether the search has succeed
Allows searching an specific amount of eigenvalues in the subintervals left and right.
Returns whether the search has succeeded
*/
#undef __FUNCT__
#define __FUNCT__ "EPSKrylovSchur_Slice"
168,20 → 161,18
PetscInt i,conv,k,l,lds,lt,nv,m,*iwork,p,j;
Vec u=eps->work[0];
PetscScalar *Q;
PetscReal *a,*b,*work,beta,*Qreal,rtmp;//
PetscReal *a,*b,*work,beta,rtmp;
PetscBool breakdown;
PetscInt count0,count1;//nconv_prev;
PetscInt count0,count1;
PetscReal theta,lambda;
shift sPres;
PetscBool complIterating;/* shows whether iterations are made for completion */
PetscBool sch0,sch1;//shows whether values are looked after on each side
PetscInt iterCompl=0,n0,n1;
//PetscReal res;
 
PetscBool complIterating;/* Shows whether iterations are made for completion */
PetscBool sch0,sch1;/* Shows whether values are looked after on each side */
PetscInt iterCompl=0,n0,n1,aux,auxc;
SR sr;
 
PetscFunctionBegin;
/* spectrum slicing data */
/* Spectrum slicing data */
sr = (SR)eps->data;
sPres = sr->sPres;
complIterating =PETSC_FALSE;
193,15 → 184,19
lt = PetscMin(eps->nev+eps->mpd,eps->ncv);
ierr = PetscMalloc(lt*sizeof(PetscReal),&a);CHKERRQ(ierr);
ierr = PetscMalloc(lt*sizeof(PetscReal),&b);CHKERRQ(ierr);
count0=0;count1=0; // found on both sides
count0=0;count1=0; /* Found on both sides */
 
/* filling in values for the monitor */
for(i=0;i<sr->indexEig;i++){
sr->monit[i]=1.0/(sr->eig[i] - sPres->value);
}
/* Get the starting Lanczos vector */
ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);
l = 0;
/* Restart loop */
while (eps->reason == EPS_CONVERGED_ITERATING) {
eps->its++;
 
eps->its++; sr->itsKs++;
/* Compute an nv-step Lanczos factorization */
m = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 
211,52 → 206,36
 
/* Solve projected problem and compute residual norm estimates */
ierr = EPSProjectedKSSym(eps,nv,l,a,b,eps->eigr+eps->nconv,Q,work,iwork);CHKERRQ(ierr);
/* Check convergence */
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv,PETSC_NULL,nv,Q,eps->V+eps->nconv,nv,beta,1.0,&k,PETSC_NULL);CHKERRQ(ierr);
if(allKs ==1){//option for accepting all converging values
Qreal = (PetscReal*)Q;//
conv=k=j=0;
for(i=0;i<nv;i++)if(eps->errest[eps->nconv+i] < eps->tol)conv++;
for(i=0;i<nv;i++){
if(eps->errest[eps->nconv+i] < eps->tol){
iwork[j++]=i;
}else iwork[conv+k++]=i;
}
for(i=0;i<nv;i++) a[i]=PetscRealPart(eps->eigr[eps->nconv+i]);
for(i=0;i<nv;i++){
eps->eigr[eps->nconv+i] = a[iwork[i]];
}
for( i=0;i<nv;i++){
p=iwork[i];
if(p!=i){
j=i+1;
while(iwork[j]!=i)j++;
iwork[j]=p;iwork[i]=i;
for(k=0;k<nv;k++){
rtmp=Qreal[k+p*nv];Qreal[k+p*nv]=Qreal[k+i*nv];Qreal[k+i*nv]=rtmp;
//rtmp=Q[k+p*nv];Q[k+p*nv]=Q[k+i*nv];Q[k+i*nv]=rtmp;
}
}
}
k=eps->nconv+conv;
/* Residual */
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,PETSC_TRUE,eps->nconv,nv,PETSC_NULL,nv,Q,eps->V+eps->nconv,nv,beta,1.0,&k,PETSC_NULL);CHKERRQ(ierr);
/* Check convergence */
conv=k=j=0;
for(i=0;i<nv;i++)if(eps->errest[eps->nconv+i] < eps->tol)conv++;
for(i=0;i<nv;i++){
if(eps->errest[eps->nconv+i] < eps->tol){
iwork[j++]=i;
}else iwork[conv+k++]=i;
}
/*checking proximity to an eigenvalue*/
 
if(deg==1){
for(i=0;i < k; i++){
theta = PetscRealPart(eps->eigr[i]);
if(PetscAbsReal(theta*sPres->value*eps->tol*10)>1){
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"DEGENERATED SHIFT\n");CHKERRQ(ierr);}
sr->nDeg++;
sPres->deg = PETSC_TRUE;
}else break;
for(i=0;i<nv;i++){
a[i]=PetscRealPart(eps->eigr[eps->nconv+i]);
b[i]=eps->errest[eps->nconv+i];
}
}
/*
if(deg == 1 && sr->nDeg > 0){
eps->reason = EPS_CONVERGED_TOL;
}else{
*/
for(i=0;i<nv;i++){
eps->eigr[eps->nconv+i] = a[iwork[i]];
eps->errest[eps->nconv+i] = b[iwork[i]];
}
for( i=0;i<nv;i++){
p=iwork[i];
if(p!=i){
j=i+1;
while(iwork[j]!=i)j++;
iwork[j]=p;iwork[i]=i;
for(k=0;k<nv;k++){
rtmp=Q[k+p*nv];Q[k+p*nv]=Q[k+i*nv];Q[k+i*nv]=rtmp;
}
}
}
k=eps->nconv+conv;
/* Checking values obtained for completing */
count0=count1=0;
for(i=0;i<k;i++){
266,19 → 245,19
if( ((sr->dir)*theta > 0) && ((sr->dir)*(sPres->ext[1] - lambda) > 0))count1++;
}
/* checks completion */
/* Checks completion */
if( (!sch0||count0 >= sPres->nsch[0]) && (!sch1 ||count1 >= sPres->nsch[1]) ) {
eps->reason = EPS_CONVERGED_TOL;
}else {
if(!complIterating && eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
if(complIterating){
if(--iterCompl <= 0) eps->reason = EPS_DIVERGED_ITS;
}else if (k >= eps->nev) {//eps->reason = EPS_CONVERGED_TOL;
}else if (k >= eps->nev) {
n0 = sPres->nsch[0]-count0;
n1 = sPres->nsch[1]-count1;
if( sr->iterCompl>0 && ( (n0>0&& n0<= sr->nMAXCompl)||(n1>0&&n1<=sr->nMAXCompl) )){
/* Iterating for completion*/
complIterating = PETSC_TRUE;
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"iterating for completion nMAXComp=%d\n",sr->nMAXCompl);CHKERRQ(ierr);}
if(n0 >sr->nMAXCompl)sch0 = PETSC_FALSE;
if(n1 >sr->nMAXCompl)sch1 = PETSC_FALSE;
iterCompl = sr->iterCompl;
313,17 → 292,27
if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
ierr = VecAXPBY(eps->V[k+l],1.0/beta,0.0,u);CHKERRQ(ierr);
}
 
ierr = EPSMonitor(eps,eps->its,k,eps->eigr,eps->eigi,eps->errest,nv+eps->nconv);CHKERRQ(ierr);
//nconv_prev = eps->nconv;//
aux = auxc = 0;
for(i=0;i<nv+eps->nconv;i++){
theta = PetscRealPart(eps->eigr[i]);
lambda = sPres->value + 1/theta;
if( ((sr->dir)*(lambda - sPres->ext[0]) > 0)&& ((sr->dir)*(sPres->ext[1] - lambda) > 0)){
sr->monit[sr->indexEig+aux]=eps->eigr[i];
sr->errest[sr->indexEig+aux]=eps->errest[i];
aux++;
if(eps->errest[i] < eps->tol)auxc++;
}
}
ierr = EPSMonitor(eps,eps->its,auxc+sr->indexEig,sr->monit,sr->eigi,sr->errest,sr->indexEig+aux);CHKERRQ(ierr);
eps->nconv = k;
}
/* check for completion */
/* Check for completion */
for(i=0;i< eps->nconv; i++){
if( (sr->dir)*PetscRealPart(eps->eigr[i])>0 )sPres->nconv[1]++;
else sPres->nconv[0]++;
}
sPres->comp[0] = (count0 >= sPres->nsch[0])?PETSC_TRUE:PETSC_FALSE;
sPres->comp[1] = (count1 >= sPres->nsch[1])?PETSC_TRUE:PETSC_FALSE;
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD," found count0=%d(of %d) and count1=%d(of %d)\n",count0,sPres->nsch[0],count1,sPres->nsch[1]);CHKERRQ(ierr);}
sr->itsKs += eps->its;
 
ierr = PetscFree(Q);CHKERRQ(ierr);
ierr = PetscFree(a);CHKERRQ(ierr);
ierr = PetscFree(b);CHKERRQ(ierr);
333,12 → 322,12
}
 
/*
obtains value of subsequent shift
Obtains value of subsequent shift
*/
#undef __FUNCT__
#define __FUNCT__ "EPSGetNewShiftValue"
static PetscErrorCode EPSGetNewShiftValue(EPS eps,PetscInt side,PetscReal *newS){
PetscReal lambda,d_prev;//pert,d,d_avg;
PetscReal lambda,d_prev;
PetscInt i,idxP;
SR sr;
shift sPres;
346,75 → 335,72
PetscFunctionBegin;
sr = (SR)eps->data;
sPres = sr->sPres;
/*
pert = 0;
if(sPres->neigs >0){
idxP=0;//number of computed eigenvalues previous to sPres->value
for(i=sPres->index;i< sPres->index + sPres->neigs;i++){
lambda = PetscRealPart(sr->eig[i]);
if((sr->dir)*(lambda - sPres->value) < 0)idxP++;
else break;
}
//middle point between shift and previous/posterior value
pert = PetscAbs(sr->eig[sPres->index+idxP]- sr->sPres->value)/2;
}
*/
if( sPres->neighb[side]){
/* completing a previous interval */
*newS=(sPres->value + sPres->neighb[side]->value)/2;
}else{ //(only for side=1). creating a new interval.
if(sPres->neigs==0){// no value has been accepted
/* Completing a previous interval */
if(!sPres->neighb[side]->neighb[side] && sPres->neighb[side]->nconv[side]==0){ /* One of the ends might be too far from eigenvalues */
if(side) *newS = (sPres->value + PetscRealPart(sr->eig[sr->perm[sr->indexEig-1]]))/2;
else *newS = (sPres->value + PetscRealPart(sr->eig[sr->perm[0]]))/2;
}else *newS=(sPres->value + sPres->neighb[side]->value)/2;
}else{ /* (Only for side=1). Creating a new interval. */
if(sPres->neigs==0){/* No value has been accepted*/
if(sPres->neighb[0]){
// multiplying by 10 the previous distance
/* Multiplying by 10 the previous distance */
*newS = sPres->value + 10*(sr->dir)*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
}else {//first shift
sr->nleap++;
/* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
if( !sr->hasEnd && sr->nleap > 5)SETERRQ(((PetscObject)eps)->comm,1,"Unable to compute the wanted eigenvalues with open interval");
}else {/* First shift */
if(eps->nconv != 0){
//unaccepted values give information for next shift
idxP=0;//number of values left from shift
/* Unaccepted values give information for next shift */
idxP=0;/* Number of values left from shift */
for(i=0;i<eps->nconv;i++){
lambda = PetscRealPart(eps->eigr[i]);
if( (sr->dir)*(lambda - sPres->value) <0)idxP++;
else break;
}
//avoiding substraction of eigenvalues (might be the same).
/* Avoiding subtraction of eigenvalues (might be the same).*/
if(idxP>0){
d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[0]))/(idxP+0.3);
}else {
d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[eps->nconv-1]))/(eps->nconv+0.3);
}
*newS = sPres->value + ((sr->dir)*d_prev*eps->nev)/2;
}else{//no values found, no information for next shift
// changing the end of the interval
}else{/* No values found, no information for next shift */
SETERRQ(((PetscObject)eps)->comm,1,"First shift renders no information");
}
}
}else{// accepted values found
//average distance of values in previous subinterval
}else{/* Accepted values found */
sr->nleap = 0;
/* Average distance of values in previous subinterval */
shift s = sPres->neighb[0];
while(s && PetscAbs(s->inertia - sPres->inertia)==0){
s = s->neighb[0];//looking for previous shifts with eigenvalues within
s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
}
if(s){
d_prev = PetscAbsReal( (sPres->value - s->value)/(sPres->inertia - s->inertia));
}else{//first shift. average distance obtained with values in this shift
d_prev = PetscAbsReal( PetscRealPart(sr->eig[sPres->index+sPres->neigs-1]) - sPres->value)/(sPres->neigs+0.3);
}else{/* First shift. Average distance obtained with values in this shift */
/* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
if( (sr->dir)*(sr->eig[0]-sPres->value)>0 && PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0]))/PetscRealPart(sr->eig[0])) > PetscSqrtReal(eps->tol) ){
d_prev = PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0])))/(sPres->neigs+0.3);
}else{
d_prev = PetscAbsReal( PetscRealPart(sr->eig[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
}
}
// average distance is used for next shift by adding it to value on the right or to shift
/* Average distance is used for next shift by adding it to value on the right or to shift */
if( (sr->dir)*(PetscRealPart(sr->eig[sPres->index + sPres->neigs -1]) - sPres->value) >0){
*newS = PetscRealPart(sr->eig[sPres->index + sPres->neigs -1])+ ((sr->dir)*d_prev*(eps->nev))/2;
}else{//last accepted value is on the left of shift. adding to shift.
}else{/* Last accepted value is on the left of shift. Adding to shift */
*newS = sPres->value + ((sr->dir)*d_prev*(eps->nev))/2;
}
//}
}
//end of interval can not be surpassed
if(sr->hasEnd && ((sr->dir)*(*newS - sr->int1) > 0))*newS=sr->int1;
/* End of interval can not be surpassed */
if((sr->dir)*( sr->int1 - *newS) < 0) *newS = sr->int1;
}//of neighb[side]==null
PetscFunctionReturn(0);
}
 
/*
function for sorting an array of real values
Function for sorting an array of real values
*/
#undef __FUNCT__
#define __FUNCT__ "sortRealEigenvalues"
425,7 → 411,7
PetscFunctionBegin;
if(!prev) for (i=0; i < nr; i++) { perm[i] = i; }
/* insertion sort */
/* Insertion sort */
for (i=1; i < nr; i++) {
re = PetscRealPart(r[perm[i]]);
j = i-1;
442,7 → 428,7
PetscErrorCode EPSStoreEigenpairs(EPS eps)
{
PetscErrorCode ierr;
PetscReal lambda,error;
PetscReal lambda,err,norm;
PetscInt i,count;
PetscBool cond;
SR sr;
453,32 → 439,32
sPres = sr->sPres;
sPres->index = sr->indexEig;
count = sr->indexEig;
error=0;
/* backtransform */
/* Backtransform */
for(i=0;i < eps->nconv; i++) eps->eigr[i] = sPres->value + 1.0/(eps->eigr[i]);
/* sort eigenvalues */
/* Sort eigenvalues */
ierr = sortRealEigenvalues(eps->eigr,eps->perm,eps->nconv,PETSC_FALSE,sr->dir);
/* values stored in global array */
// condition for avoiding comparing whith a non-existing end.
/* Values stored in global array */
/* Condition for avoiding comparing with a non-existing end */
cond = (!sPres->neighb[1] && !sr->hasEnd)?PETSC_TRUE:PETSC_FALSE;
for( i=0; i < eps->nconv ;i++ ){
lambda = PetscRealPart(eps->eigr[eps->perm[i]]);
if(db>1){ierr = EPSComputeRelativeError(eps,eps->perm[i],&error);CHKERRQ(ierr);}
if( ( ((sr->dir)*(lambda - sPres->ext[0]) > 0) && ( cond || ((sr->dir)*(lambda - sPres->ext[1]) < 0)) ) ){
if(count>=sr->numEigs){//Error found
ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of reserved values exceeded lambda=%.14g\n",lambda);
break;
err = eps->errest[eps->perm[i]];
if( (sr->dir)*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0 ){/* Valid value */
if(count>=sr->numEigs){/* Error found */
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Unexpected error in Spectrum Slicing");
}
sr->eig[count] = lambda;
ierr = VecCopy(eps->V[eps->perm[i]], sr->V[count]);CHKERRQ(ierr);
if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"i=%d perm[i]=%d lambda=%.14g error=%g indexEig=%d\n",i,eps->perm[i],lambda,error,count);CHKERRQ(ierr);}
sr->errest[count] = err;
/* Purification */
ierr = STApply(eps->OP,eps->V[eps->perm[i]],sr->V[count]);CHKERRQ(ierr);
ierr = IPNorm(eps->ip,sr->V[count],&norm);CHKERRQ(ierr);
ierr = VecScale(sr->V[count],1.0/norm);CHKERRQ(ierr);
count++;
}else if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"i=%d perm[i]=%d lambda=%.14g NOT VALID\n",i,eps->perm[i],lambda);CHKERRQ(ierr);}
}
}
sPres->neigs = count - sr->indexEig;
if(db>=1){PetscPrintf(PETSC_COMM_WORLD," stored between %d and %d\n",sr->indexEig,count);CHKERRQ(ierr);}
sr->indexEig = count;
/* Global ordering array updating */
ierr = sortRealEigenvalues(sr->eig,sr->perm,count,PETSC_TRUE,sr->dir);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
487,41 → 473,36
#define __FUNCT__ "EPSLookForDeflation"
PetscErrorCode EPSLookForDeflation(EPS eps)
{
PetscErrorCode ierr;
PetscReal val,lambda,lambda2;
PetscReal val;
PetscInt i,count0=0,count1=0;
shift sPres;
PetscInt ini,fin,defMin,k,idx0,idx1;
PetscBool consec;
PetscInt ini,fin,k,idx0,idx1;
SR sr;
 
PetscFunctionBegin;
sr = (SR)(eps->data);
sPres = sr->sPres;
defMin = sr->defMin;
 
if(sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
else ini = 0;
fin = sr->indexEig;
// selection of ends for searching new values
// later modified with version def=0
if(!sPres->neighb[0]) sPres->ext[0] = sr->int0;//first shift
/* Selection of ends for searching new values */
if(!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
else sPres->ext[0] = sPres->neighb[0]->value;
if(!sPres->neighb[1]) {
if(sr->hasEnd) sPres->ext[1] = sr->int1;
else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
}else sPres->ext[1] = sPres->neighb[1]->value;
if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"ext0=%g ext1=%g\n",sPres->ext[0],sPres->ext[1]);CHKERRQ(ierr);}
//selection of values between right and left ends
/* Selection of values between right and left ends */
for(i=ini;i<fin;i++){
val=PetscRealPart(sr->eig[sr->perm[i]]);
//values to the right of left shift
/* Values to the right of left shift */
if( (sr->dir)*(val - sPres->ext[1]) < 0 ){
if((sr->dir)*(val - sPres->value) < 0)count0++;
else count1++;
}else break;
}
// the number of values on each side are found
/* The number of values on each side are found */
if(sPres->neighb[0])
sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
else sPres->nsch[0] = 0;
530,65 → 511,14
sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);
//completing vector of indexes for deflation
if(def==0 && !sPres->neighb[1]){//new interval && no deflation
if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"def=0 y neig1=null\n");CHKERRQ(ierr);}
k=0;
for(i=fin-1;i>ini;i--){
k++;
lambda = PetscRealPart(sr->eig[sr->perm[i]]);
lambda2 = PetscRealPart(sr->eig[sr->perm[i-1]]);
if( PetscAbsReal((lambda - lambda2)/lambda) > sr->tolDeg){//relative tolerance
break;
}
}
// if i!=ini values for i and i-1 more than toldeg apart
if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"lookDef i=%d ini=%d\n",i,ini);CHKERRQ(ierr);}
if(i<=ini){
sPres->ext[0] = sPres->value;
}else{//middle point
sPres->ext[0] = (PetscRealPart(sr->eig[sr->perm[i]])+PetscRealPart(sr->eig[sr->perm[i-1]]))/2;
}
idx0=ini+count0-k;
idx1=ini+count0;
if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"ext0=%g ext1=%g idx0=%d idx1=%d count0=%d k=%d\n",sPres->ext[0],sPres->ext[1],idx0,idx1,count0,k);CHKERRQ(ierr);}
}else{ //completing a subinterval or without deflation
k = PetscMax(0,defMin-count0);
idx0 = PetscMax(0,ini-k);
if(def==0 && sPres->nsch[0]==0){//no deflation towards 0
idx0 = ini + count0;
sPres->ext[0] = sPres->value;
}
k = PetscMax(0,defMin-count1);
idx1 = PetscMin(sr->indexEig,ini+count0+count1+k);
if(def==0 && sPres->nsch[1]==0){//no deflation towards 1
idx1 = ini + count0;
sPres->ext[1] = sPres->value;
}
}
/* Completing vector of indexes for deflation */
idx0 = ini;
idx1 = ini+count0+count1;
k=0;
for(i=idx0;i<idx1;i++)sr->idxDef[k++]=sr->perm[i];
///// info
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD," deflated %d in [0] and %d in [1]",count0,count1);CHKERRQ(ierr);}
/////
 
// check for consecutives
consec=PETSC_TRUE;
for(i=1;i<k;i++)if(sr->idxDef[i]!=sr->idxDef[i-1]+1){consec = PETSC_FALSE; break;}
// if not consecutives, copied in array
//if(consec){
// V o which
//else{
for(i=0;i<k;i++)sr->VDef[i]=sr->V[sr->idxDef[i]];
eps->DS = sr->VDef;
//}
eps->nds = k;
//////info
if(db>=1){
if(consec){ ierr = PetscPrintf(PETSC_COMM_WORLD," (%d consecutive values)\n",k);CHKERRQ(ierr);}
else{ ierr = PetscPrintf(PETSC_COMM_WORLD," (%d non consecutive values)\n",k);CHKERRQ(ierr);}
}
PetscFunctionReturn(0);
}
 
597,71 → 527,48
PetscErrorCode EPSSolve_KrylovSchur_Slice(EPS eps)
{
PetscErrorCode ierr;
PetscInt i,j;
PetscInt imax,jmax;
PetscInt i;
PetscReal newS;
KSP ksp;
PC pc;
Mat B,F;
PetscScalar *eigi;
Vec t,w;
Mat F;
PetscReal *errest_left;
Vec t;
SR sr;
PetscReal orthMax;
PetscScalar inerd;
double t1,t2;
PetscFunctionBegin;
eps->trackall = PETSC_TRUE;
allKs = 0;
def = 1;
deg=0;
db = 0;
ierr = PetscOptionsGetInt(PETSC_NULL,"-db",&db,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-deg",&deg,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-def",&def,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-allKs",&allKs,PETSC_NULL);CHKERRQ(ierr);
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"Options: allKs=%d, def=%d, deg=%d \n",allKs,def,deg);CHKERRQ(ierr);}
ierr = PetscMalloc(sizeof(struct _n_SR),&sr);CHKERRQ(ierr);
eps->data = sr;
sr->tolDeg = PetscSqrtReal(eps->tol);//default
ierr = PetscOptionsGetReal(PETSC_NULL,"-toldeg",&sr->tolDeg,PETSC_NULL);CHKERRQ(ierr);
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"toldeg=%g\n",sr->tolDeg);CHKERRQ(ierr);}
sr->defMin = 0;
ierr = PetscOptionsGetInt(PETSC_NULL,"-defMin",&sr->defMin,PETSC_NULL);CHKERRQ(ierr);
if(def==0)sr->defMin =0;
//checking presence of ends and finding direction
sr->itsKs = 0;
sr->nleap = 0;
sr->nMAXCompl = eps->nev/4;
sr->iterCompl = eps->max_it/4;
/* Checking presence of ends and finding direction */
if( eps->inta > PETSC_MIN_REAL){
sr->int0 = eps->inta;
sr->int1 = eps->intb;
sr->dir = 1;
if(eps->intb >= PETSC_MAX_REAL){ /* right-open interval */
if(eps->intb >= PETSC_MAX_REAL){ /* Right-open interval */
sr->hasEnd = PETSC_FALSE;
sr->inertia1 = eps->n;
}else sr->hasEnd = PETSC_TRUE;
}else{ /* left-open interval */
}else{ /* Left-open interval */
sr->int0 = eps->intb;
sr->int1 = eps->inta;
sr->dir = -1;
sr->hasEnd = PETSC_FALSE;
sr->inertia1 = 0;
}
if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"dir=%d int0=%g\n",sr->dir,sr->int0);CHKERRQ(ierr);}
sr->nMAXCompl = 0;
ierr = PetscOptionsGetInt(PETSC_NULL,"-comp",&sr->nMAXCompl,PETSC_NULL);
sr->iterCompl = sr->nMAXCompl+5;//=======
//i = PetscMin(eps->mpd,eps->ncv);//=======
//ierr = PetscMalloc(i*sizeof(PetscReal),&sr->aprox);CHKERRQ(ierr);//======
// array of pending shifts
sr->maxPend = 100;//initial size;
/* Array of pending shifts */
sr->maxPend = 100;/* Initial size */
ierr = PetscMalloc((sr->maxPend)*sizeof(shift),&sr->pending);CHKERRQ(ierr);
if(sr->hasEnd){
ierr = STGetKSP(eps->OP, &ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr);
/* not looking for values in b (just inertia).*/
/* Not looking for values in b (just inertia).*/
ierr = MatGetInertia(F,&sr->inertia1,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
}
sr->nShifts = 0;
sr->nPend = 0;
ierr = EPSCreateShift(eps,sr->int0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSExtractShift(eps);
669,111 → 576,80
sr->inertia0 = sr->s0->inertia;
sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
sr->indexEig = 0;
sr->itsKs = 0;
sr->nDeg = 0;
if(db>=1){
ierr = PetscPrintf(PETSC_COMM_WORLD,"dir=%d inertia in 0(=%g) %d and in 1(=%g) %d\n",sr->dir,sr->int0,sr->inertia0,sr->int1,sr->inertia1);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"numEigs=%d\n\n",sr->numEigs);
}
/* only with eigenvalues present in the interval ...*/
/* Only with eigenvalues present in the interval ...*/
if(sr->numEigs==0){
eps->reason = EPS_CONVERGED_TOL;
ierr = PetscFree(sr->s0);CHKERRQ(ierr);
ierr = PetscFree(sr->pending);CHKERRQ(ierr);
ierr = PetscFree(sr);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
/* memory reservation for eig, V and perm */
/* Memory reservation for eig, V and perm */
ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eig);CHKERRQ(ierr);
ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&eigi);CHKERRQ(ierr);
for(i=0;i<sr->numEigs;i++)eigi[i]=0;
ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eigi);CHKERRQ(ierr);
ierr = PetscMalloc((sr->numEigs+eps->ncv) *sizeof(PetscReal),&sr->errest);CHKERRQ(ierr);
ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscReal),&errest_left);CHKERRQ(ierr);
ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscReal),&sr->monit);CHKERRQ(ierr);
for(i=0;i<sr->numEigs;i++){sr->eigi[i]=0;sr->eig[i] = 0;}
for(i=0;i<sr->numEigs+eps->ncv;i++){errest_left[i]=0;sr->errest[i]=0;sr->monit[i]=0;}
ierr = VecCreateMPI(((PetscObject)eps)->comm,eps->nloc,PETSC_DECIDE,&t);CHKERRQ(ierr);
ierr = VecDuplicateVecs(t,sr->numEigs,&sr->V);CHKERRQ(ierr);
ierr = VecDestroy(&t);CHKERRQ(ierr);
// vector for maintaining order of eigenvalues
/* Vector for maintaining order of eigenvalues */
ierr = PetscMalloc((sr->numEigs)*sizeof(PetscInt),&sr->perm);CHKERRQ(ierr);
for(i=0;i< sr->numEigs;i++)sr->perm[i]=i;
// vectors for deflation
ierr = PetscMalloc((sr->numEigs+sr->defMin)*sizeof(PetscInt),&sr->idxDef);CHKERRQ(ierr);
/* Vectors for deflation */
ierr = PetscMalloc((sr->numEigs)*sizeof(PetscInt),&sr->idxDef);CHKERRQ(ierr);
ierr = PetscMalloc((sr->numEigs)*sizeof(Vec),&sr->VDef);CHKERRQ(ierr);
sr->indexEig = 0;
 
t1 = MPI_Wtime();
while(sr->sPres){
 
//////////info
if(db>=1){
if(sr->sPres->neighb[1]){ierr = PetscPrintf(PETSC_COMM_WORLD,"Completing ");CHKERRQ(ierr);}
else {ierr = PetscPrintf(PETSC_COMM_WORLD,"New ");CHKERRQ(ierr);}
ierr = PetscPrintf(PETSC_COMM_WORLD,"shift: %.14g (s0=",sr->sPres->value);CHKERRQ(ierr);
if (sr->sPres->neighb[0]){ierr = PetscPrintf(PETSC_COMM_WORLD,"%g",sr->sPres->neighb[0]->value);CHKERRQ(ierr);}
ierr = PetscPrintf(PETSC_COMM_WORLD," s1=");CHKERRQ(ierr);
if (sr->sPres->neighb[1]){ierr = PetscPrintf(PETSC_COMM_WORLD,"%g",sr->sPres->neighb[1]->value);CHKERRQ(ierr);}
ierr = PetscPrintf(PETSC_COMM_WORLD,")\n");CHKERRQ(ierr);
}
///////////
 
/* search for deflation */
/* Search for deflation */
ierr = EPSLookForDeflation(eps);CHKERRQ(ierr);
/* krylovSchur */
/* KrylovSchur */
ierr = EPSKrylovSchur_Slice(eps);CHKERRQ(ierr);
ierr = EPSStoreEigenpairs(eps);CHKERRQ(ierr);
/* select new shift */
/* Select new shift */
if(!sr->sPres->comp[1]){
ierr = EPSGetNewShiftValue(eps,1,&newS);CHKERRQ(ierr);
ierr = EPSCreateShift(eps,newS,sr->sPres,sr->sPres->neighb[1]);
}
if(!sr->sPres->comp[0]){
// completing earlier interval
/* Completing earlier interval */
ierr = EPSGetNewShiftValue(eps,0,&newS);CHKERRQ(ierr);
ierr = EPSCreateShift(eps,newS,sr->sPres->neighb[0],sr->sPres);
}
/* preparing for a new search of values */
/* Preparing for a new search of values */
ierr = EPSExtractShift(eps);CHKERRQ(ierr);
}
t2 = MPI_Wtime();
/* checking orthogonality */
if(db>=1){
ierr = STGetOperators(eps->OP,PETSC_NULL,&B);CHKERRQ(ierr);
orthMax=0;
imax=jmax=-1;
ierr = VecDuplicate(sr->V[0],&w);CHKERRQ(ierr);
for(i=0;i < sr->indexEig; i++){
ierr = MatMult(B,sr->V[i],w);CHKERRQ(ierr);
for(j=0;j < sr->indexEig;j++){
if(i != j) {
ierr = VecDot(w,sr->V[j],&inerd);CHKERRQ(ierr);
if(PetscRealPart(inerd)>orthMax){orthMax = PetscRealPart(inerd); imax = i; jmax = j;}
}
}
}
ierr = VecDestroy(&w);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nStored indexEig=%d (of %d)\n (orthog max %g in i=%d j=%d)\n\n",sr->indexEig,sr->numEigs,orthMax,imax,jmax);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," time %g\n",t2-t1);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," number of shifts: %d\n",sr->nShifts);CHKERRQ(ierr);
}
/* updating eps values prior to exit */
//ierr = EPSFreeSolution(eps);
 
/* Updating eps values prior to exit */
ierr = VecDestroyVecs(eps->allocated_ncv,&eps->V);CHKERRQ(ierr);
eps->V = sr->V;
ierr = PetscFree(eps->eigr);CHKERRQ(ierr);
ierr = PetscFree(eps->eigi);CHKERRQ(ierr);
eps->eigr = sr->eig;
eps->eigi = eigi;
eps->its = sr->itsKs;
eps->ncv = eps->allocated_ncv = sr->numEigs;
ierr = PetscFree(eps->errest);CHKERRQ(ierr);
ierr = PetscFree(eps->errest_left);CHKERRQ(ierr);
ierr = PetscMalloc(eps->ncv*sizeof(PetscReal),&eps->errest);CHKERRQ(ierr);
ierr = PetscMalloc(eps->ncv*sizeof(PetscReal),&eps->errest_left);CHKERRQ(ierr);
ierr = PetscFree(eps->perm);CHKERRQ(ierr);
eps->eigr = sr->eig;
eps->eigi = sr->eigi;
eps->errest = sr->errest;
eps->errest_left = errest_left;
eps->perm = sr->perm;
eps->ncv = eps->allocated_ncv = sr->numEigs;
eps->nconv = sr->indexEig;
eps->reason = EPS_CONVERGED_TOL;
eps->its = sr->itsKs;
eps->nds = 0;
eps->DS = PETSC_NULL;
eps->evecsavailable = PETSC_TRUE;
ierr = PetscFree(sr->VDef);CHKERRQ(ierr);
ierr = PetscFree(sr->idxDef);CHKERRQ(ierr);
ierr = PetscFree(sr->pending);CHKERRQ(ierr);
// reviewing list of shifts to free memmory
ierr = PetscFree(sr->monit);CHKERRQ(ierr);
/* Reviewing list of shifts to free memory */
shift s = sr->s0;
if(s){
while(s->neighb[1]){
785,4 → 661,3
ierr = PetscFree(sr);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
/branches/slepc-3_2-proj/src/eps/impls/krylov/krylovschur/krylovschur.c
50,6 → 50,7
PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
{
PetscErrorCode ierr;
PetscBool issinv;
 
PetscFunctionBegin;
/* spectrum slicing requires special treatment of default values */
60,15 → 61,16
if (!((PetscObject)(eps->OP))->type_name) { /* default to shift-and-invert */
ierr = STSetType(eps->OP,STSINVERT);CHKERRQ(ierr);
}
if(eps->intb >= PETSC_MAX_REAL){/* right-open interval */
if(eps->inta <= PETSC_MIN_REAL) SETERRQ(((PetscObject)eps)->comm,1,"The defined computational interval should have at least one of their sides bounded");
ierr = PetscTypeCompare((PetscObject)eps->OP,STSINVERT,&issinv);CHKERRQ(ierr);
if (!issinv) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Shift-and-invert ST is needed for spectrum slicing");
if (eps->intb >= PETSC_MAX_REAL) { /* right-open interval */
if (eps->inta <= PETSC_MIN_REAL) SETERRQ(((PetscObject)eps)->comm,1,"The defined computational interval should have at least one of their sides bounded");
ierr = STSetDefaultShift(eps->OP,eps->inta);CHKERRQ(ierr);
}else ierr = STSetDefaultShift(eps->OP,eps->intb);CHKERRQ(ierr);
}
else { ierr = STSetDefaultShift(eps->OP,eps->intb);CHKERRQ(ierr); }
 
if (eps->nev==1) eps->nev = 20; /* nev not set, use default value */
if (eps->nev==1) eps->nev = 40; /* nev not set, use default value */
if (eps->nev<10) SETERRQ(((PetscObject)eps)->comm,1,"nev cannot be less than 10 in spectrum slicing runs");
eps->mpd = 2*eps->nev;
eps->ncv = 2*eps->nev;
eps->ops->backtransform = PETSC_NULL;
}
 
191,7 → 193,7
ierr = EPSProjectedKSNonsym(eps,l,S,eps->ncv,Q,nv);CHKERRQ(ierr);
 
/* Check convergence */
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,S,eps->ncv,Q,eps->V,nv,beta,1.0,&k,work);CHKERRQ(ierr);
ierr = EPSKrylovConvergence(eps,PETSC_FALSE,eps->trackall,eps->nconv,nv-eps->nconv,S,eps->ncv,Q,eps->V,nv,beta,1.0,&k,work);CHKERRQ(ierr);
if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
if (k >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
/branches/slepc-3_2-proj/src/svd/interface/svdsetup.c
111,7 → 111,6
PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
if (svd->setupcalled) PetscFunctionReturn(0);
ierr = PetscLogEventBegin(SVD_SetUp,svd,0,0,0);CHKERRQ(ierr);
if (!svd->ip) { ierr = SVDGetIP(svd,&svd->ip);CHKERRQ(ierr); }
 
/* Set default solver type (SVDSetFromOptions was not called) */
if (!((PetscObject)svd)->type_name) {
/branches/slepc-3_2-proj/bin/matlab/classes/Slepc.m
8,8 → 8,30
% See also: help Petsc
%
% Each of the following have their own help
% PetscEPS - Eigenvalue solver class
% SlepcEPS - Eigenvalue solver class
% SlepcST - Spectral transformation class
% SlepcSVD - Singular value solver class
% SlepcQEP - Quadratic eigenvalue solver class
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
properties (Constant)
end
end
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcFinalize.m
1,5 → 1,24
function err = SlepcFinalize()
% SlepcFinalize: clean up SLEPc classes in MATLAB.
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
err = calllib('libslepc', 'SlepcFinalize');PetscCHKERRQ(err);
 
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcInitialize.m
1,5 → 1,5
function err = SlepcInitialize(args,argfile,arghelp)
%
% SlepcInitialize: start to use SLEPc classes in MATLAB.
% In order to use the SLEPc MATLAB classes, PETSc must have been configured with
% specific options. See: help PetscInitialize.
%
7,7 → 7,25
% ${SLEPC_DIR}/bin/matlab/classes
%
% In MATLAB use help Slepc to get started using SLEPc from MATLAB
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
if exist('PetscInitialize')~=2
error('Must add ${PETSC_DIR}/bin/matlab/classes to your MATLAB path')
/branches/slepc-3_2-proj/bin/matlab/classes/slepcmatlabheader.h
1,9 → 1,28
 
/*
This is used by bin/matlab/classes/SlepcInitialize() to define to Matlab all the functions available in the
SLEPc shared library. We cannot simply use the regular SLEPc include files because they are too complicated for
Matlab to parse.
This is used by bin/matlab/classes/SlepcInitialize() to define to Matlab all
the functions available in the SLEPc shared library. We cannot simply use
the regular SLEPc include files because they are too complicated for Matlab
to parse.
 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
int SlepcInitializeMatlab(int,char **,const char*,const char*);
int SlepcInitializedMatlab(void);
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcQEP.m
8,7 → 8,26
% QEP.SetOperators(M,C,K);
% (optional) QEP.SetProblemType(...);
% QEP.SetFromOptions();
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
properties (Constant)
GENERAL=1;
HERMITIAN=2;
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcST.m
4,7 → 4,26
%
% Creation from an EPS:
% st = eps.GetST();
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
properties (Constant)
HEP=1;
GHEP=2;
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcEPS.m
8,7 → 8,26
% eps.SetOperators(A,B);
% (optional) eps.SetProblemType(...);
% eps.SetFromOptions;
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
properties (Constant)
HEP=1;
GHEP=2;
/branches/slepc-3_2-proj/bin/matlab/classes/SlepcSVD.m
7,7 → 7,26
% svd.SetType('cross');
% svd.SetOperator(A);
% svd.SetFromOptions;
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
properties (Constant)
SVD_TRANSPOSE_EXPLICIT=0;
SVD_TRANSPOSE_IMPLICIT=1;
/branches/slepc-3_2-proj/bin/matlab/classes/examples/tutorials/exEPS.m
4,6 → 4,25
% User creates directly a PETSc Mat
%
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
%%
% Set the Matlab path and initialize SLEPc
%
33,7 → 52,6
end
mat.AssemblyBegin(PetscMat.FINAL_ASSEMBLY);
mat.AssemblyEnd(PetscMat.FINAL_ASSEMBLY);
%mat.View;
 
%%
% Create the eigensolver, pass the matrix and solve the problem
56,7 → 74,6
fprintf(' %12f%+12f %12g\n',real(lambda),imag(lambda),relerr)
end
end
%eps.View();
 
%%
% Free objects and shutdown SLEPc
/branches/slepc-3_2-proj/bin/matlab/classes/examples/tutorials/exSVD.m
4,6 → 4,25
% User creates directly a PETSc Mat
%
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
%%
% Set the Matlab path and initialize SLEPc
%
33,7 → 52,6
end
mat.AssemblyBegin(PetscMat.FINAL_ASSEMBLY);
mat.AssemblyEnd(PetscMat.FINAL_ASSEMBLY);
%mat.View;
 
%%
% Create the solver, pass the matrix and solve the problem
51,7 → 69,6
relerr = svd.ComputeRelativeError(i);
fprintf(' %12f %12g\n',sigma,relerr)
end
%svd.View();
 
%%
% Free objects and shutdown SLEPc
/branches/slepc-3_2-proj/bin/matlab/classes/examples/tutorials/exQEP.m
4,6 → 4,25
% User creates directly the three PETSc Mat
%
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
%%
% Set the Matlab path and initialize SLEPc
%
73,7 → 92,6
fprintf(' %9f%+9f %12g\n',real(lambda),imag(lambda),relerr)
end
end
%qep.View();
 
%%
% Free objects and shutdown SLEPc
/branches/slepc-3_2-proj/bin/matlab/classes/examples/tutorials/exST.m
4,6 → 4,25
% with the matrices loaded from file
%
 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
%%
% Set the Matlab path and initialize SLEPc
%
21,11 → 40,11
% Load matrices from file
%
SLEPC_DIR = getenv('SLEPC_DIR');
viewer = PetscViewer([SLEPC_DIR '/src/mat/examples/bfw62a.petsc'],Petsc.FILE_MODE_READ);
viewer = PetscViewer([SLEPC_DIR '/share/slepc/datafiles/matrices/bfw62a.petsc'],Petsc.FILE_MODE_READ);
A = PetscMat();
A.Load(viewer);
viewer.Destroy;
viewer = PetscViewer([SLEPC_DIR '/src/mat/examples/bfw62b.petsc'],Petsc.FILE_MODE_READ);
viewer = PetscViewer([SLEPC_DIR '/share/slepc/datafiles/matrices/bfw62b.petsc'],Petsc.FILE_MODE_READ);
B = PetscMat();
B.Load(viewer);
viewer.Destroy;
/branches/slepc-3_2-proj/docs/manual/exEPS.m New file
0,0 → 1,0
link ../../bin/matlab/classes/examples/tutorials/exEPS.m
slepc-3_2-proj/docs/manual/exEPS.m Property changes : Added: svn:special
+ *
/branches/slepc-3_2-proj/docs/manual/extra.tex
40,7 → 40,7
 
Methods implemented in \ident{EPS} merely require vector operations and matrix-vector products. In \petsc, mathematical objects such as vectors and matrices have an interface that is independent of the underlying data structures. \slepc manipulates vectors and matrices via this interface and, therefore, it can be used with any of the matrix representations provided by \petsc, including dense, sparse, and symmetric formats, either sequential or parallel.
 
The above statement must be reconsidered when using \ident{EPS} in combination with \ident{ST}. As explained in chapter \ref{cap:st}, in many cases the operator associated to a spectral transformation not only consists in pure matrix-vector products but also other operations may be required as well, most notably a linear system solve (see Table \ref{tab:op}). In this case, the limitation is that there must be support for the requested operation for the selected matrix representation. For instance, if one wants to use \texttt{cholesky} for the solution of the linear systems, then it may be necessary to work with a symmetric matrix format such as \texttt{MATSEQSBAIJ}.
The above statement must be reconsidered when using \ident{EPS} in combination with \ident{ST}. As explained in chapter \ref{cap:st}, in many cases the operator associated to a spectral transformation not only consists in pure matrix-vector products but also other operations may be required as well, most notably a linear system solve (see Table \ref{tab:op}). In this case, the limitation is that there must be support for the requested operation for the selected matrix representation. For instance, if one wants to use \texttt{cholesky} for the solution of the linear systems, then it may be necessary to work with a symmetric matrix format such as \texttt{MATSBAIJ}.
 
\paragraph{Shell Matrices.}
 
57,7 → 57,7
 
All these cases can be easily handled in \slepc by means of shell matrices. These are matrices that do not require explicit storage of the component values. Instead, the user must provide subroutines for all the necessary matrix operations, typically only the application of the linear operator to a vector.
 
Shell matrices, also called matrix-free matrices, are created in \petsc{} with the command \ident{MatCreateShell}. Then, the function \ident{MatShellSetOperation} is used to provide any user-defined shell matrix operations (see the \petsc{} documentation for additional details). Several examples are available in \slepc that illustrate how to solve a matrix-free eigenvalue problem.
Shell matrices, also called matrix-free matrices, are created in \petsc with the command \ident{MatCreateShell}. Then, the function \ident{MatShellSetOperation} is used to provide any user-defined shell matrix operations (see the \petsc{} documentation for additional details). Several examples are available in \slepc that illustrate how to solve a matrix-free eigenvalue problem.
 
In the simplest case, defining matrix-vector product operations (\ident{MATOP\_MULT}) is enough for using \ident{EPS} with shell matrices. However, in the case of generalized problems, if matrix $B$ is also a shell matrix then it may be necessary to define other operations in order to be able to solve the linear system successfully, for example \ident{MATOP\_GET\_DIAGONAL} to use an iterative linear solver with Jacobi preconditioning. On the other hand, if the shift-and-invert \ident{ST} is to be used, then in addition it may also be necessary to define \ident{MATOP\_SHIFT} or \ident{MATOP\_AXPY} (see \S\ref{sec:explicit} for discussion).
 
67,16 → 67,27
\section{GPU Computing}
\label{sec:gpu}
 
Experimental support for graphics processor unit (GPU) computing is included in \slepc 3.2. This is related to the previous section because GPU support in \petsc is based on using special types of \texttt{Mat} and \texttt{Vec}. Currently, GPU support in \slepc has been tested only in the Krylov-Schur \ident{EPS} solver, although it may work in other solvers as well. Regarding \petsc, all iterative linear solvers are prepared to run on the GPU, but this is not the case for direct solvers and preconditioners (see \petsc documentation for details). The user must not expect a spectacular performance boost, but in general moderate gains can be achieved by running the eigensolver on the GPU instead of the CPU (in some cases a 10-fold improvement).
 
PETSc's GPU support currently relies on NVIDIA CUDA Toolkit 4.0\footnote{\url{http://developer.nvidia.com/cuda-toolkit-40}}, that provides a C/C++ compiler with CUDA extensions and the Thrust library, together with the Cusp library\footnote{\url{http://code.google.com/p/cusp-library}}, that implements sparse linear algebra and graph computations on CUDA based on Thrust data structures. For instance, to configure \petsc with GPU support in single precision arithmetic use the following options:
\begin{Verbatim}[fontsize=\small]
$ ./configure --with-precision=single --with-cuda --with-thrust
--with-cusp --with-cusp-dir=/path/to/cusp-library
\end{Verbatim}
 
\ident{VECCUSP} and \ident{MATAIJCUSP} are currently the mechanism in \petsc to run a computation on the GPU. \ident{VECCUSP} is a special type of \texttt{Vec} whose array is mirrored in the GPU (and similarly for \ident{MATAIJCUSP}). \petsc takes care of keeping memory coherence between the two copies of the array, and performs the computation on the GPU when possible, trying to avoid unnecessary copies between the host and the device. For maximum efficiency, the user has to make sure that all vectors and matrices are of these types. If they are created in the standard way (\texttt{VecCreate} plus \texttt{VecSetFromOptions}) then it is sufficient to run the program with
\begin{Verbatim}[fontsize=\small]
$ ./program -vec_type cusp -mat_type aijcusp
\end{Verbatim}
 
%---------------------------------------------------
\section{Extending SLEPc}
\label{sec:shell}
 
Shell matrices are a simple mechanism of extensibility, in the sense that the package is extended with new user-defined matrix objects. Once the new matrix has been defined, it can be used by \slepc in the same way as the rest of the matrices as long as the required operations are provided.
Shell matrices, presented in \S\ref{sec:supported}, are a simple mechanism of extensibility, in the sense that the package is extended with new user-defined matrix objects. Once the new matrix has been defined, it can be used by \slepc in the same way as the rest of the matrices as long as the required operations are provided.
 
A similar mechanism is available in \slepc also for extending the system incorporating new spectral transformations (\ident{ST}). This is done by using the \ident{STSHELL} spectral transformation type, in a similar way as shell matrices or shell preconditioners. In this case, the user defines how the operator is applied to a vector and optionally how the computed eigenvalues are transformed back to the solution of the original problem (see \S\ref{sec:shell} for details). This tool is intended for simple spectral transformations. For more sophisticated transformations, the user should register a new \ident{ST} type (see below).
 
At least, user-defined spectral transformations have to define how the operator is to be applied to a vector. Optionally, they can also specify the way in which computed eigenvalues must be transformed back to the solution of the original eigenproblem. An example program is provided in the \slepc distribution in order to illustrate the use of shell transformations.
 
The function
\findex{STShellSetApply}
\begin{Verbatim}[fontsize=\small]
87,20 → 98,20
\begin{Verbatim}[fontsize=\small]
STShellSetBackTransform(ST,PetscErrorCode(*)(ST,PetscInt,PetscScalar*,PetscScalar*));
\end{Verbatim}
can be used optionally to specify the routine for the back-transformation of eigenvalues. The two functions provided by the user can make use of any required user-defined information via a context that can be retrieved with \ident{STShellGetContext}.
can be used optionally to specify the routine for the back-transformation of eigenvalues. The two functions provided by the user can make use of any required user-defined information via a context that can be retrieved with \ident{STShellGetContext}. An example program is provided in the \slepc distribution in order to illustrate the use of shell transformations.
 
\slepc further supports extensibility by allowing application programmers to code their own subroutines for unimplemented features such as new eigensolvers or new spectral transformations. It is possible to register these new methods to the system and use them as the rest of standard subroutines. For example, to implement a variant of the Subspace Iteration method, one could copy the \slepc code associated to the \texttt{subspace} solver, modify it and register a new \ident{EPS} type with the following line of code
\begin{Verbatim}[fontsize=\small]
EPSRegister("newsubspace",0,"EPSCreate_NEWSUB",EPSCreate_NEWSUB);
\end{Verbatim}
After this call, the new solver could be used in the same way as the rest of \slepc solvers. For instance,
After this call, the new solver could be used in the same way as the rest of \slepc solvers:
\begin{Verbatim}[fontsize=\small]
$ ./program -eps_type newsubspace
\end{Verbatim}
 
\ident{EPSRegister} can be used to register new types whose code is linked into the executable. To register types in a dynamic library use \ident{EPSRegisterDynamic}.
 
A similar mechanism is available for registering new types of classes \ident{ST} and \ident{SVD}.
A similar mechanism is available for registering new types of classes \ident{ST}, \ident{SVD} and \ident{QEP}.
 
%---------------------------------------------------
\section{Directory Structure}
164,7 → 175,7
\item Use the runtime option \Verb!-eps_type <type>! to select the solver.
\end{enumerate}
 
Exceptions to the above rule are \lapack and \blopex, which should be enabled during \petsc's configuration.
Exceptions to the above rule are \lapack, which should be enabled during \petsc's configuration, and \blopex, that must be installed with \Verb!--download-blopex! in \slepc's configure.
 
\subsection*{\underline{\lapack}}
\begin{description}
177,7 → 188,7
\slepc explicitly creates the operator matrix in dense form and then the appropriate \lapack driver routine is invoked. Therefore, this interface should be used only for testing and validation purposes and not in a production code. The operator matrix is created by applying the operator to the columns of the identity matrix.
 
\item[Installation.]
The \slepc interface to \lapack can be used directly. If \slepc's configure script complains about missing \lapack functions, then configure \petsc with option \texttt{-{}-download-c-blas-lapack}.
The \slepc interface to \lapack can be used directly. If \slepc's configure script complains about missing \lapack functions, then configure \petsc with option \texttt{-{}-download-f2cblaslapack}.
\end{description}
 
\subsection*{\underline{\arpack}}
221,8 → 232,6
\blzpack\ uses a combination of partial and selective re-orthogonalization strategies. It can be run in either sequential or parallel mode, by means of MPI or PVM interfaces, and it uses the reverse communication strategy.
\item[Installation.] For the compilation of the \texttt{libblzpack.a} library, first check the appropriate architecture file in the directory \texttt{sys/MACROS} and then type \texttt{creator -mpi}.
\item[Specific options.] The \slepc interface to this package allows the user to specify the block size with the function \ident{EPSBlzpackSetBlockSize} or at run time with the option \Verb!-eps_blzpack_! \Verb!block_size <size>!. Also, the function \ident{EPSBlzpackSetNSteps} can be used to set the maximum number of steps per run (also with \Verb!-eps_blzpack_nsteps!).
 
For the spectrum slicing feature, \slepc allows the programmer to provide the computational interval with the option \Verb!-eps_blzpack_interval!, or with the function \ident{EPSBlzpackSetInterval} in the program source.
\end{description}
 
\subsection*{\underline{\trlan}}
230,10 → 239,11
\setlength{\itemsep}{0pt}
\item[References.]\citep{Wu:2000:TLM}.
\item[Website.] \url{http://crd.lbl.gov/\~kewu/trlan.html}.
\item[Version.] 201009.
\item[Summary.] This package provides a Fortran 90 implementation of the dynamic thick-restart Lanczos algorithm. This is a specialized version of Lanczos that targets only the case in which one wants both eigenvalues and eigenvectors of a large real symmetric eigenvalue problem that cannot use the shift-and-invert scheme. In this case the standard non-restarted Lanczos algorithm requires to store a large number of Lanczos vectors, what can cause storage problems and make each iteration of the method very expensive.
 
\trlan{} requires the user to provide a matrix-vector multiplication routine. The parallel version uses MPI as the message passing layer.
\item[Installation.] To install this package, it is necessary to have access to a Fortran 90 compiler. The compiler name and the options used are specified in the file called \texttt{Make.inc}. To generate the library, type \texttt{make libtrlan\_mpi.a} in the \texttt{TRLan} directory.
\item[Installation.] To install this package, it is necessary to have access to a Fortran 90 compiler. The compiler name and the options used are specified in the file called \texttt{Make.inc}. To generate the library, type \texttt{make plib} in the \texttt{TRLan} directory.
\end{description}
 
\subsection*{\underline{\blopex}}
242,7 → 252,7
\item[References.]\citep{Knyazev:2007:BLO}.
\item[Website.] \url{http://code.google.com/p/blopex}.
\item[Summary.] \blopex is a package that implements the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) method for computing several extreme eigenpairs of symmetric positive generalized eigenproblems. Numerical comparisons suggest that this method is a genuine analog for eigenproblems of the standard preconditioned conjugate gradient method for symmetric linear systems.
\item[Installation.] In order to use \blopex from \slepc, it is a prior requirement that \petsc has been configured with \Verb!./configure --download-hypre --download-blopex!.
\item[Installation.] In order to use \blopex from \slepc, it necessary to install it during \slepc's configuration: \Verb!./configure --download-blopex!.
\item[Specific options.] Since BLOPEX contains preconditioned solvers, the \slepc interface uses \ident{STPRECOND}, as described in \ref{sec:precond}.
\end{description}
 
259,3 → 269,13
\section{Matlab Interface}
\label{sec:matlab}
 
\slepc 3.2 includes an interface intended to make most of \slepc's functionality available from Matlab. It is experimental and needs further development, so users planning to use it seriously are recommended to contact the authors. Below are some guidelines for using this interface.
 
First of all, \petsc must have been configured with the Matlab interface enabled. This can be done as follows (check \petsc documentation for details):
\begin{Verbatim}[fontsize=\small]
$ ./configure --with-matlab --with-matlab-engine --with-shared-libraries
\end{Verbatim}
 
Once the \petsc and \slepc libraries have been built, one has to set Matlab's path to point to the directories containing Matlab classes: \Verb!$SLEPC_DIR/bin/matlab/classes! and \Verb!$PETSC_DIR/bin/matlab/classes!. Below we show a simple Matlab example (included in \slepc's distribution) that does this, and then solves a simple eigenproblem.
\MyVerbatimInput{exEPS.m}
 
/branches/slepc-3_2-proj/docs/manual/qep.tex
10,7 → 10,7
 
\section{\label{sec:qep}Overview of Quadratic Eigenproblems}
 
In this section, we review some basic properties of quadratic eigenvalue problems. The main goal is to set up the notation as well as describe the linearization approaches that will be employed for solving via the \ident{EPS} object. For additional background material about the quadratic eigenproblem, the reader is referred to \citep{Tisseur:2001:QEP}. As always, some details of the implemented methods can be found in the \slepc \hyperlink{str}{technical reports}.
In this section, we review some basic properties of quadratic eigenvalue problems. The main goal is to set up the notation as well as to describe the linearization approaches that will be employed for solving via the \ident{EPS} object. For additional background material about the quadratic eigenproblem, the reader is referred to \citep{Tisseur:2001:QEP}. As always, some details of the implemented methods can be found in the \slepc \hyperlink{str}{technical reports}.
 
In many applications, e.g., problems arising from second-order differential equations such as the analysis of damped vibrating systems, the eigenproblem to be solved is quadratic,
\begin{equation}
20,7 → 20,7
 
It is important to point out some outstanding differences with respect to the linear eigenproblem. In the quadratic eigenproblem, the number of eigenvalues is $2n$, and the corresponding eigenvectors do not form a linearly independent set. If $M$ is singular, some eigenvalues are infinite. Even when the three matrices are symmetric and positive definite, there is no guarantee that the eigenvalues are real, but still methods can exploit symmetry to some extent. Furthermore, numerical difficulties are more likely than in the linear case, so the computed solution can sometimes be untrustworthy.
 
If Eq.\ \ref{eq:eigquad} is written as $Q(\lambda)x=0$, where $Q$ is the matrix polynomial, then multiplication by $\lambda^{-2}$ results in $R(\lambda^{-1})x=0$, where $R$ is a matrix polynomial with the coefficients in the reverse order. In other words, if a method is available for computing the largest eigenvalues, then reversing the roles of $M$ and $K$ results in the computation of the smallest eigenvalues. In general, it is also possible to formulate different spectral trasformation for computing eigenvalues closest to a given target.
If Eq.\ \ref{eq:eigquad} is written as $Q(\lambda)x=0$, where $Q$ is the matrix polynomial, then multiplication by $\lambda^{-2}$ results in $R(\lambda^{-1})x=0$, where $R$ is a matrix polynomial with the coefficients in the reverse order. In other words, if a method is available for computing the largest eigenvalues, then reversing the roles of $M$ and $K$ results in the computation of the smallest eigenvalues. In general, it is also possible to formulate different spectral transformations for computing eigenvalues closest to a given target.
 
\paragraph{Problem Types.}
 
/branches/slepc-3_2-proj/docs/manual/st.tex
153,7 → 153,7
 
\subsection{Shift-and-invert}
 
\begin{figure}
\begin{figure}[t]
\centering
% Requires gnuplot
\begin{tikzpicture}
434,17 → 434,21
The \slepc interface hides all the complications of the algorithm. However, the user must be aware of all the restrictions for this technique to be employed:
\begin{itemize}
\item This is currently implemented only in Krylov-Schur, for \ident{EPS\_GHEP} problems.
\item The implementation is based on shift-and-invert, so \ident{STSINVERT} must be selected. Furthermore, direct linear solvers must be used.
\item The direct linear solver must provide the matrix inertia (see \petsc's \ident{MatGetInertia}). In particular, a symmetric factorization must be used (\texttt{cholesky}), which requires a symmetric matrix storage (\texttt{sbaij}).
\item The implementation is based on shift-and-invert, so \ident{STSINVERT} must be used (it is selected by default). Furthermore, direct linear solvers must be used.
\item The direct linear solver must provide the matrix inertia (see \petsc's \ident{MatGetInertia}). In particular, a symmetric factorization must be used (\texttt{cholesky}).
\end{itemize}
 
An example command-line that sets up all the required options is:
\begin{Verbatim}[fontsize=\small]
$ ./ex13 -eps_interval 0.4,0.8 -mat_type sbaij
-st_ksp_type preonly -st_pc_type cholesky
$ ./ex13 -eps_interval 0.4,0.8 -st_ksp_type preonly -st_pc_type cholesky
\end{Verbatim}
 
Note that \petsc's Cholesky factorization is not parallel, so for doing spectrum slicing in parallel it is required to use an external solver that supports inertia, e.g., MUMPS (see \S\ref{sec:lin} on how to use external linear solvers).
Note that \petsc's Cholesky factorization is not parallel, so for doing spectrum slicing in parallel it is required to use an external solver that supports inertia, e.g., MUMPS (see \S\ref{sec:lin} on how to use external linear solvers):
\begin{Verbatim}[fontsize=\small]
$ ./ex13 -eps_interval 0.4,0.8 -st_ksp_type preonly -st_pc_type cholesky
-st_pc_factor_mat_solver_package mumps -mat_mumps_icntl_13 1
\end{Verbatim}
The last option is required by MUMPS to compute the inertia.
 
Apart from the above recommendations, the following must be taken into account:
\begin{itemize}
slepc-3_2-proj Property changes : Modified: svn:mergeinfo
Merged /trunk:r2579-2602