<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Dear All,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
I am seeking to use slepc/petsc to solve a generalised eigenvalue problem</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
that arises in a hydrodynamic stability problem. The code is a parallelisation of an existing
<br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
serial Fortran code. Before I get to grips with this target problem, I would of course like to get</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
some relevant examples working. I have installed petsc/ slepc seemingly without any problem,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
and the provided slepc example fortran program ex1f.F, which solves a regular eigenvalue</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
problem Ax = lambda x, seemed to compile and run correctly.</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
I have now written a short program to instead solve the complex generalised <br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
problem Ax = lambda B x (see below) . This code compiles and runs w/out errors but</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
for some reason hangs when calling EPSSolve - we enter EPSSolve but never leave.</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
The matrices appear to be correctly assembled -all the values are correct in the Matview</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
printout, so I am not quite sure where I have gone wrong, can anyone spot my mistake?<br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
( Note that for the actual problem I wish to solve I have already written the code to construct the matrix,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
which distributes the rows across the processes and it is fully tested and working. Hence I want to specify</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
the distribution of rows and not leave it up to a PETS_DECIDE .)</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
I would be very grateful if someone can point out what is wrong with this small example code (see below),</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Many thanks,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Dan.</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
! this program MUST be run with NGLOBAL = 10, MY_NO_ROWS = 5
<div>! and two MPI processes since row distribution is hard-baked into code </div>
<div>!</div>
<div> program main</div>
<div>#include <slepc/finclude/slepceps.h></div>
<div> use slepceps</div>
<div> implicit none</div>
<div><br>
</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! Declarations</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>!</div>
<div>! Variables:</div>
<div>! A , B double complex operator matrices</div>
<div>! eps eigenproblem solver context</div>
<div><br>
</div>
<div> Mat A,B</div>
<div> EPS eps</div>
<div> EPSType tname</div>
<div> PetscReal tol, error</div>
<div> PetscScalar kr, ki</div>
<div> Vec xr, xi</div>
<div> PetscInt NGLOBAL , MY_NO_ROWS, NL3, owner</div>
<div> PetscInt nev, maxit, its, nconv</div>
<div> PetscInt i,j,ic,jc</div>
<div> PetscReal reala, imaga, realb, imagb, di, dj</div>
<div> PetscScalar a_entry, b_entry</div>
<div> PetscMPIInt rank</div>
<div> PetscErrorCode ierr</div>
<div> PetscInt,parameter :: zero = 0, one = 1, two = 2, three = 3 </div>
<div> </div>
<div> PetscInt M1, N1, mm1, nn1</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! Beginning of program</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div><br>
</div>
<div> call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)</div>
<div> call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)</div>
<div>! make sure you set NGLOBAL = 10, MY_NO_ROWS = 5 and run with two processes</div>
<div> NGLOBAL = 10</div>
<div> MY_NO_ROWS = 5</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! Compute the operator matrices that define the eigensystem, Ax=kBx</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>!!!!!!!!! Setup A matrix</div>
<div><br>
</div>
<div> call MatCreate(PETSC_COMM_WORLD,A,ierr) </div>
<div> call MatSetsizes(A,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) </div>
<div> call MatSetFromOptions(A,ierr)</div>
<div> call MatGetSize(A,M1,N1,ierr)</div>
<div> write(*,*)'Rank [',rank,']: global size of A is ',M1, N1</div>
<div> call MatGetLocalSize(A,mm1,nn1,ierr)</div>
<div> write(*,*)'Rank [',rank,']: my local size of A is ',mm1, nn1</div>
<div> call MatMPIAIJSetPreallocation(A,three, PETSC_NULL_INTEGER,one, &</div>
<div> & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation</div>
<div><br>
</div>
<div>!!!!!!!!! Setup B matrix</div>
<div> call MatCreate(PETSC_COMM_WORLD,B,ierr) </div>
<div> call MatSetsizes(B,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) </div>
<div> call MatSetFromOptions(B,ierr)</div>
<div> call MatGetSize(B,M1,N1,ierr)</div>
<div> write(*,*)'Rank [',rank,']: global size of B is ',M1, N1</div>
<div> call MatGetLocalSize(B,mm1,nn1,ierr)</div>
<div> write(*,*)'Rank [',rank,']: my local size of B is ',mm1, nn1</div>
<div> </div>
<div> call MatMPIAIJSetPreallocation(B,three, PETSC_NULL_INTEGER,one, &</div>
<div> & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation</div>
<div> </div>
<div>! initalise</div>
<div> call MatZeroEntries(A,ierr)</div>
<div> call MatZeroEntries(B,ierr)</div>
<div> </div>
<div>! Fill in values of A, B and assemble matrices</div>
<div>! Both matrices are tridiagonal with</div>
<div>! Aij = cmplx( (i-j)**2, (i+j)**2)</div>
<div>! Bij = cmplx( ij/i + j, (i/j)**2)</div>
<div>! (numbering from 1 )</div>
<div> </div>
<div> do i = 1, NGLOBAL</div>
<div> ! a rather crude way to distribute rows </div>
<div> if (i < 6) owner = 0</div>
<div> if (i >= 6) owner = 1</div>
<div> if (rank /= owner) cycle</div>
<div> do j = 1, NGLOBAL</div>
<div> if ( abs(i-j) < 2 ) then</div>
<div> write(*,*)rank,' : Setting ',i,j</div>
<div> di = dble(i) ; dj = dble(j)</div>
<div> </div>
<div> reala = (di - dj)**2 ; imaga = (di + dj)**2</div>
<div> a_entry = dcmplx(reala, imaga)</div>
<div> realb = (di*dj)/(di + dj) ; imagb = di**2/dj**2</div>
<div> b_entry = dcmplx(realb, imagb) <br>
</div>
<div> </div>
<div> ic = i -1 ; jc = j-1 ! convert to C indexing</div>
<div> call MatSetValue(A, ic, jc, a_entry, ADD_VALUES,ierr)</div>
<div> call MatSetValue(B, ic, jc, b_entry, ADD_VALUES,ierr)</div>
<div> endif</div>
<div> enddo</div>
<div> enddo</div>
<div><br>
</div>
<div> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) </div>
<div> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)</div>
<div> call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) </div>
<div> call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)</div>
<div> </div>
<div>! Check matrices </div>
<div> write(*,*)'A matrix ... '</div>
<div> call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)</div>
<div> write(*,*)'B matrix ... '</div>
<div> call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)</div>
<div> </div>
<div> call MatCreateVecs(A,PETSC_NULL_VEC,xr)</div>
<div> call MatCreateVecs(A, PETSC_NULL_VEC,xi)</div>
<div><br>
</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! Create the eigensolver and display info</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! ** Create eigensolver context</div>
<div> call EPSCreate(PETSC_COMM_WORLD,eps,ierr)</div>
<div><br>
</div>
<div>! ** Set operators.for general problem Ax = lambda B x</div>
<div> call EPSSetOperators(eps,A, B, ierr) </div>
<div> call EPSSetProblemType(eps,EPS_GNHEP,ierr)</div>
<div><br>
</div>
<div>! ** Set solver parameters at runtime</div>
<div> call EPSSetFromOptions(eps,ierr)</div>
<div> </div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div>! Solve the eigensystem</div>
<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>
<div><br>
</div>
<div> write(*,*)'rank',rank, 'entering solver ...'</div>
<div> call EPSSolve(eps,ierr)</div>
<div> </div>
<div>! ** Free work space</div>
<div> call EPSDestroy(eps,ierr)</div>
<div> call MatDestroy(A,ierr)</div>
<div> call MatDestroy(B,ierr)</div>
<div> </div>
<div> call VecDestroy(xr,ierr)</div>
<div> call VecDestroy(xi,ierr)</div>
<div><br>
</div>
<div> call SlepcFinalize(ierr)</div>
<div> end</div>
<br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
</body>
</html>