/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
This file is part of SLEPc.
SLEPc is free software: you can redistribute it and/or modify it under the
terms of version 3 of the GNU Lesser General Public License as published by
the Free Software Foundation.
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
static char help[] = "Tests B-orthonormality of eigenvectors in a GHEP problem.\n\n";
#include <slepceps.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv )
{
Mat A,B; /* matrices */
EPS eps; /* eigenproblem solver context */
Vec *X,v;
PetscReal lev,tol=1000*PETSC_MACHINE_EPSILON;
PetscInt N,n=45,m,Istart,Iend,II,i,j,nconv;
PetscBool flag;
PetscErrorCode ierr;
SlepcInitialize(&argc,&argv,(char*)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);CHKERRQ(ierr);
if(!flag) m=n;
N = n*m;
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized Symmetric Eigenproblem, N=%D (%Dx%D grid)\n\n",N,n,m);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrices that define the eigensystem, Ax=kBx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(B);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
for( II=Istart; II<Iend; II++ ) {
i = II/n; j = II-i*n;
if(i>0) { ierr = MatSetValue(A,II,II-n,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if(i<m-1) { ierr = MatSetValue(A,II,II+n,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if(j>0) { ierr = MatSetValue(A,II,II-1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if(j<n-1) { ierr = MatSetValue(A,II,II+1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
ierr = MatSetValue(A,II,II,4.0,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValue(B,II,II,2.0/PetscLogScalar(II+2),INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatGetVecs(B,&v,PETSC_NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
ierr = EPSSetOperators(eps,A,B);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps,EPS_GHEP);CHKERRQ(ierr);
ierr = EPSSetTolerances(eps,tol,PETSC_DECIDE);CHKERRQ(ierr);
ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the eigensystem
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = EPSSolve(eps);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = EPSPrintSolution(eps,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
if (nconv>0) {
ierr = VecDuplicateVecs(v,nconv,&X);CHKERRQ(ierr);
for( i=0; i<nconv; i++ ) {
ierr = EPSGetEigenvector(eps,i,X[i],PETSC_NULL);CHKERRQ(ierr);
}
ierr = SlepcCheckOrthogonality(X,nconv,PETSC_NULL,nconv,B,&lev);CHKERRQ(ierr);
if (lev<tol) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n");CHKERRQ(ierr);
} else {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %G\n",lev);CHKERRQ(ierr);
}
ierr = VecDestroyVecs(nconv,&X);CHKERRQ(ierr);
}
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = MatDestroy(&B);CHKERRQ(ierr);
ierr = VecDestroy(&v);CHKERRQ(ierr);
ierr = SlepcFinalize();CHKERRQ(ierr);
return 0;
}