| Line 176... |
Line 176... |
}
|
}
|
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
|
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
|
ierr = VecDestroy(&w);CHKERRQ(ierr);
|
ierr = VecDestroy(&w);CHKERRQ(ierr);
|
}
|
}
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
/* reorder conjugate eigenvalues (positive imaginary first) */
|
/* reorder conjugate eigenvalues (positive imaginary first) */
|
for (i=0; i<eps->nconv-1; i++) {
|
for (i=0; i<eps->nconv-1; i++) {
|
if (eps->eigi[i] != 0) {
|
if (eps->eigi[i] != 0) {
|
if (eps->eigi[i] < 0) {
|
if (eps->eigi[i] < 0) {
|
eps->eigi[i] = -eps->eigi[i];
|
eps->eigi[i] = -eps->eigi[i];
|
| Line 609... |
Line 609... |
if (i<0 || i>=eps->nconv) {
|
if (i<0 || i>=eps->nconv) {
|
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#if defined(PETSC_USE_COMPLEX)
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigi) *eigi = 0;
|
if (eigi) *eigi = 0;
|
#else
|
#else
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigr) *eigr = eps->eigr[k];
|
if (eigi) *eigi = eps->eigi[k];
|
if (eigi) *eigi = eps->eigi[k];
|
| Line 670... |
Line 670... |
if (!eps->evecsavailable && (Vr || Vi) ) {
|
if (!eps->evecsavailable && (Vr || Vi) ) {
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
}
|
}
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#if defined(PETSC_USE_COMPLEX)
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
if (Vi) { ierr = VecSet(Vi,0.0); CHKERRQ(ierr); }
|
#else
|
#else
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
if (Vr) { ierr = VecCopy(eps->V[k], Vr); CHKERRQ(ierr); }
|
| Line 742... |
Line 742... |
if (!eps->evecsavailable && (Wr || Wi) ) {
|
if (!eps->evecsavailable && (Wr || Wi) ) {
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
ierr = (*eps->ops->computevectors)(eps);CHKERRQ(ierr);
|
}
|
}
|
if (!eps->perm) k = i;
|
if (!eps->perm) k = i;
|
else k = eps->perm[i];
|
else k = eps->perm[i];
|
#ifdef PETSC_USE_COMPLEX
|
#if defined(PETSC_USE_COMPLEX)
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
if (Wi) { ierr = VecSet(Wi,0.0); CHKERRQ(ierr); }
|
#else
|
#else
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (eps->eigi[k] > 0) { /* first value of conjugate pair */
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
if (Wr) { ierr = VecCopy(eps->W[k], Wr); CHKERRQ(ierr); }
|
| Line 855... |
Line 855... |
PetscErrorCode EPSComputeResidualNorm_Private(EPS eps, PetscScalar kr, PetscScalar ki, Vec xr, Vec xi, PetscReal *norm)
|
PetscErrorCode EPSComputeResidualNorm_Private(EPS eps, PetscScalar kr, PetscScalar ki, Vec xr, Vec xi, PetscReal *norm)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec u, w;
|
Vec u, w;
|
Mat A, B;
|
Mat A, B;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
Vec v;
|
Vec v;
|
PetscReal ni, nr;
|
PetscReal ni, nr;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->V[0],&u);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->V[0],&u);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = MatMult(A,xr,u);CHKERRQ(ierr); /* u=A*x */
|
ierr = MatMult(A,xr,u);CHKERRQ(ierr); /* u=A*x */
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (eps->isgeneralized) { ierr = MatMult(B,xr,w);CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMult(B,xr,w);CHKERRQ(ierr); }
|
else { ierr = VecCopy(xr,w);CHKERRQ(ierr); } /* w=B*x */
|
else { ierr = VecCopy(xr,w);CHKERRQ(ierr); } /* w=B*x */
|
ierr = VecAXPY(u,-kr,w);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);
|
ierr = VecNorm(u,NORM_2,norm);CHKERRQ(ierr);
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
} else {
|
} else {
|
ierr = VecDuplicate(eps->V[0],&v); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->V[0],&v); CHKERRQ(ierr);
|
ierr = MatMult(A,xr,u);CHKERRQ(ierr); /* u=A*xr */
|
ierr = MatMult(A,xr,u);CHKERRQ(ierr); /* u=A*xr */
|
if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
|
if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
|
if (eps->isgeneralized) { ierr = MatMult(B,xr,v);CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMult(B,xr,v);CHKERRQ(ierr); }
|
| Line 981... |
Line 981... |
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec u, v, w, xr, xi;
|
Vec u, v, w, xr, xi;
|
Mat A, B;
|
Mat A, B;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
PetscReal ni, nr;
|
PetscReal ni, nr;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
|
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
|
| Line 999... |
Line 999... |
ierr = VecDuplicate(eps->W[0],&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenvalue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetEigenvalue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetEigenvectorLeft(eps,i,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenvectorLeft(eps,i,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*x */
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*x */
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, w ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, w ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B'*x */
|
else { ierr = VecCopy( xr, w ); CHKERRQ(ierr); } /* w=B'*x */
|
ierr = VecAXPY( u, -kr, w); 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);
|
ierr = VecNorm( u, NORM_2, norm); CHKERRQ(ierr);
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
} else {
|
} else {
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*xr */
|
ierr = MatMultTranspose( A, xr, u ); CHKERRQ(ierr); /* u=A'*xr */
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, v ); CHKERRQ(ierr); }
|
if (eps->isgeneralized) { ierr = MatMultTranspose( B, xr, v ); CHKERRQ(ierr); }
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B'*xr */
|
else { ierr = VecCopy( xr, v ); CHKERRQ(ierr); } /* v=B'*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
|
ierr = VecAXPY( u, -kr, v ); CHKERRQ(ierr); /* u=A'*xr-kr*B'*xr */
|
| Line 1046... |
Line 1046... |
*/
|
*/
|
PetscErrorCode EPSComputeRelativeError_Private(EPS eps, PetscScalar kr, PetscScalar ki, Vec xr, Vec xi, PetscReal *error)
|
PetscErrorCode EPSComputeRelativeError_Private(EPS eps, PetscScalar kr, PetscScalar ki, Vec xr, Vec xi, PetscReal *error)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal norm, er;
|
PetscReal norm, er;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
PetscReal ei;
|
PetscReal ei;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = EPSComputeResidualNorm_Private(eps,kr,ki,xr,xi,&norm);CHKERRQ(ierr);
|
ierr = EPSComputeResidualNorm_Private(eps,kr,ki,xr,xi,&norm);CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
if (ki == 0) {
|
if (ki == 0) {
|
#endif
|
#endif
|
ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
|
ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
} else {
|
} else {
|
ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
|
ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
|
ierr = VecNorm(xi,NORM_2,&ei);CHKERRQ(ierr);
|
ierr = VecNorm(xi,NORM_2,&ei);CHKERRQ(ierr);
|
er = SlepcAbsEigenvalue(er,ei);
|
er = SlepcAbsEigenvalue(er,ei);
|
}
|
}
|
| Line 1135... |
Line 1135... |
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec xr, xi;
|
Vec xr, xi;
|
PetscScalar kr, ki;
|
PetscScalar kr, ki;
|
PetscReal norm, er;
|
PetscReal norm, er;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
Vec u;
|
Vec u;
|
PetscReal ei;
|
PetscReal ei;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
| Line 1148... |
Line 1148... |
ierr = VecDuplicate(eps->W[0],&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xr); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xi); CHKERRQ(ierr);
|
ierr = VecDuplicate(eps->W[0],&xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenvalue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetEigenvalue(eps,i,&kr,&ki); CHKERRQ(ierr);
|
ierr = EPSGetEigenvectorLeft(eps,i,xr,xi); CHKERRQ(ierr);
|
ierr = EPSGetEigenvectorLeft(eps,i,xr,xi); CHKERRQ(ierr);
|
|
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
if (ki == 0 ||
|
if (ki == 0 ||
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
|
#endif
|
#endif
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(xr, NORM_2, &er); CHKERRQ(ierr);
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
|
*error = norm / (PetscAbsScalar(kr) * er);
|
*error = norm / (PetscAbsScalar(kr) * er);
|
} else {
|
} else {
|
*error = norm / er;
|
*error = norm / er;
|
}
|
}
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
} else {
|
} else {
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecDuplicate(xi, &u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecCopy(xi, u); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecAXPBY(u, kr, -ki, xr); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
ierr = VecNorm(u, NORM_2, &er); CHKERRQ(ierr);
|
| Line 1217... |
Line 1217... |
/* insertion sort */
|
/* insertion sort */
|
for (i=n-1; i>=0; i--) {
|
for (i=n-1; i>=0; i--) {
|
re = eigr[perm[i]];
|
re = eigr[perm[i]];
|
im = eigi[perm[i]];
|
im = eigi[perm[i]];
|
j = i + 1;
|
j = i + 1;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
if (im != 0) {
|
if (im != 0) {
|
/* complex eigenvalue */
|
/* complex eigenvalue */
|
i--;
|
i--;
|
im = eigi[perm[i]];
|
im = eigi[perm[i]];
|
}
|
}
|
#endif
|
#endif
|
while (j<n) {
|
while (j<n) {
|
ierr = EPSCompareEigenvalues(eps,re,im,eigr[perm[j]],eigi[perm[j]],&result);CHKERRQ(ierr);
|
ierr = EPSCompareEigenvalues(eps,re,im,eigr[perm[j]],eigi[perm[j]],&result);CHKERRQ(ierr);
|
if (result < 0) break;
|
if (result < 0) break;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
/* keep together every complex conjugated eigenpair */
|
/* keep together every complex conjugated eigenpair */
|
if (im == 0) {
|
if (im == 0) {
|
if (eigi[perm[j]] == 0) {
|
if (eigi[perm[j]] == 0) {
|
#endif
|
#endif
|
tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
|
tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = tmp;
|
j++;
|
j++;
|
#ifndef PETSC_USE_COMPLEX
|
#if !defined(PETSC_USE_COMPLEX)
|
} else {
|
} else {
|
tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
|
tmp = perm[j-1]; perm[j-1] = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp;
|
j+=2;
|
j+=2;
|
}
|
}
|
} else {
|
} else {
|