| Line 1201... |
Line 1201... |
/* Else orthogonalize first against cX and then against V */
|
/* Else orthogonalize first against cX and then against V */
|
ierr = IPOrthogonalize(ip, size_cX, cX, i, PETSC_NULL, V,
|
ierr = IPOrthogonalize(ip, size_cX, cX, i, PETSC_NULL, V,
|
V[i], auxS0, &norm, &lindep); CHKERRQ(ierr);
|
V[i], auxS0, &norm, &lindep); CHKERRQ(ierr);
|
}
|
}
|
if(!lindep && (norm > PETSC_MACHINE_EPSILON)) break;
|
if(!lindep && (norm > PETSC_MACHINE_EPSILON)) break;
|
ierr = PetscInfo1(ip, "Orthonormalization problems adding the vector %d to the searching subspace\n", i);
|
ierr = PetscInfo1(ip, "Orthonormalization problems adding the vector %D to the searching subspace\n", i);
|
CHKERRQ(ierr);
|
CHKERRQ(ierr);
|
}
|
}
|
if(lindep || (norm < PETSC_MACHINE_EPSILON)) {
|
if(lindep || (norm < PETSC_MACHINE_EPSILON)) {
|
SETERRQ(((PetscObject)ip)->comm,1, "Error during orthonormalization of eigenvectors");
|
SETERRQ(((PetscObject)ip)->comm,1, "Error during orthonormalization of eigenvectors");
|
}
|
}
|
| Line 1280... |
Line 1280... |
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
auxR = (PetscReal*)auxS; auxS = (PetscScalar*)(auxR+2*n); size_auxS-= 2*n;
|
auxR = (PetscReal*)auxS; auxS = (PetscScalar*)(auxR+2*n); size_auxS-= 2*n;
|
if (size_auxS < 2*n)
|
if (size_auxS < 2*n)
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xTGEVC");
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xTGEVC");
|
LAPACKtgevc_(side,howmny,PETSC_NULL,&n,Sc,&n,Tc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,auxR,&info);
|
LAPACKtgevc_(side,howmny,PETSC_NULL,&n,Sc,&n,Tc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,auxR,&info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTGEVC %i",info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTGEVC %d",info);
|
#else
|
#else
|
alphar = auxS; auxS+= n; size_auxS-= n;
|
alphar = auxS; auxS+= n; size_auxS-= n;
|
alphai = auxS; auxS+= n; size_auxS-= n;
|
alphai = auxS; auxS+= n; size_auxS-= n;
|
beta = auxS; auxS+= n; size_auxS-= n;
|
beta = auxS; auxS+= n; size_auxS-= n;
|
if (doProd) {
|
if (doProd) {
|
| Line 1299... |
Line 1299... |
of T that represent complex conjugate eigenpairs to be zero */
|
of T that represent complex conjugate eigenpairs to be zero */
|
n1 = size_auxS;
|
n1 = size_auxS;
|
if (size_auxS < 8*n)
|
if (size_auxS < 8*n)
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xGGEV");
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xGGEV");
|
LAPACKggev_(pY?"V":"N",pX?"V":"N",&n,Sc,&n,Tc,&n,alphar,alphai,beta,pB,&ldpB,pA,&ldpA,auxS,&n1,&info);
|
LAPACKggev_(pY?"V":"N",pX?"V":"N",&n,Sc,&n,Tc,&n,alphar,alphai,beta,pB,&ldpB,pA,&ldpA,auxS,&n1,&info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGGEV %i",info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGGEV %d",info);
|
if (doProd) {
|
if (doProd) {
|
if (pX) {
|
if (pX) {
|
/* pX <- pX * pA */
|
/* pX <- pX * pA */
|
ierr = SlepcDenseCopy(Sc, n, pX, ldpX, n, n); CHKERRQ(ierr);
|
ierr = SlepcDenseCopy(Sc, n, pX, ldpX, n, n); CHKERRQ(ierr);
|
ierr = SlepcDenseMatProd(pX, ldpX, 0.0, 1.0,
|
ierr = SlepcDenseMatProd(pX, ldpX, 0.0, 1.0,
|
| Line 1327... |
Line 1327... |
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xTREVC");
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Insufficient work space for xTREVC");
|
LAPACKtrevc_(side,howmny,PETSC_NULL,&n,Sc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,auxR,&info);
|
LAPACKtrevc_(side,howmny,PETSC_NULL,&n,Sc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,auxR,&info);
|
#else
|
#else
|
LAPACKtrevc_(side,howmny,PETSC_NULL,&n,Sc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,&info);
|
LAPACKtrevc_(side,howmny,PETSC_NULL,&n,Sc,&n,pY,&ldpY,pX,&ldpX,&n,&nout,auxS,&info);
|
#endif
|
#endif
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTREVC %d",info);
|
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
| Line 1364... |
Line 1364... |
if (i<n-1 && S[i*ldS+i+1] != 0.0) {
|
if (i<n-1 && S[i*ldS+i+1] != 0.0) {
|
ierr = SlepcDenseCopy(Sc, 2, &S[i*ldS+i], ldS, 2, 2); CHKERRQ(ierr);
|
ierr = SlepcDenseCopy(Sc, 2, &S[i*ldS+i], ldS, 2, 2); CHKERRQ(ierr);
|
ierr = SlepcDenseCopy(Tc, 2, &T[i*ldT+i], ldT, 2, 2); CHKERRQ(ierr);
|
ierr = SlepcDenseCopy(Tc, 2, &T[i*ldT+i], ldT, 2, 2); CHKERRQ(ierr);
|
LAPACKggev_("N","N",&two,Sc,&two,Tc,&two,&eigr[i],&eigi[i],beta,
|
LAPACKggev_("N","N",&two,Sc,&two,Tc,&two,&eigr[i],&eigi[i],beta,
|
PETSC_NULL, &two,PETSC_NULL,&two,auxS,&size_auxS,&info);
|
PETSC_NULL, &two,PETSC_NULL,&two,auxS,&size_auxS,&info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGGEV %i",info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGGEV %d",info);
|
eigr[i] /= beta[0]; eigi[i] /= beta[0];
|
eigr[i] /= beta[0]; eigi[i] /= beta[0];
|
eigr[i+1]/= beta[1]; eigi[i+1]/= beta[1];
|
eigr[i+1]/= beta[1]; eigi[i+1]/= beta[1];
|
i++;
|
i++;
|
} else
|
} else
|
#endif
|
#endif
|