Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 827 → Rev 828

/trunk/src/st/interface/stsolve.c
196,7 → 196,7
PetscErrorCode ierr;
Vec in,out;
int i,M,m,*rows,start,end;
PetscScalar *array,zero = 0.0,one = 1.0;
PetscScalar *array,one = 1.0;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(st,ST_COOKIE,1);
212,7 → 212,7
ierr = MatCreateMPIDense(st->comm,m,m,M,M,PETSC_NULL,mat);CHKERRQ(ierr);
 
for (i=0; i<M; i++) {
ierr = VecSet(&zero,in);CHKERRQ(ierr);
ierr = VecSet(in,0.0);CHKERRQ(ierr);
ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
ierr = VecAssemblyEnd(in);CHKERRQ(ierr);
/trunk/src/st/interface/shellmat.c
12,21 → 12,18
{
PetscErrorCode ierr;
ST ctx;
PetscScalar alpha;
 
PetscFunctionBegin;
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
 
ierr = MatMult(ctx->A,x,y);CHKERRQ(ierr);
if (ctx->sigma != 0.0) {
alpha = -ctx->sigma;
if (ctx->B) { /* y = (A - sB) x */
ierr = MatMult(ctx->B,x,ctx->w);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,ctx->w,y);CHKERRQ(ierr);
ierr = VecAXPY(y,-ctx->sigma,ctx->w);CHKERRQ(ierr);
} else { /* y = (A - sI) x */
ierr = VecAXPY(y,-ctx->sigma,x);CHKERRQ(ierr);
}
else { /* y = (A - sI) x */
ierr = VecAXPY(&alpha,x,y);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
}
37,21 → 34,18
{
PetscErrorCode ierr;
ST ctx;
PetscScalar alpha;
 
PetscFunctionBegin;
ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
 
ierr = MatMultTranspose(ctx->A,x,y);CHKERRQ(ierr);
if (ctx->sigma != 0.0) {
alpha = -ctx->sigma;
if (ctx->B) { /* y = (A - sB) x */
ierr = MatMultTranspose(ctx->B,x,ctx->w);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,ctx->w,y);CHKERRQ(ierr);
ierr = VecAXPY(y,-ctx->sigma,ctx->w);CHKERRQ(ierr);
} else { /* y = (A - sI) x */
ierr = VecAXPY(y,-ctx->sigma,x);CHKERRQ(ierr);
}
else { /* y = (A - sI) x */
ierr = VecAXPY(&alpha,x,y);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
}
62,7 → 56,6
{
PetscErrorCode ierr;
ST ctx;
PetscScalar alpha;
Vec diagb;
 
PetscFunctionBegin;
70,16 → 63,14
 
ierr = MatGetDiagonal(ctx->A,diag);CHKERRQ(ierr);
if (ctx->sigma != 0.0) {
alpha = -ctx->sigma;
if (ctx->B) {
ierr = VecDuplicate(diag,&diagb);CHKERRQ(ierr);
ierr = MatGetDiagonal(ctx->B,diagb);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,diagb,diag);CHKERRQ(ierr);
ierr = VecAXPY(diag,-ctx->sigma,diagb);CHKERRQ(ierr);
ierr = VecDestroy(diagb);CHKERRQ(ierr);
} else {
ierr = VecShift(diag,-ctx->sigma);CHKERRQ(ierr);
}
else {
ierr = VecShift(&alpha,diag);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
}
/trunk/src/st/impls/fold/fold.c
14,36 → 14,30
PetscErrorCode STApply_Fold(ST st,Vec x,Vec y)
{
PetscErrorCode ierr;
PetscScalar alpha;
ST_FOLD *ctx = (ST_FOLD *) st->data;
 
PetscFunctionBegin;
if (st->B) {
/* generalized eigenproblem: y = (B^-1 A + sI)^2 x */
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,ctx->w2);CHKERRQ(ierr);
if (st->sigma != 0.0) {
alpha = - st->sigma;
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,ctx->w2);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,x,ctx->w2);CHKERRQ(ierr);
ierr = MatMult(st->A,ctx->w2,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,ctx->w2,y);CHKERRQ(ierr);
} else {
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
ierr = MatMult(st->A,y,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(ctx->w2,-st->sigma,x);CHKERRQ(ierr);
}
ierr = MatMult(st->A,ctx->w2,st->w);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
if (st->sigma != 0.0) {
ierr = VecAXPY(y,-st->sigma,ctx->w2);CHKERRQ(ierr);
}
} else {
/* standard eigenproblem: y = (A + sI)^2 x */
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
alpha = - st->sigma;
if (st->sigma != 0.0) {
ierr = VecAXPY(&alpha,x,st->w);CHKERRQ(ierr);
ierr = VecAXPY(st->w,-st->sigma,x);CHKERRQ(ierr);
}
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
if (st->sigma != 0.0) {
ierr = VecAXPY(&alpha,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(y,-st->sigma,st->w);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
54,36 → 48,30
PetscErrorCode STApplyTranspose_Fold(ST st,Vec x,Vec y)
{
PetscErrorCode ierr;
PetscScalar alpha;
ST_FOLD *ctx = (ST_FOLD *) st->data;
 
PetscFunctionBegin;
if (st->B) {
/* generalized eigenproblem: y = (A^T B^-T + sI)^2 x */
ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,ctx->w2);CHKERRQ(ierr);
if (st->sigma != 0.0) {
alpha = - st->sigma;
ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,ctx->w2);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,x,ctx->w2);CHKERRQ(ierr);
ierr = STAssociatedKSPSolveTranspose(st,ctx->w2,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(&alpha,ctx->w2,y);CHKERRQ(ierr);
} else {
ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
ierr = STAssociatedKSPSolveTranspose(st,y,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(ctx->w2,-st->sigma,x);CHKERRQ(ierr);
}
ierr = STAssociatedKSPSolveTranspose(st,ctx->w2,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
if (st->sigma != 0.0) {
ierr = VecAXPY(y,-st->sigma,ctx->w2);CHKERRQ(ierr);
}
} else {
/* standard eigenproblem: y = (A^T + sI)^2 x */
ierr = MatMultTranspose(st->A,x,st->w);CHKERRQ(ierr);
alpha = - st->sigma;
if (st->sigma != 0.0) {
ierr = VecAXPY(&alpha,x,st->w);CHKERRQ(ierr);
ierr = VecAXPY(st->w,-st->sigma,x);CHKERRQ(ierr);
}
ierr = MatMultTranspose(st->A,st->w,y);CHKERRQ(ierr);
if (st->sigma != 0.0) {
ierr = VecAXPY(&alpha,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(y,-st->sigma,st->w);CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
/trunk/src/st/impls/cayley/cayley.c
24,13 → 24,13
/* generalized eigenproblem: y = (A - sB)^-1 (A + tB)x */
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->B,x,ctx->w2);CHKERRQ(ierr);
ierr = VecAXPY(&tau,ctx->w2,st->w);CHKERRQ(ierr);
ierr = VecAXPY(st->w,tau,ctx->w2);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
}
else {
/* standard eigenproblem: y = (A - sI)^-1 (A + tI)x */
ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
ierr = VecAXPY(&tau,x,st->w);CHKERRQ(ierr);
ierr = VecAXPY(st->w,tau,x);CHKERRQ(ierr);
ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
52,13 → 52,13
ierr = STAssociatedKSPSolve(st,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
ierr = MatMult(st->B,st->w,ctx->w2);CHKERRQ(ierr);
ierr = VecAXPY(&tau,ctx->w2,y);CHKERRQ(ierr);
ierr = VecAXPY(y,tau,ctx->w2);CHKERRQ(ierr);
}
else {
/* standard eigenproblem: y = (A + tI)^T (A - sI)^-T x */
ierr = STAssociatedKSPSolve(st,x,st->w);CHKERRQ(ierr);
ierr = MatMult(st->A,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(&tau,st->w,y);CHKERRQ(ierr);
ierr = VecAXPY(y,tau,st->w);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
89,12 → 89,12
/* generalized eigenproblem: y = (A + tB)x */
ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
ierr = MatMult(st->B,x,ctx->w2);CHKERRQ(ierr);
ierr = VecAXPY(&tau,ctx->w2,y);CHKERRQ(ierr);
ierr = VecAXPY(y,tau,ctx->w2);CHKERRQ(ierr);
}
else {
/* standard eigenproblem: y = (A + tI)x */
ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
ierr = VecAXPY(&tau,x,y);CHKERRQ(ierr);
ierr = VecAXPY(y,tau,x);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
132,13 → 132,14
PetscErrorCode STPost_Cayley(ST st)
{
PetscErrorCode ierr;
PetscScalar alpha;
 
PetscFunctionBegin;
if (st->shift_matrix == STMATMODE_INPLACE) {
alpha = st->sigma;
if( st->B ) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha, st->A); CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,st->sigma); CHKERRQ(ierr);
}
st->setupcalled = 0;
}
PetscFunctionReturn(0);
150,7 → 151,6
{
PetscErrorCode ierr;
ST_CAYLEY *ctx = (ST_CAYLEY *) st->data;
PetscScalar alpha;
 
PetscFunctionBegin;
 
164,11 → 164,10
case STMATMODE_INPLACE:
st->mat = PETSC_NULL;
if (st->sigma != 0.0) {
alpha = -st->sigma;
if (st->B) {
ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr);
ierr = MatAXPY(st->A,-st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(&alpha,st->A);CHKERRQ(ierr);
ierr = MatShift(st->A,-st->sigma);CHKERRQ(ierr);
}
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
182,11 → 181,10
default:
ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
if (st->sigma != 0.0) {
alpha = -st->sigma;
if (st->B) {
ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr);
ierr = MatAXPY(st->mat,-st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr);
ierr = MatShift(st->mat,-st->sigma);CHKERRQ(ierr);
}
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
206,7 → 204,6
PetscErrorCode STSetShift_Cayley(ST st,PetscScalar newshift)
{
PetscErrorCode ierr;
PetscScalar alpha;
ST_CAYLEY *ctx = (ST_CAYLEY *) st->data;
 
PetscFunctionBegin;
218,15 → 215,19
case STMATMODE_INPLACE:
/* Undo previous operations */
if (st->sigma != 0.0) {
alpha = st->sigma;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,st->sigma);CHKERRQ(ierr);
}
}
/* Apply new shift */
if (newshift != 0.0) {
alpha = -newshift;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->A,-newshift,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,-newshift);CHKERRQ(ierr);
}
}
ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
break;
236,9 → 237,8
default:
ierr = MatCopy(st->A, st->mat,SUBSET_NONZERO_PATTERN); CHKERRQ(ierr);
if (newshift != 0.0) {
alpha = -newshift;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr); }
if (st->B) { ierr = MatAXPY(st->mat,-newshift,st->B,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(st->mat,-newshift);CHKERRQ(ierr); }
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
* improve performance when solving a number of related eigenproblems */
/trunk/src/st/impls/sinvert/sinvert.c
96,13 → 96,14
PetscErrorCode STPost_Sinvert(ST st)
{
PetscErrorCode ierr;
PetscScalar alpha;
 
PetscFunctionBegin;
if (st->shift_matrix == STMATMODE_INPLACE) {
alpha = st->sigma;
if( st->B ) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift( &alpha, st->A ); CHKERRQ(ierr); }
if( st->B ) {
ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,st->sigma); CHKERRQ(ierr);
}
st->setupcalled = 0;
}
PetscFunctionReturn(0);
113,7 → 114,6
PetscErrorCode STSetUp_Sinvert(ST st)
{
PetscErrorCode ierr;
PetscScalar alpha;
 
PetscFunctionBegin;
123,11 → 123,10
case STMATMODE_INPLACE:
st->mat = PETSC_NULL;
if (st->sigma != 0.0) {
alpha = -st->sigma;
if (st->B) {
ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr);
ierr = MatAXPY(st->A,-st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(&alpha,st->A);CHKERRQ(ierr);
ierr = MatShift(st->A,-st->sigma);CHKERRQ(ierr);
}
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
141,11 → 140,10
default:
ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
if (st->sigma != 0.0) {
alpha = -st->sigma;
if (st->B) {
ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr);
ierr = MatAXPY(st->mat,-st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr);
ierr = MatShift(st->mat,-st->sigma);CHKERRQ(ierr);
}
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
161,7 → 159,6
PetscErrorCode STSetShift_Sinvert(ST st,PetscScalar newshift)
{
PetscErrorCode ierr;
PetscScalar alpha;
 
PetscFunctionBegin;
 
172,15 → 169,19
case STMATMODE_INPLACE:
/* Undo previous operations */
if (st->sigma != 0.0) {
alpha = st->sigma;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,st->sigma);CHKERRQ(ierr);
}
}
/* Apply new shift */
if (newshift != 0.0) {
alpha = -newshift;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->A,-newshift,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->A,-newshift);CHKERRQ(ierr);
}
}
ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
break;
190,9 → 191,11
default:
ierr = MatCopy(st->A, st->mat,SUBSET_NONZERO_PATTERN); CHKERRQ(ierr);
if (newshift != 0.0) {
alpha = -newshift;
if (st->B) { ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr); }
else { ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr); }
if (st->B) {
ierr = MatAXPY(st->mat,-newshift,st->B,st->str);CHKERRQ(ierr);
} else {
ierr = MatShift(st->mat,-newshift);CHKERRQ(ierr);
}
}
/* In the following line, the SAME_NONZERO_PATTERN flag has been used to
* improve performance when solving a number of related eigenproblems */
/trunk/src/st/impls/shift/shift.c
21,7 → 21,7
ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
}
if (st->sigma != 0.0) {
ierr = VecAXPY(&st->sigma,x,y);CHKERRQ(ierr);
ierr = VecAXPY(y,st->sigma,x);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
43,7 → 43,7
ierr = MatMultTranspose(st->A,x,y);CHKERRQ(ierr);
}
if (st->sigma != 0.0) {
ierr = VecAXPY(&st->sigma,x,y);CHKERRQ(ierr);
ierr = VecAXPY(y,st->sigma,x);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
/trunk/src/eps/interface/orthog.c
43,7 → 43,6
{
PetscErrorCode ierr;
int k;
PetscScalar alpha;
PetscReal norm;
PetscTruth lindep;
61,8 → 60,7
ierr = SlepcVecSetRandom(V[k]);CHKERRQ(ierr);
ierr = STNorm(eps->OP,V[k],&norm);CHKERRQ(ierr);
}
alpha = 1.0/norm;
ierr = VecScale(&alpha,V[k]);CHKERRQ(ierr);
ierr = VecScale(V[k],1.0/norm);CHKERRQ(ierr);
if (R) R[k+ldr*k] = norm;
 
}
79,8 → 77,7
{
PetscErrorCode ierr;
int j;
PetscScalar shh[100],*lhh,
zero = 0.0,minus = -1.0;
PetscScalar shh[100],*lhh;
Vec w;
 
PetscFunctionBegin;
97,9 → 94,9
/* h = W^* v */
/* q = v - V h */
ierr = STMInnerProduct(eps->OP,n,v,W,H);CHKERRQ(ierr);
ierr = VecSet(&zero,w);CHKERRQ(ierr);
ierr = VecMAXPY(n,H,w,V);CHKERRQ(ierr);
ierr = VecAXPY(&minus,w,v);CHKERRQ(ierr);
ierr = VecSet(w,0.0);CHKERRQ(ierr);
ierr = VecMAXPY(w,n,H,V);CHKERRQ(ierr);
ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);
/* compute hnorm */
if (hnorm) {
130,9 → 127,9
for (j=0;j<n;j++) {
H[j] += lhh[j];
}
ierr = VecSet(&zero,w);CHKERRQ(ierr);
ierr = VecMAXPY(n,lhh,w,V);CHKERRQ(ierr);
ierr = VecAXPY(&minus,w,v);CHKERRQ(ierr);
ierr = VecSet(w,0.0);CHKERRQ(ierr);
ierr = VecMAXPY(w,n,lhh,V);CHKERRQ(ierr);
ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);
 
if (hnorm) *hnorm = *norm;
}
168,8 → 165,7
/* store coefficients if requested */
H[j] = alpha;
/* v <- v - alpha v_j */
alpha = -alpha;
ierr = VecAXPY(&alpha,V[j],v);CHKERRQ(ierr);
ierr = VecAXPY(v,-alpha,V[j]);CHKERRQ(ierr);
}
/* compute hnorm */
200,8 → 196,7
/* store coefficients if requested */
H[j] += alpha;
/* v <- v - alpha v_j */
alpha = -alpha;
ierr = VecAXPY(&alpha,V[j],v);CHKERRQ(ierr);
ierr = VecAXPY(v,-alpha,V[j]);CHKERRQ(ierr);
}
if (hnorm) *hnorm = *norm;
}
/trunk/src/eps/interface/solve.c
81,12 → 81,11
#ifndef PETSC_USE_COMPLEX
/* reorder conjugate eigenvalues (positive imaginary first) */
for (i=0; i<eps->nconv-1; i++) {
PetscScalar minus = -1.0;
if (eps->eigi[i] != 0) {
if (eps->eigi[i] < 0) {
eps->eigi[i] = -eps->eigi[i];
eps->eigi[i+1] = -eps->eigi[i+1];
ierr = VecScale(&minus, eps->V[i+1]); CHKERRQ(ierr);
ierr = VecScale(eps->V[i+1],-1.0); CHKERRQ(ierr);
}
i++;
}
504,10 → 503,6
{
PetscErrorCode ierr;
int k;
PetscScalar zero = 0.0;
#ifndef PETSC_USE_COMPLEX
PetscScalar minus = -1.0;
#endif
 
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
525,7 → 520,7
else k = eps->perm[i];
#ifdef PETSC_USE_COMPLEX
if (Vr) { ierr = VecCopy(eps->AV[k], Vr); CHKERRQ(ierr); }
if (Vi) { ierr = VecSet(&zero, Vi); CHKERRQ(ierr); }
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
#else
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
if (Vr) { ierr = VecCopy(eps->AV[k], Vr); CHKERRQ(ierr); }
534,11 → 529,11
if (Vr) { ierr = VecCopy(eps->AV[k-1], Vr); CHKERRQ(ierr); }
if (Vi) {
ierr = VecCopy(eps->AV[k], Vi); CHKERRQ(ierr);
ierr = VecScale(&minus, Vi); CHKERRQ(ierr);
ierr = VecScale(Vi,-1.0); CHKERRQ(ierr);
}
} else { /* real eigenvalue */
if (Vr) { ierr = VecCopy(eps->AV[k], Vr); CHKERRQ(ierr); }
if (Vi) { ierr = VecSet(&zero, Vi); CHKERRQ(ierr); }
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
}
#endif
579,10 → 574,6
{
PetscErrorCode ierr;
int k;
PetscScalar zero = 0.0;
#ifndef PETSC_USE_COMPLEX
PetscScalar minus = -1.0;
#endif
 
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
604,7 → 595,7
else k = eps->perm[i];
#ifdef PETSC_USE_COMPLEX
if (Wr) { ierr = VecCopy(eps->AW[k], Wr); CHKERRQ(ierr); }
if (Wi) { ierr = VecSet(&zero, Wi); CHKERRQ(ierr); }
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
#else
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
if (Wr) { ierr = VecCopy(eps->AW[k], Wr); CHKERRQ(ierr); }
613,11 → 604,11
if (Wr) { ierr = VecCopy(eps->AW[k-1], Wr); CHKERRQ(ierr); }
if (Wi) {
ierr = VecCopy(eps->AW[k], Wi); CHKERRQ(ierr);
ierr = VecScale(&minus, Wi); CHKERRQ(ierr);
ierr = VecScale(Wi,-1.0); CHKERRQ(ierr);
}
} else { /* real eigenvalue */
if (Wr) { ierr = VecCopy(eps->AW[k], Wr); CHKERRQ(ierr); }
if (Wi) { ierr = VecSet(&zero, Wi); CHKERRQ(ierr); }
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
}
#endif
736,7 → 727,7
PetscErrorCode ierr;
Vec u, v, w, xr, xi;
Mat A, B;
PetscScalar alpha, kr, ki;
PetscScalar kr, ki;
#ifndef PETSC_USE_COMPLEX
PetscReal ni, nr;
#endif
759,8 → 750,7
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
if (eps->isgeneralized) { ierr = MatMult( B, xr, w ); CHKERRQ(ierr); }
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B*x */
alpha = -kr;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A*x-k*B*x */
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*x-k*B*x */
}
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
#ifndef PETSC_USE_COMPLEX
768,18 → 758,14
ierr = MatMult( A, xr, u ); CHKERRQ(ierr); /* u=A*xr */
if (eps->isgeneralized) { ierr = MatMult( B, xr, v ); CHKERRQ(ierr); }
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B*xr */
alpha = -kr;
ierr = VecAXPY( &alpha, v, u ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr */
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr */
if (eps->isgeneralized) { ierr = MatMult( B, xi, w ); CHKERRQ(ierr); }
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B*xi */
alpha = ki;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr+ki*B*xi */
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A*xr-kr*B*xr+ki*B*xi */
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
ierr = MatMult( A, xi, u ); CHKERRQ(ierr); /* u=A*xi */
alpha = -kr;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi */
alpha = -ki;
ierr = VecAXPY( &alpha, v, u ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi-ki*B*xr */
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi */
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A*xi-kr*B*xi-ki*B*xr */
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
*norm = SlepcAbsEigenvalue( nr, ni );
}
824,7 → 810,7
PetscErrorCode ierr;
Vec u, v, w, xr, xi;
Mat A, B;
PetscScalar alpha, kr, ki;
PetscScalar kr, ki;
#ifndef PETSC_USE_COMPLEX
PetscReal ni, nr;
#endif
848,8 → 834,7
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, w ); CHKERRQ(ierr); }
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B'*x */
alpha = -kr;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A'*x-k*B'*x */
ierr = VecAXPY( u, -kr, w); CHKERRQ(ierr); /* u=A'*x-k*B'*x */
}
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
#ifndef PETSC_USE_COMPLEX
857,18 → 842,14
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*xr */
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, v ); CHKERRQ(ierr); }
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B'*xr */
alpha = -kr;
ierr = VecAXPY( &alpha, v, u ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xi, w ); CHKERRQ(ierr); }
else { ierr = VecCopy( xi, w ); CHKERRQ(ierr); } /* w=B'*xi */
alpha = ki;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr+ki*B'*xi */
ierr = VecAXPY( u, ki, w ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr+ki*B'*xi */
ierr = VecNorm( u, NORM_2, &nr ); CHKERRQ(ierr);
ierr = MatMultTranspose( A, xi, u ); CHKERRQ(ierr); /* u=A'*xi */
alpha = -kr;
ierr = VecAXPY( &alpha, w, u ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi */
alpha = -ki;
ierr = VecAXPY( &alpha, v, u ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi-ki*B'*xr */
ierr = VecAXPY( u, -kr, w ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi */
ierr = VecAXPY( u, -ki, v ); CHKERRQ(ierr); /* u=A'*xi-kr*B'*xi-ki*B'*xr */
ierr = VecNorm( u, NORM_2, &ni ); CHKERRQ(ierr);
*norm = SlepcAbsEigenvalue( nr, ni );
}
911,7 → 892,6
PetscReal norm, er;
#ifndef PETSC_USE_COMPLEX
Vec u;
PetscScalar alpha;
PetscReal ei;
#endif
927,7 → 907,7
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
#endif
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
ierr = VecScale(&kr, xr); CHKERRQ(ierr);
ierr = VecScale(xr, kr); CHKERRQ(ierr);
}
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
*error = norm / er;
935,10 → 915,9
} else {
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
ierr = VecCopy(xi, u); CHKERRQ(ierr);
alpha = -ki;
ierr = VecAXPBY(&kr, &alpha, xr, u); CHKERRQ(ierr);
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
ierr = VecAXPBY(&kr, &ki, xr, xi); CHKERRQ(ierr);
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
ierr = VecDestroy(u); CHKERRQ(ierr);
*error = norm / SlepcAbsEigenvalue(er, ei);
980,7 → 959,6
PetscReal norm, er;
#ifndef PETSC_USE_COMPLEX
Vec u;
PetscScalar alpha;
PetscReal ei;
#endif
997,7 → 975,7
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
#endif
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
ierr = VecScale(&kr, xr); CHKERRQ(ierr);
ierr = VecScale(xr, kr); CHKERRQ(ierr);
}
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
*error = norm / er;
1005,10 → 983,9
} else {
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
ierr = VecCopy(xi, u); CHKERRQ(ierr);
alpha = -ki;
ierr = VecAXPBY(&kr, &alpha, xr, u); CHKERRQ(ierr);
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
ierr = VecAXPBY(&kr, &ki, xr, xi); CHKERRQ(ierr);
ierr = VecAXPBY(xi, kr, ki, xr); CHKERRQ(ierr);
ierr = VecNorm(xi, NORM_2, &ei); CHKERRQ(ierr);
ierr = VecDestroy(u); CHKERRQ(ierr);
*error = norm / SlepcAbsEigenvalue(er, ei);
1044,13 → 1021,12
{
PetscErrorCode ierr;
int i;
PetscScalar zero = 0.0;
PetscFunctionBegin;
ierr = PetscLogEventBegin(EPS_ReverseProjection,eps,0,0,0);CHKERRQ(ierr);
for (i=k;i<m;i++) {
ierr = VecSet(&zero,work[i]);CHKERRQ(ierr);
ierr = VecMAXPY(m,S+m*i,work[i],V);CHKERRQ(ierr);
ierr = VecSet(work[i],0.0);CHKERRQ(ierr);
ierr = VecMAXPY(work[i],m,S+m*i,V);CHKERRQ(ierr);
}
for (i=k;i<m;i++) {
ierr = VecCopy(work[i],V[i]);CHKERRQ(ierr);
1187,7 → 1163,6
{
PetscErrorCode ierr;
PetscTruth breakdown;
PetscScalar t;
PetscReal norm;
Vec w;
1213,8 → 1188,7
if (i==0) { SETERRQ(1,"Initial vector is zero or belongs to the deflation space"); }
else { SETERRQ(1,"Unable to generate more start vectors"); }
}
t = 1 / norm;
ierr = VecScale(&t,vec);CHKERRQ(ierr);
ierr = VecScale(vec,1/norm);CHKERRQ(ierr);
 
if (i!=0) {
ierr = VecDestroy(w);CHKERRQ(ierr);
1255,7 → 1229,6
{
PetscErrorCode ierr;
PetscTruth breakdown;
PetscScalar t;
PetscReal norm;
Vec w;
1281,8 → 1254,7
if (i==0) { SETERRQ(1,"Left initial vector is zero"); }
else { SETERRQ(1,"Unable to generate more left start vectors"); }
}
t = 1 / norm;
ierr = VecScale(&t,vec);CHKERRQ(ierr);
ierr = VecScale(vec,1/norm);CHKERRQ(ierr);
 
if (i!=0) {
ierr = VecDestroy(w);CHKERRQ(ierr);
/trunk/src/eps/impls/subspace/subspace.c
163,16 → 163,15
{
PetscErrorCode ierr;
int i;
PetscScalar zero = 0.0,minus = -1.0;
#if defined(PETSC_USE_COMPLEX)
PetscScalar t;
#endif
 
PetscFunctionBegin;
for (i=l;i<m;i++) {
ierr = VecSet(&zero,eps->work[0]);CHKERRQ(ierr);
ierr = VecMAXPY(m,T+ldt*i,eps->work[0],V);CHKERRQ(ierr);
ierr = VecWAXPY(&minus,eps->work[0],AV[i],eps->work[1]);CHKERRQ(ierr);
ierr = VecSet(eps->work[0],0.0);CHKERRQ(ierr);
ierr = VecMAXPY(eps->work[0],m,T+ldt*i,V);CHKERRQ(ierr);
ierr = VecWAXPY(eps->work[1],-1.0,eps->work[0],AV[i]);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
ierr = VecDot(eps->work[1],eps->work[1],rsd+i);CHKERRQ(ierr);
#else
205,7 → 204,7
PetscErrorCode ierr;
int i,j,ilo,lwork,info,ngrp,nogrp,*itrsd,*itrsdold,
nxtsrr,idsrr,*iwork,idort,nxtort,ncv = eps->ncv;
PetscScalar *T=eps->T,*U,*tau,*work,t;
PetscScalar *T=eps->T,*U,*tau,*work;
PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,*rsdold,norm,tcond;
PetscTruth breakdown;
/* Parameters */
347,8 → 346,7
for (i=eps->nconv;i<ncv;i++) {
ierr = VecCopy(eps->AV[i],eps->V[i]);CHKERRQ(ierr);
ierr = VecNorm(eps->V[i],NORM_INFINITY,&norm);CHKERRQ(ierr);
t = 1 / norm;
ierr = VecScale(&t,eps->V[i]);CHKERRQ(ierr);
ierr = VecScale(eps->V[i],1/norm);CHKERRQ(ierr);
}
eps->its++;
360,8 → 358,7
ierr = SlepcVecSetRandom(eps->V[i]);CHKERRQ(ierr);
ierr = EPSOrthogonalize(eps,i+eps->nds,eps->DSV,eps->V[i],PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
}
t = 1 / norm;
ierr = VecScale(&t,eps->V[i]);CHKERRQ(ierr);
ierr = VecScale(eps->V[i],1/norm);CHKERRQ(ierr);
}
nxtort = PetscMin(eps->its+idort,nxtsrr);
} while (eps->its<nxtsrr);
/trunk/src/eps/impls/arnoldi/arnoldi.c
76,7 → 76,6
PetscErrorCode ierr;
int j;
PetscReal norm;
PetscScalar t;
PetscTruth breakdown;
 
PetscFunctionBegin;
89,11 → 88,9
if (breakdown) {
PetscLogInfo((eps,"Breakdown in Arnoldi method (norm=%g)\n",norm));
ierr = EPSGetStartVector(eps,j,V[j+1]);CHKERRQ(ierr);
} else {
ierr = VecScale(V[j+1],1/norm);CHKERRQ(ierr);
}
else {
t = 1 / norm;
ierr = VecScale(&t,V[j+1]);CHKERRQ(ierr);
}
}
ierr = STApply(eps->OP,V[m-1],f);CHKERRQ(ierr);
eps->its++;
/trunk/src/eps/impls/lanczos/lanczos.c
89,7 → 89,6
PetscErrorCode ierr;
int j,m = *M;
PetscReal norm;
PetscScalar t;
PetscTruth breakdown;
 
PetscFunctionBegin;
105,8 → 104,7
ierr = EPSGetStartVector(eps,j+1,V[j+1]);CHKERRQ(ierr);
} else {
T[m*j+j+1] = norm;
t = 1 / norm;
ierr = VecScale(&t,V[j+1]);CHKERRQ(ierr);
ierr = VecScale(V[j+1],1/norm);CHKERRQ(ierr);
}
}
ierr = STApply(eps->OP,V[m-1],f);CHKERRQ(ierr);
131,7 → 129,6
PetscErrorCode ierr;
int j,m = *M;
PetscReal norm;
PetscScalar t;
PetscTruth breakdown = PETSC_FALSE;
 
PetscFunctionBegin;
150,8 → 147,7
PetscLogInfo((eps,"Breakdown in Lanczos method (norm=%g)\n",norm));
ierr = EPSGetStartVector(eps,j+1,V[j+1]);CHKERRQ(ierr);
} else {
t = 1 / norm;
ierr = VecScale(&t,V[j+1]);CHKERRQ(ierr);
ierr = VecScale(V[j+1],1/norm);CHKERRQ(ierr);
}
}
ierr = STApply(eps->OP,V[m-1],f);CHKERRQ(ierr);
196,7 → 192,6
PetscErrorCode ierr;
int i,j,n,m = *M,info;
PetscReal *omega,*oldomega,*a,*b,delta,eps1,*Y,*work;
PetscScalar t;
PetscTruth breakdown = PETSC_FALSE;
 
