[petsc-users] Speedup calculation of a high number of generalized eigenvalue problems

[Jose E. Roman]

El 08/08/2012, a las 11:52, Thomas Hisch escribió:

> Hello list,
> 
> What can be done to speedup the calculation of the eigenvals and
> eigenvecs of many (~1000)!! generalized EVP, where the operator A is
> constant for all GEVPs and only B is different (but the nonzero
> pattern, size, type,.. are the same for all B). For now I create an
> EPS object for every pair (A,B) like in the following code: (Note: B
> is updated in other functions)
> 
> 
> void SlepcEigenSolver::Solve(ngbla::Matrix<ngla::Complex> &evecs,
>                             ngbla::Vector<ngla::Complex> &evals)
> {
>  ....
> 
>  Vec xev;
>  MatGetVecs(A_,PETSC_NULL,&xev);
>  PetscPrintf(PETSC_COMM_WORLD,"Create EPS object and solve eigensystem\n");
>  EPSCreate(PETSC_COMM_WORLD, &eps_);
>  EPSSetOperators(eps_,A_,B_);
>  EPSSetDimensions(eps_, GetNev(), PETSC_DECIDE, PETSC_DECIDE);
>  EPSSetFromOptions(eps_);
> 
>  PetscScalar target = GetShift();
>  target *= target;
> 
>  EPSSetTarget(eps_,target);
>  EPSSetWhichEigenpairs(eps_,EPS_TARGET_MAGNITUDE);
> 
>  ST myst;
>  EPSGetST(eps_,&myst);
>  STSetType(myst,STSINVERT);
> 
>  PetscLogDouble t1,t2,t3;
>  PetscGetTime(&t1);
>  EPSSetUp(eps_);
>  PetscGetTime(&t2);
>  EPSSolve(eps_);
>  PetscGetTime(&t3);
>  PetscPrintf(PETSC_COMM_WORLD," -> EPS Solving Done\n\n");
>  PetscPrintf(PETSC_COMM_WORLD,"EPS: Elapsed time: %G (setup), %G
> (solve)\n",t2-t1,t3-t2);
> 
>  PetscInt its;
>  EPSGetIterationNumber(eps_,&its);
> 
>  PetscPrintf(PETSC_COMM_WORLD, "number of iterations of the method %D\n", its);
> 
>  PetscInt nconv;
>  EPSGetConverged(eps_,&nconv);
>  PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
> 
>  if (nconv>0) {
>    evals.SetSize(nconv);
>    if(!rank_) {
>      evecs.SetSize(matsize,nconv);
>    }
> 
>    PetscScalar lambda; // = k
>    PetscReal error,re,im;
> 
>    for (PetscInt i=0;i<nconv;i++) {
>      // Get converged eigenpairs: i-th eigenvalue is stored in kr
> (real part) and
>      // ki (imaginary part)
>      EPSGetEigenpair(eps_,i,&lambda,PETSC_NULL,xev,PETSC_NULL);
> 
>      PetscScalar* tmpptr;
> 
>      // begin ugly inefficient(?) eigenvector gather code
>      {
>        VecScatter ctx;
>        Vec vout;
>        VecScatterCreateToZero(xev,&ctx,&vout);
>        // scatter as many times as you need
>        VecScatterBegin(ctx,xev,vout,INSERT_VALUES,SCATTER_FORWARD);
>        VecScatterEnd(ctx,xev,vout,INSERT_VALUES,SCATTER_FORWARD);
> 
>        int locsize;
>        VecGetLocalSize(vout,&locsize);
>        VecGetArray(vout,&tmpptr);
>        if(!rank_) {
>          const ngla::Complex *evecptr = static_cast<ngbla::Complex*>(tmpptr);
> 
>          for(size_t j = 0; j < matsize; j++) {
>            evecs(j,i) = evecptr[j];
>          }
>        }
> 
>        VecRestoreArray(xev, &tmpptr);
>        // destroy scatter context and local vector when no longer needed
>        VecScatterDestroy(&ctx);
>        VecDestroy(&vout);
>      }
>      /* end ugly inefficient(?) vector gather code */
> 
>      evals(i) = static_cast<ngla::Complex>(lambda);
>    }
>    PetscPrintf(PETSC_COMM_WORLD,"\n");
>  }
> 
>  VecDestroy(&xev);
>  EPSDestroy(&eps_);
> 
>  if(!nconv)
>    return;
> 
>  if(rank_)
>    evecs.SetSize(matsize,nconv);
> 
>  //every process already owns all the eigenvalues
>  //but not all eigenstates -> only bcast all eigenstates
>  MPI_Bcast(&evecs(0), nconv*matsize, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
> }
> 
> Here are some facts about my EVP:
> The dimension of the Operators A & B is ~90000.
> For nev=1 the setup time (EPSSetUp) is 18sec, and the solve time
> 11sec. (on 16 nodes)
>  -> number of iterations of the method 2
>      number of converged eigenpairs: 6
> 
> For nev=150 the setup time (EPSSetUp) is 18sec, and the solve time
> 65sec. (on 16 nodes)
>  -> number of iterations of the method 1
>      number of converged eigenpairs: 155
>  (The output of -log_summary is attached)
> 
> 
> Do you think that there is any benefit in creating the EPS object at
> the beginning of my program and only call EPSSetOperators in my Solve
> function?
> 
> Regards
> Thomas
> <out_nodebugging_full.log>

Yes, you can get some modest benefit. I would strongly suggest using a parallel linear solver such as MUMPS. See section 3.4.1 of SLEPc's manual.
Jose