| - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| */ |
| #include "src/st/stimpl.h" |
| #include "src/eps/epsimpl.h" |
| #include "src/contrib/blopex/petsc-interface/petsc-interface.h" |
| #include "lobpcg.h" |
| KSP ksp; |
| } EPS_BLOPEX; |
| #undef __FUNCT__ |
| #define __FUNCT__ "Precond_FnSingleVector" |
| static void Precond_FnSingleVector(void *data,void *x,void *y) |
| PetscErrorCode ierr; |
| EPS eps = (EPS)data; |
| EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data; |
| PetscInt lits; |
| PetscFunctionBegin; |
| ierr = KSPSolve(blopex->ksp,(Vec)x,(Vec)y); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| ierr = KSPGetIterationNumber(blopex->ksp, &lits); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| eps->OP->lineariterations+= lits; |
| PetscFunctionReturnVoid(); |
| } |
| { |
| PetscErrorCode ierr; |
| EPS eps = (EPS)data; |
| Mat A; |
| PetscFunctionBegin; |
| ierr = STApply(eps->OP,(Vec)x,(Vec)y); |
| ierr = IPOrthogonalize(eps->ip,eps->nds,PETSC_NULL,eps->DS,(Vec)y,PETSC_NULL,PETSC_NULL,PETSC_NULL,eps->work[0]); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| ierr = STGetOperators(eps->OP,&A,PETSC_NULL); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| ierr = MatMult(A,(Vec)x,(Vec)y); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| //ierr = IPOrthogonalize(eps->ip,eps->nds,PETSC_NULL,eps->DS,(Vec)y,PETSC_NULL,PETSC_NULL,PETSC_NULL,eps->work[0]); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| PetscFunctionReturnVoid(); |
| } |
| } |
| #undef __FUNCT__ |
| #define __FUNCT__ "OperatorBSingleVector" |
| static void OperatorBSingleVector(void *data,void *x,void *y) |
| { |
| PetscErrorCode ierr; |
| EPS eps = (EPS)data; |
| Mat B; |
| PetscFunctionBegin; |
| ierr = STGetOperators(eps->OP,PETSC_NULL,&B); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| ierr = MatMult(B,(Vec)x,(Vec)y); CHKERRABORT(PETSC_COMM_WORLD,ierr); |
| PetscFunctionReturnVoid(); |
| } |
| #undef __FUNCT__ |
| #define __FUNCT__ "OperatorBMultiVector" |
| static void OperatorBMultiVector(void *data,void *x,void *y) |
| { |
| EPS eps = (EPS)data; |
| EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data; |
| PetscFunctionBegin; |
| blopex->ii.Eval(OperatorBSingleVector,data,x,y); |
| PetscFunctionReturnVoid(); |
| } |
| #undef __FUNCT__ |
| #define __FUNCT__ "EPSSetUp_BLOPEX" |
| PetscErrorCode EPSSetUp_BLOPEX(EPS eps) |
| { |
| PetscErrorCode ierr; |
| EPS_BLOPEX *blopex = (EPS_BLOPEX *)eps->data; |
| Mat A; |
| Mat A,B; |
| int N; |
| PetscTruth isShift; |
| PetscFunctionBegin; |
| if (!eps->ishermitian || eps->isgeneralized) { |
| if (!eps->ishermitian) { |
| SETERRQ(PETSC_ERR_SUP,"blopex only works for standard symmetric problems"); |
| } |
| ierr = PetscTypeCompare((PetscObject)eps->OP,STSHIFT,&isShift);CHKERRQ(ierr); |
| if (eps->which!=EPS_SMALLEST_REAL) { |
| SETERRQ(1,"Wrong value of eps->which"); |
| } |
| ierr = STGetOperators(eps->OP,&A,PETSC_NULL);CHKERRQ(ierr); |
| ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr); |
| ierr = KSPSetOperators(blopex->ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); |
| ierr = KSPSetUp(blopex->ksp);CHKERRQ(ierr); |
| KSPView(blopex->ksp, PETSC_VIEWER_STDOUT_WORLD); |
| ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr); |
| eps->ncv = eps->nev = PetscMin(eps->nev,N); |
| PetscFunctionBegin; |
| info = lobpcg_solve(blopex->eigenvectors,eps,OperatorAMultiVector, |
| PETSC_NULL,PETSC_NULL,eps,Precond_FnMultiVector,NULL, |
| blopex->blap_fn,blopex->tol,eps->max_it,0,&eps->its, |
| eps->isgeneralized?eps:PETSC_NULL,eps->isgeneralized?OperatorBMultiVector:PETSC_NULL, |
| eps,Precond_FnMultiVector,NULL, |
| blopex->blap_fn,blopex->tol,eps->max_it,1,&eps->its, |
| eps->eigr,PETSC_NULL,0,eps->errest,PETSC_NULL,0); |
| if (info>0) SETERRQ1(PETSC_ERR_LIB,"Error in blopex (code=%d)",info); |
| eps->ops->backtransform = EPSBackTransform_Default; |
| eps->ops->computevectors = EPSComputeVectors_Default; |
| eps->which = EPS_SMALLEST_REAL; |
| /* TODO: meter opciones predeterminadas: ST sinvert, KSP preonly */ |
| PetscFunctionReturn(0); |
| } |
| EXTERN_C_END |