[petsc-users] Newbie question: Something is wrong with this Slepc simple example for generalised eigenvalue problem

[dazza simplythebest]

Dear All,
            I am seeking to use slepc/petsc to solve a generalised eigenvalue problem
that arises in a  hydrodynamic stability problem. The code is a parallelisation of an existing
serial Fortran code. Before I get to grips with this target problem, I would of course like to get
some relevant examples working. I have installed petsc/ slepc seemingly without any problem,
and the provided slepc example fortran program ex1f.F, which solves a regular eigenvalue
problem Ax = lambda x, seemed to compile and run correctly.
I have now written a short program to instead solve the complex generalised
problem Ax = lambda B x (see below) . This code compiles and runs w/out errors but
for some reason hangs when calling EPSSolve - we enter EPSSolve but never leave.
The matrices appear to be correctly assembled -all the values are correct in the Matview
printout, so I am not quite sure where I have gone wrong, can anyone spot my mistake?
 ( Note that for the actual problem I wish to solve I have already written the code to construct the matrix,
which distributes the rows across the processes and it is fully tested and working. Hence I want to specify
 the distribution of rows and not leave it up to a PETS_DECIDE .)
I would be very grateful if someone can point out what is wrong with this small example code (see below),
   Many thanks,
                             Dan.

!  this program MUST be run with NGLOBAL = 10, MY_NO_ROWS = 5
!  and two MPI processes since row distribution is hard-baked into code
!
      program main
#include <slepc/finclude/slepceps.h>
      use slepceps
      implicit none

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Declarations
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!  Variables:
!     A , B    double complex operator matrices
!     eps   eigenproblem solver context

      Mat            A,B
      EPS            eps
      EPSType        tname
      PetscReal      tol, error
      PetscScalar    kr, ki
      Vec            xr, xi
      PetscInt       NGLOBAL ,  MY_NO_ROWS, NL3, owner
      PetscInt       nev, maxit, its, nconv
      PetscInt       i,j,ic,jc
      PetscReal      reala, imaga, realb, imagb, di, dj
      PetscScalar    a_entry, b_entry
      PetscMPIInt    rank
      PetscErrorCode ierr
      PetscInt,parameter :: zero = 0, one = 1, two = 2, three = 3

      PetscInt   M1, N1, mm1, nn1
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Beginning of program
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
      call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
!     make sure you set NGLOBAL = 10, MY_NO_ROWS = 5 and run with two processes
      NGLOBAL = 10
      MY_NO_ROWS = 5
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Compute the operator matrices that define the eigensystem, Ax=kBx
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!!!!!!!!!  Setup  A matrix

      call MatCreate(PETSC_COMM_WORLD,A,ierr)
      call MatSetsizes(A,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr)
      call MatSetFromOptions(A,ierr)
      call MatGetSize(A,M1,N1,ierr)
      write(*,*)'Rank [',rank,']: global size of A is ',M1, N1
      call MatGetLocalSize(A,mm1,nn1,ierr)
      write(*,*)'Rank [',rank,']: my local size of A is ',mm1, nn1
      call MatMPIAIJSetPreallocation(A,three, PETSC_NULL_INTEGER,one,    &
     &      PETSC_NULL_INTEGER,ierr)  !parallel (MPI) allocation

!!!!!!!!!  Setup  B matrix
      call MatCreate(PETSC_COMM_WORLD,B,ierr)
      call MatSetsizes(B,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr)
      call MatSetFromOptions(B,ierr)
      call MatGetSize(B,M1,N1,ierr)
       write(*,*)'Rank [',rank,']: global size of B is ',M1, N1
      call MatGetLocalSize(B,mm1,nn1,ierr)
      write(*,*)'Rank [',rank,']: my local size of B is ',mm1, nn1

      call MatMPIAIJSetPreallocation(B,three, PETSC_NULL_INTEGER,one,    &
     &      PETSC_NULL_INTEGER,ierr)  !parallel (MPI) allocation

! initalise
      call MatZeroEntries(A,ierr)
      call MatZeroEntries(B,ierr)

!  Fill in values of A, B and assemble matrices
!  Both matrices are tridiagonal with
!  Aij = cmplx( (i-j)**2, (i+j)**2)
!  Bij = cmplx( ij/i + j, (i/j)**2)
!  (numbering from 1 )

      do i = 1, NGLOBAL
        ! a rather crude way to distribute rows
        if (i < 6) owner = 0
        if (i >= 6) owner = 1
        if (rank /= owner) cycle
        do j = 1, NGLOBAL
            if ( abs(i-j) < 2 ) then
                 write(*,*)rank,' : Setting ',i,j
                 di = dble(i)         ; dj = dble(j)

                 reala = (di - dj)**2 ; imaga = (di + dj)**2
                 a_entry = dcmplx(reala, imaga)
                 realb = (di*dj)/(di + dj) ; imagb = di**2/dj**2
                 b_entry = dcmplx(realb, imagb)

                 ic = i -1 ; jc = j-1  ! convert to C indexing
                 call MatSetValue(A, ic, jc, a_entry, ADD_VALUES,ierr)
                 call MatSetValue(B, ic, jc, b_entry, ADD_VALUES,ierr)
            endif
        enddo
      enddo

      call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
      call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
      call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
      call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)

!     Check matrices
      write(*,*)'A matrix ... '
      call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)
      write(*,*)'B matrix ... '
      call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)

      call MatCreateVecs(A,PETSC_NULL_VEC,xr)
      call MatCreateVecs(A, PETSC_NULL_VEC,xi)

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Create the eigensolver and display info
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     ** Create eigensolver context
      call EPSCreate(PETSC_COMM_WORLD,eps,ierr)

!     ** Set operators.for general problem  Ax = lambda B x
      call EPSSetOperators(eps,A, B, ierr)
      call EPSSetProblemType(eps,EPS_GNHEP,ierr)

!     ** Set solver parameters at runtime
      call EPSSetFromOptions(eps,ierr)

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Solve the eigensystem
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      write(*,*)'rank',rank, 'entering solver ...'
      call EPSSolve(eps,ierr)

!     ** Free work space
      call EPSDestroy(eps,ierr)
      call MatDestroy(A,ierr)
      call MatDestroy(B,ierr)

      call VecDestroy(xr,ierr)
      call VecDestroy(xi,ierr)

      call SlepcFinalize(ierr)
      end


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210328/89a898c9/attachment.html>