PetscFunctionBegin;
259,9 → 254,8
}
}
}
t = 1 / *beta;
ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
ierr = VecScale(&t,V[j+1]);CHKERRQ(ierr);
ierr = VecScale(V[j+1],1 / *beta);CHKERRQ(ierr);
}
ierr = STApply(eps->OP,V[m-1],f);CHKERRQ(ierr);
eps->its++;
280,7 → 274,7
PetscErrorCode ierr;
int i,j,l,n,m = *M;
PetscReal norm,*omega,*oldomega,*a,*b,delta,eps1;
PetscScalar *H,t;
PetscScalar *H;
PetscTruth breakdown,reorthog;
 
PetscFunctionBegin;
326,8 → 320,7
ierr = STNorm(eps->OP,V[j+1],&norm);CHKERRQ(ierr);
} else {
T[m*j+j+1] = norm;
t = 1 / norm;
ierr = VecScale(&t,V[j+1]);CHKERRQ(ierr);
ierr = VecScale(V[j+1],1/norm);CHKERRQ(ierr);
}
}
442,7 → 435,7
int nconv,i,j,n,m,info,N,*perm,
ncv=eps->ncv;
Vec f=eps->work[ncv];
PetscScalar *T=eps->T,*Y,*W,ts;
PetscScalar *T=eps->T,*Y,*W;
PetscReal *ritz,*bnd,*E,*work,norm,anorm,beta;
 
