Discussion:
[petsc-users] number of converged eigenpairs
Soumya Mukherjee
2015-09-24 16:01:48 UTC
Permalink
Hi.

The following code generates nconv numbers of eigenpairs. This number, for
my case is so few. And I am afraid, there are more important eigenmodes. Is
it possible to generate all (or more) eigenpairs?

Also, I was wondering, which eigenvalues/ eigenvectors does it calculate in
this code, whose number is given by nconv?

MatCreateSeqAIJ( MPI_COMM_SELF, Natom * Norbit, Natom * Norbit,
Nneighbor * Norbit, PETSC_NULL, &H );
MatSetFromOptions( H );

MatCreateSeqAIJ( MPI_COMM_SELF, Natom * Norbit, Natom * Norbit,
Nneighbor * Norbit, PETSC_NULL, &S );
MatSetFromOptions( S );

Hamiltonian ( );
// MatView(H,PETSC_VIEWER_STDOUT_WORLD);
// VecView(H0,PETSC_VIEWER_STDOUT_WORLD);

//-------------------------------------------------------------------------------------------------------
/*Setup the Ax=cBx generalized eigenproblem*/

//-------------------------------------------------------------------------------------------------------
EPSCreate( PETSC_COMM_WORLD, &eps );
EPSSetOperators( eps, H, S);
MatGetVecs(H,NULL,&xr);
MatGetVecs(H,NULL,&xi);

EPSGetST( eps, &st);
STSetType(st,STCAYLEY);
STCayleySetAntishift(st,1);

EPSSetDimensions( eps, nev, ncv, PETSC_DECIDE);
EPSSetTarget( eps,5.0);
EPSSetFromOptions( eps );
EPSSolve( eps );
EPSPrintSolution(eps,PETSC_NULL);
EPSGetConverged( eps, &nconv );

eigenval = new complex<double>[ nconv ];
eigenvec = new complex<double>*[ nconv ];
for (int i = 0; i < nconv; i ++ ) {
eigenvec[ i ] = new complex<double>[ Norbit * Natom ];
}

for ( jjj = 0; jjj < nconv; jjj ++ ) {
/* Get converged eigenpairs: j-th eigenvalue is stored in kr (real
part) and ki (imaginary part) */
EPSGetEigenpair(eps, jjj, &kr, &ki, xr, xi);
eigenval[ jjj ] = complex
<double>(PetscRealPart(kr),PetscImaginaryPart(kr));

VecGetArray(xr, &Xr);
for ( int i = 0; i < Norbit * Natom; i ++ ) {
eigenvec[ jjj ][ i ] = complex
<double>(PetscRealPart(Xr[i]),PetscImaginaryPart(Xr[i]));
}
VecRestoreArray(xr, &Xr);
}
//--------------------------
/*finalizing*/
//--------------------------
outputs( true );
VecDestroy(&xr);
VecDestroy(&xi);
EPSDestroy(&eps);
MatDestroy(&H);
MatDestroy(&S);
SlepcFinalize();
PetscFinalize();
return 0;
}

Regards,
Soumya
Jose E. Roman
2015-09-24 20:50:25 UTC
Permalink
By default only 1 eigenpair is computed. Use EPSSetDimensions() to specify how many eigenvalues you want. This is explained in the manual in basic examples.
http://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetDimensions.html

SLEPc is not intended for computing all eigenvalues. If you need all eigenvalues you should use another software.

In your code, the computed eigenvalues are those closest to the target.

Jose
Hi.
The following code generates nconv numbers of eigenpairs. This number, for my case is so few. And I am afraid, there are more important eigenmodes. Is it possible to generate all (or more) eigenpairs?
Also, I was wondering, which eigenvalues/ eigenvectors does it calculate in this code, whose number is given by nconv?
MatCreateSeqAIJ( MPI_COMM_SELF, Natom * Norbit, Natom * Norbit, Nneighbor * Norbit, PETSC_NULL, &H );
MatSetFromOptions( H );
MatCreateSeqAIJ( MPI_COMM_SELF, Natom * Norbit, Natom * Norbit, Nneighbor * Norbit, PETSC_NULL, &S );
MatSetFromOptions( S );
Hamiltonian ( );
// MatView(H,PETSC_VIEWER_STDOUT_WORLD);
// VecView(H0,PETSC_VIEWER_STDOUT_WORLD);
//-------------------------------------------------------------------------------------------------------
/*Setup the Ax=cBx generalized eigenproblem*/
//-------------------------------------------------------------------------------------------------------
EPSCreate( PETSC_COMM_WORLD, &eps );
EPSSetOperators( eps, H, S);
MatGetVecs(H,NULL,&xr);
MatGetVecs(H,NULL,&xi);
EPSGetST( eps, &st);
STSetType(st,STCAYLEY);
STCayleySetAntishift(st,1);
EPSSetDimensions( eps, nev, ncv, PETSC_DECIDE);
EPSSetTarget( eps,5.0);
EPSSetFromOptions( eps );
EPSSolve( eps );
EPSPrintSolution(eps,PETSC_NULL);
EPSGetConverged( eps, &nconv );
eigenval = new complex<double>[ nconv ];
eigenvec = new complex<double>*[ nconv ];
for (int i = 0; i < nconv; i ++ ) {
eigenvec[ i ] = new complex<double>[ Norbit * Natom ];
}
for ( jjj = 0; jjj < nconv; jjj ++ ) {
/* Get converged eigenpairs: j-th eigenvalue is stored in kr (real part) and ki (imaginary part) */
EPSGetEigenpair(eps, jjj, &kr, &ki, xr, xi);
eigenval[ jjj ] = complex <double>(PetscRealPart(kr),PetscImaginaryPart(kr));
VecGetArray(xr, &Xr);
for ( int i = 0; i < Norbit * Natom; i ++ ) {
eigenvec[ jjj ][ i ] = complex <double>(PetscRealPart(Xr[i]),PetscImaginaryPart(Xr[i]));
}
VecRestoreArray(xr, &Xr);
}
//--------------------------
/*finalizing*/
//--------------------------
outputs( true );
VecDestroy(&xr);
VecDestroy(&xi);
EPSDestroy(&eps);
MatDestroy(&H);
MatDestroy(&S);
SlepcFinalize();
PetscFinalize();
return 0;
}
Regards,
Soumya
Loading...