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

[Thomas Hisch]

On Wed, Aug 8, 2012 at 1:43 PM, Jose E. Roman <jroman at dsic.upv.es> wrote:
>
> 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

Thx. This does wonders!! The Setup time went down from 18sec to 1 sec
and the Solve time from 65sec to 10 sec.

Regards
Thomas

>