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

[Thomas Hisch]

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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: out_nodebugging_full.log
Type: application/octet-stream
Size: 11312 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120808/68a719bb/attachment.obj>