| Line 63... |
Line 63... |
ierr = QEPSetType(qep,QEPLINEAR);CHKERRQ(ierr);
|
ierr = QEPSetType(qep,QEPLINEAR);CHKERRQ(ierr);
|
}
|
}
|
|
|
/* Check matrices */
|
/* Check matrices */
|
if (!qep->M || !qep->C || !qep->K)
|
if (!qep->M || !qep->C || !qep->K)
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "QEPSetOperators must be called first");
|
SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_WRONGSTATE, "QEPSetOperators must be called first");
|
|
|
/* Set problem dimensions */
|
/* Set problem dimensions */
|
ierr = MatGetSize(qep->M,&qep->n,PETSC_NULL);CHKERRQ(ierr);
|
ierr = MatGetSize(qep->M,&qep->n,PETSC_NULL);CHKERRQ(ierr);
|
ierr = MatGetLocalSize(qep->M,&qep->nloc,PETSC_NULL);CHKERRQ(ierr);
|
ierr = MatGetLocalSize(qep->M,&qep->nloc,PETSC_NULL);CHKERRQ(ierr);
|
|
|
| Line 94... |
Line 94... |
|
|
/* Call specific solver setup */
|
/* Call specific solver setup */
|
ierr = (*qep->ops->setup)(qep);CHKERRQ(ierr);
|
ierr = (*qep->ops->setup)(qep);CHKERRQ(ierr);
|
|
|
if (qep->ncv > 2*qep->n)
|
if (qep->ncv > 2*qep->n)
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"ncv must be twice the problem size at most");
|
SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_OUTOFRANGE,"ncv must be twice the problem size at most");
|
if (qep->nev > qep->ncv)
|
if (qep->nev > qep->ncv)
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"nev bigger than ncv");
|
SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nev bigger than ncv");
|
|
|
/* Free memory for previous solution */
|
/* Free memory for previous solution */
|
if (qep->eigr) {
|
if (qep->eigr) {
|
ierr = PetscFree(qep->eigr);CHKERRQ(ierr);
|
ierr = PetscFree(qep->eigr);CHKERRQ(ierr);
|
ierr = PetscFree(qep->eigi);CHKERRQ(ierr);
|
ierr = PetscFree(qep->eigi);CHKERRQ(ierr);
|
| Line 126... |
Line 126... |
}
|
}
|
|
|
/* process initial vectors */
|
/* process initial vectors */
|
if (qep->nini<0) {
|
if (qep->nini<0) {
|
qep->nini = -qep->nini;
|
qep->nini = -qep->nini;
|
if (qep->nini>qep->ncv) SETERRQ(1,"The number of initial vectors is larger than ncv")
|
if (qep->nini>qep->ncv) SETERRQ(((PetscObject)qep)->comm,1,"The number of initial vectors is larger than ncv")
|
k = 0;
|
k = 0;
|
for (i=0;i<qep->nini;i++) {
|
for (i=0;i<qep->nini;i++) {
|
ierr = VecCopy(qep->IS[i],qep->V[k]);CHKERRQ(ierr);
|
ierr = VecCopy(qep->IS[i],qep->V[k]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->IS[i]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->IS[i]);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(qep->ip,0,PETSC_NULL,k,PETSC_NULL,qep->V,qep->V[k],PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(qep->ip,0,PETSC_NULL,k,PETSC_NULL,qep->V,qep->V[k],PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
|
| Line 145... |
Line 145... |
}
|
}
|
if (qep->ninil<0) {
|
if (qep->ninil<0) {
|
if (!qep->leftvecs) PetscInfo(qep,"Ignoring initial left vectors\n");
|
if (!qep->leftvecs) PetscInfo(qep,"Ignoring initial left vectors\n");
|
else {
|
else {
|
qep->ninil = -qep->ninil;
|
qep->ninil = -qep->ninil;
|
if (qep->ninil>qep->ncv) SETERRQ(1,"The number of initial left vectors is larger than ncv")
|
if (qep->ninil>qep->ncv) SETERRQ(((PetscObject)qep)->comm,1,"The number of initial left vectors is larger than ncv")
|
k = 0;
|
k = 0;
|
for (i=0;i<qep->ninil;i++) {
|
for (i=0;i<qep->ninil;i++) {
|
ierr = VecCopy(qep->ISL[i],qep->W[k]);CHKERRQ(ierr);
|
ierr = VecCopy(qep->ISL[i],qep->W[k]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->ISL[i]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->ISL[i]);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(qep->ip,0,PETSC_NULL,k,PETSC_NULL,qep->W,qep->W[k],PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
|
ierr = IPOrthogonalize(qep->ip,0,PETSC_NULL,k,PETSC_NULL,qep->W,qep->W[k],PETSC_NULL,&norm,&lindep);CHKERRQ(ierr);
|
| Line 204... |
Line 204... |
PetscCheckSameComm(qep,1,C,3);
|
PetscCheckSameComm(qep,1,C,3);
|
PetscCheckSameComm(qep,1,K,4);
|
PetscCheckSameComm(qep,1,K,4);
|
|
|
/* Check for square matrices */
|
/* Check for square matrices */
|
ierr = MatGetSize(M,&m,&n);CHKERRQ(ierr);
|
ierr = MatGetSize(M,&m,&n);CHKERRQ(ierr);
|
if (m!=n) { SETERRQ(1,"M is a non-square matrix"); }
|
if (m!=n) { SETERRQ(((PetscObject)qep)->comm,1,"M is a non-square matrix"); }
|
m0=m;
|
m0=m;
|
ierr = MatGetSize(C,&m,&n);CHKERRQ(ierr);
|
ierr = MatGetSize(C,&m,&n);CHKERRQ(ierr);
|
if (m!=n) { SETERRQ(1,"C is a non-square matrix"); }
|
if (m!=n) { SETERRQ(((PetscObject)qep)->comm,1,"C is a non-square matrix"); }
|
if (m!=m0) { SETERRQ(1,"Dimensions of M and C do not match"); }
|
if (m!=m0) { SETERRQ(((PetscObject)qep)->comm,1,"Dimensions of M and C do not match"); }
|
ierr = MatGetSize(K,&m,&n);CHKERRQ(ierr);
|
ierr = MatGetSize(K,&m,&n);CHKERRQ(ierr);
|
if (m!=n) { SETERRQ(1,"K is a non-square matrix"); }
|
if (m!=n) { SETERRQ(((PetscObject)qep)->comm,1,"K is a non-square matrix"); }
|
if (m!=m0) { SETERRQ(1,"Dimensions of M and K do not match"); }
|
if (m!=m0) { SETERRQ(((PetscObject)qep)->comm,1,"Dimensions of M and K do not match"); }
|
|
|
/* Store a copy of the matrices */
|
/* Store a copy of the matrices */
|
ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr);
|
ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr);
|
if (qep->M) {
|
if (qep->M) {
|
ierr = MatDestroy(qep->M);CHKERRQ(ierr);
|
ierr = MatDestroy(qep->M);CHKERRQ(ierr);
|
| Line 300... |
Line 300... |
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
|
PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
|
if (n<0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
|
if (n<0) SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
|
|
|
/* free previous non-processed vectors */
|
/* free previous non-processed vectors */
|
if (qep->nini<0) {
|
if (qep->nini<0) {
|
for (i=0;i<-qep->nini;i++) {
|
for (i=0;i<-qep->nini;i++) {
|
ierr = VecDestroy(qep->IS[i]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->IS[i]);CHKERRQ(ierr);
|
| Line 360... |
Line 360... |
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i;
|
PetscInt i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
|
PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
|
if (n<0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
|
if (n<0) SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
|
|
|
/* free previous non-processed vectors */
|
/* free previous non-processed vectors */
|
if (qep->ninil<0) {
|
if (qep->ninil<0) {
|
for (i=0;i<-qep->ninil;i++) {
|
for (i=0;i<-qep->ninil;i++) {
|
ierr = VecDestroy(qep->ISL[i]);CHKERRQ(ierr);
|
ierr = VecDestroy(qep->ISL[i]);CHKERRQ(ierr);
|