PetscFunctionBegin;
564,8 → 557,7
ierr = STNorm(eps->OP,eps->V[nconv],&norm);CHKERRQ(ierr);
eps->errest[nconv] = eps->errest[nconv] / norm;
if (eps->errest[nconv]<eps->tol) {
ts = 1 / norm;
ierr = VecScale(&ts,eps->V[nconv]);CHKERRQ(ierr);
ierr = VecScale(eps->V[nconv],1/norm);CHKERRQ(ierr);
nconv++;
}
}
/trunk/src/eps/impls/lapack/lapack.c
40,7 → 40,7
} else la->B = PETSC_NULL;
ierr = STGetShift(eps->OP,&shift);CHKERRQ(ierr);
if (shift != 0.0) {
ierr = MatShift(&shift,la->A);CHKERRQ(ierr);
ierr = MatShift(la->A,shift);CHKERRQ(ierr);
}
} else {
PetscLogInfo((eps,"EPSSetup_LAPACK: Using slow explicit operator\n"));
/trunk/src/eps/impls/power/power.c
86,7 → 86,7
Vec e, z;
Mat A;
PetscReal norm, rt1, rt2, cs1;
PetscScalar alpha, alpha1, alpha2, beta1, sn1;
PetscScalar alpha1, alpha2, beta1, sn1;
 
PetscFunctionBegin;
e = eps->work[0];
103,8 → 103,7
SETERRQ(PETSC_ERR_SUP,"LAEV2 - Lapack routine is unavailable.");
#else
/* beta1 is the norm of the residual associated to R(v) */
alpha = -alpha1;
ierr = VecAXPY(&alpha,v,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-alpha1,v);CHKERRQ(ierr);
ierr = STNorm(eps->OP,e,&norm);CHKERRQ(ierr);
beta1 = norm;
137,7 → 136,7
Vec v, y, e, *SV;
Mat A;
PetscReal relerr, norm, rt1, rt2, cs1, anorm;
PetscScalar theta, alpha, rho, delta, sigma, alpha2, beta1, sn1;
PetscScalar theta, rho, delta, sigma, alpha2, beta1, sn1;
 
PetscFunctionBegin;
v = eps->V[0];
177,8 → 176,7
 
/* compute relative error as ||y-theta v||_2/|theta| */
ierr = VecCopy(y,e);CHKERRQ(ierr);
alpha = -theta;
ierr = VecAXPY(&alpha,v,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-theta,v);CHKERRQ(ierr);
ierr = VecNorm(e,NORM_2,&norm);CHKERRQ(ierr);
relerr = norm / PetscAbsScalar(theta);
 
203,10 → 201,8
SETERRQ(PETSC_ERR_SUP,"LAEV2 - Lapack routine is unavailable.");
#else
/* beta1 is the norm of the residual associated to R(v) */
alpha = -theta/(delta*delta);
ierr = VecAXPY(&alpha,y,v);CHKERRQ(ierr);
alpha = 1.0/delta;
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
ierr = VecAXPY(v,-theta/(delta*delta),y);CHKERRQ(ierr);
ierr = VecScale(v,1.0/delta);CHKERRQ(ierr);
ierr = STNorm(eps->OP,v,&norm);CHKERRQ(ierr);
beta1 = norm;
245,8 → 241,7
 
/* v = y/||y||_B */
ierr = VecCopy(y,v);CHKERRQ(ierr);
alpha = 1.0/norm;
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
ierr = VecScale(v,1.0/norm);CHKERRQ(ierr);
 
/* if relerr<tol, accept eigenpair */
if (relerr<eps->tol) {
316,14 → 311,12
 
/* compute relative errors (right and left) */
ierr = VecCopy(y,e);CHKERRQ(ierr);
alpha = -theta;
ierr = VecAXPY(&alpha,v,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-theta,v);CHKERRQ(ierr);
ierr = VecNorm(e,NORM_2,&norm);CHKERRQ(ierr);
relerr = norm / PetscAbsScalar(theta);
eps->errest[eps->nconv] = relerr;
ierr = VecCopy(z,e);CHKERRQ(ierr);
alpha = -theta;
ierr = VecAXPY(&alpha,w,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-theta,w);CHKERRQ(ierr);
ierr = VecNorm(e,NORM_2,&norm);CHKERRQ(ierr);
relerr = norm / PetscAbsScalar(theta);
eps->errest_left[eps->nconv] = relerr;
353,10 → 346,8
SETERRQ(PETSC_ERR_SUP,"LAEV2 - Lapack routine is unavailable.");
#else
/* beta1 is the norm of the residual associated to R(v,w) */
alpha = -theta/(delta*delta);
ierr = VecAXPY(&alpha,y,v);CHKERRQ(ierr);
alpha = 1.0/delta;
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
ierr = VecAXPY(v,-theta/(delta*delta),y);CHKERRQ(ierr);
ierr = VecScale(v,1.0/delta);CHKERRQ(ierr);
ierr = STNorm(eps->OP,v,&norm);CHKERRQ(ierr);
beta1 = norm;
392,13 → 383,12
if (alpha==0.0) SETERRQ(1,"Breakdown in two-sided Power/RQI");
if (alpha>0.0) {
alpha = 1.0/PetscSqrtScalar(alpha);
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
ierr = VecScale(&alpha,w);CHKERRQ(ierr);
ierr = VecScale(v,alpha);CHKERRQ(ierr);
ierr = VecScale(w,alpha);CHKERRQ(ierr);
} else {
alpha = 1.0/PetscSqrtScalar(-alpha);
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
alpha = -alpha;
ierr = VecScale(&alpha,w);CHKERRQ(ierr);
ierr = VecScale(v,alpha);CHKERRQ(ierr);
ierr = VecScale(w,-alpha);CHKERRQ(ierr);
}
 
/* if relerr<tol (both right and left), accept eigenpair */
459,13 → 449,12
if (alpha==0.0) SETERRQ(1,"Breakdown in two-sided Power/RQI");
if (alpha>0.0) {
alpha = 1.0/PetscSqrtScalar(alpha);
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
ierr = VecScale(&alpha,w);CHKERRQ(ierr);
ierr = VecScale(v,alpha);CHKERRQ(ierr);
ierr = VecScale(w,alpha);CHKERRQ(ierr);
} else {
alpha = 1.0/PetscSqrtScalar(-alpha);
ierr = VecScale(&alpha,v);CHKERRQ(ierr);
alpha = -alpha;
ierr = VecScale(&alpha,w);CHKERRQ(ierr);
ierr = VecScale(v,alpha);CHKERRQ(ierr);
ierr = VecScale(w,-alpha);CHKERRQ(ierr);
}
 
/* y = OP v */
477,14 → 466,12
 
/* compute residual norm */
ierr = VecCopy(y,e);CHKERRQ(ierr);
alpha = -theta;
ierr = VecAXPY(&alpha,v,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-theta,v);CHKERRQ(ierr);
ierr = VecNorm(e,NORM_2,&relerr);CHKERRQ(ierr);
relerr = relerr / PetscAbsScalar(theta);
eps->errest[eps->nconv] = relerr;
ierr = VecCopy(z,e);CHKERRQ(ierr);
alpha = -theta;
ierr = VecAXPY(&alpha,w,e);CHKERRQ(ierr);
ierr = VecAXPY(e,-theta,w);CHKERRQ(ierr);
ierr = VecNorm(e,NORM_2,&relerr);CHKERRQ(ierr);
relerr = relerr / PetscAbsScalar(theta);
eps->errest_left[eps->nconv] = relerr;
/trunk/src/examples/ex2.c
31,7 → 31,8
Compute the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/trunk/src/examples/ex5.c
23,7 → 23,6
PetscReal error, tol, re, im;
PetscScalar kr, ki;
int N, m=15, nev, ierr, maxit, i, its, nconv;
PetscScalar alpha;
 
SlepcInitialize(&argc,&argv,(char*)0,help);
 
35,7 → 34,8
Compute the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatMarkovModel( m, A );CHKERRQ(ierr);
 
64,8 → 64,7
vector is set to random values
*/
ierr = MatGetVecs(A,&v0,PETSC_NULL);CHKERRQ(ierr);
alpha = 1.0;
ierr = VecSet(&alpha,v0);CHKERRQ(ierr);
ierr = VecSet(v0,1.0);CHKERRQ(ierr);
ierr = EPSSetInitialVector(eps,v0);CHKERRQ(ierr);
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
/trunk/src/examples/ex8.c
73,7 → 73,8
Generate the matrix
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
 
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/trunk/src/examples/ex9.c
81,7 → 81,8
/*
Create matrix T
*/
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&ctx->T);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&ctx->T);CHKERRQ(ierr);
ierr = MatSetSizes(ctx->T,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(ctx->T);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(ctx->T,&Istart,&Iend);CHKERRQ(ierr);
240,7 → 241,7
int MatBrussel_Mult(Mat A,Vec x,Vec y)
{
int n, ierr;
PetscScalar alpha, *px, *py;
PetscScalar *px, *py;
MPI_Comm comm;
CTX_BRUSSEL *ctx;
 
255,18 → 256,14
ierr = VecPlaceArray(ctx->y2,py+n);CHKERRQ(ierr);
 
ierr = MatMult(ctx->T,ctx->x1,ctx->y1);CHKERRQ(ierr);
ierr = VecScale(&ctx->tau1,ctx->y1);CHKERRQ(ierr);
alpha = ctx->beta - 1.0 + ctx->sigma;
ierr = VecAXPY(&alpha,ctx->x1,ctx->y1);CHKERRQ(ierr);
alpha = ctx->alpha * ctx->alpha;
ierr = VecAXPY(&alpha,ctx->x2,ctx->y1);CHKERRQ(ierr);
ierr = VecScale(ctx->y1,ctx->tau1);CHKERRQ(ierr);
ierr = VecAXPY(ctx->y1,ctx->beta - 1.0 + ctx->sigma,ctx->x1);CHKERRQ(ierr);
ierr = VecAXPY(ctx->y1,ctx->alpha * ctx->alpha,ctx->x2);CHKERRQ(ierr);
 
ierr = MatMult(ctx->T,ctx->x2,ctx->y2);CHKERRQ(ierr);
ierr = VecScale(&ctx->tau2,ctx->y2);CHKERRQ(ierr);
alpha = -ctx->beta;
ierr = VecAXPY(&alpha,ctx->x1,ctx->y2);CHKERRQ(ierr);
alpha = -ctx->alpha * ctx->alpha + ctx->sigma;
ierr = VecAXPY(&alpha,ctx->x2,ctx->y2);CHKERRQ(ierr);
ierr = VecScale(ctx->y2,ctx->tau2);CHKERRQ(ierr);
ierr = VecAXPY(ctx->y2,-ctx->beta,ctx->x1);CHKERRQ(ierr);
ierr = VecAXPY(ctx->y2,-ctx->alpha * ctx->alpha + ctx->sigma,ctx->x2);CHKERRQ(ierr);
 
ierr = VecRestoreArray(x,&px);CHKERRQ(ierr);
ierr = VecRestoreArray(y,&py);CHKERRQ(ierr);
292,7 → 289,7
{
Vec d1, d2;
int n, ierr;
PetscScalar alpha, *pd;
PetscScalar *pd;
MPI_Comm comm;
CTX_BRUSSEL *ctx;
 
303,10 → 300,8
ierr = VecCreateMPIWithArray(comm,n,PETSC_DECIDE,pd,&d1);CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(comm,n,PETSC_DECIDE,pd+n,&d2);CHKERRQ(ierr);
 
alpha = -2.0*ctx->tau1 + ctx->beta - 1.0 + ctx->sigma;
ierr = VecSet(&alpha,d1);CHKERRQ(ierr);
alpha = -2.0*ctx->tau2 - ctx->alpha*ctx->alpha + ctx->sigma;
ierr = VecSet(&alpha,d2);CHKERRQ(ierr);
ierr = VecSet(d1,-2.0*ctx->tau1 + ctx->beta - 1.0 + ctx->sigma);CHKERRQ(ierr);
ierr = VecSet(d2,-2.0*ctx->tau2 - ctx->alpha*ctx->alpha + ctx->sigma);CHKERRQ(ierr);
 
ierr = VecDestroy(d1);CHKERRQ(ierr);
ierr = VecDestroy(d2);CHKERRQ(ierr);
/trunk/src/examples/ex1f.F
62,8 → 62,8
! Compute the operator matrix that defines the eigensystem, Ax=kx
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,A,
& ierr)
call MatCreate(PETSC_COMM_WORLD,A,ierr)
call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
call MatSetFromOptions(A,ierr)
 
call MatGetOwnershipRange(A,Istart,Iend,ierr)
/trunk/src/examples/ex10.c
44,7 → 44,8
Compute the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/trunk/src/examples/ex11.c
19,7 → 19,7
PetscReal error, tol, re, im;
PetscScalar kr, ki;
int N, n=10, m, nev, ierr, maxit, i, j, I, J, its, nconv, Istart, Iend;
PetscScalar v, w, one = 1.0;
PetscScalar v, w;
PetscTruth flag;
 
SlepcInitialize(&argc,&argv,(char*)0,help);
36,7 → 36,8
Lii = degree of node i, Lij = -1 if edge (i,j) exists in G
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
86,7 → 87,7
nullspace, [1 1 ... 1]^T is the eigenvector of the zero eigenvalue
*/
ierr = MatGetVecs(A,&x,PETSC_NULL);CHKERRQ(ierr);
ierr = VecSet(&one,x);CHKERRQ(ierr);
ierr = VecSet(x,1.0);CHKERRQ(ierr);
ierr = EPSAttachDeflationSpace(eps,1,&x,PETSC_FALSE);CHKERRQ(ierr);
ierr = VecDestroy(x);
 
/trunk/src/examples/ex12.c
33,7 → 33,8
Compute the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatMarkovModel( m, A );CHKERRQ(ierr);
 
/trunk/src/examples/ex1.c
28,7 → 28,8
Compute the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/trunk/src/sys/slepcutil.c
75,7 → 75,6
PetscErrorCode ierr;
int M,N,m,n;
Vec x,w1,w2;
PetscScalar alpha;
MPI_Comm comm;
PetscReal norm;
PetscTruth has;
108,8 → 107,7
ierr = MatMult(A,x,w1);CHKERRQ(ierr);
ierr = MatMultTranspose(A,x,w2);CHKERRQ(ierr);
ierr = VecConjugate(w2);CHKERRQ(ierr);
alpha = -1.0;
ierr = VecAXPY(&alpha,w1,w2);CHKERRQ(ierr);
ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr);
ierr = VecNorm(w2,NORM_2,&norm);CHKERRQ(ierr);
if (norm<1.0e-6) *is = PETSC_TRUE;
ierr = VecDestroy(x);CHKERRQ(ierr);