<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 Jose,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Many thanks for your reply, it enabled me to solve the problem I described below.</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
In particular, I recompiled PETSC and SLEPC with MUMPS and SCALAPACK included <br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
(by including the following lines in the configuration file:</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
'--download-mumps', <br>
<div> '--download-scalapack',</div>
'--download-cmake',</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
)</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
and then the program (unchanged from the previous attempt) then both compiled</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
and ran straightaway with no hanging. Just in case someone is following this example later,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
the correct eigenvalues are (I have checked them against an independent lapack zggev code):</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
0 (15.7283987073479, 79.3812583009335)
<div> 1 (10.3657189951037, 65.4935496512632)</div>
<div> 2 (20.2726807729152, 60.7235264113338)</div>
<div> 3 (15.8693370278539, 54.4403170495817)</div>
<div> 4 (8.93270390707530, 42.0676105026967)</div>
<div> 5 (18.0161334989426, 31.7976217614629)</div>
<div> 6 (16.2219350827311, 26.7999463239364)</div>
<div> 7 (6.64653598248233, 19.2535093354505)</div>
<div> 8 (7.23494239184217, 4.58606776802574)</div>
<div> 9 (3.68158090200136, 1.65838104812904)</div>
<br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
The problem would then indeed appear to have been - exactly as you suggested - <br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
the fact that the generalised eigenvalue problem requires linear systems to be</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
solved as part of the solution process, but in the previous attempt I was trying</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
to solve an MPI-distributed system using a sequential solver.<br>
</div>
<div> </div>
<div>A couple of other quick points / queries:</div>
<div>1) By 'debug mode' you mean compiling the library with '--with-debugging=1'?</div>
<div>2) Out of interest, I am wondering out of interest which LU solver was used - scalapack or
<br>
</div>
<div> MUMPS ? I could see both these libraries in the linking command, is there an easy way
<br>
</div>
<div> to find out which solver was actually called ?<br>
</div>
<div>2) One slightly strange thing is that I also compiled the library with '--with-64-bit-indices=1'</div>
<div> for 64-bit memory addressing, but I noticed in the compilation commands generated
<br>
</div>
<div>by the make command there is no '-i8' flag, which is used with mpiifort to request 64 bit</div>
<div>integers. Is it the case that this is not required since everything is somehow taken care of by</div>
<div>Petsc custom data type ( PetscInt) ? As an experiment I tried manually inserting the -i8</div>
<div>and got a few errors like:</div>
<div><span style="color: rgb(12, 100, 192);">/data/work/rotplane/omega_to_zero/stability/test/tmp10/tmp3/write_and_solve8.F(38): warning #6075: The data type of the actual argument does not match the definition. [PETSC_COMM_WORLD]</span>
<div><span style="color: rgb(12, 100, 192);"> call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);if (ierr .ne. 0) th</span></div>
<br>
</div>
<div> <br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
A big thank you once again, <br>
</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Best wishes,</div>
<div style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
Dan.<br>
</div>
<div>
<div id="appendonsend"></div>
<div style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt; color:rgb(0,0,0)">
<br>
</div>
<hr tabindex="-1" style="display:inline-block; width:98%">
<div id="divRplyFwdMsg" dir="ltr"><font style="font-size:11pt" face="Calibri, sans-serif" color="#000000"><b>From:</b> Jose E. Roman <jroman@dsic.upv.es><br>
<b>Sent:</b> Sunday, March 28, 2021 4:42 PM<br>
<b>To:</b> dazza simplythebest <sayosale@hotmail.com><br>
<b>Cc:</b> PETSc users list <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Newbie question: Something is wrong with this Slepc simple example for generalised eigenvalue problem</font>
<div> </div>
</div>
<div class="BodyFragment"><font size="2"><span style="font-size:11pt">
<div class="PlainText">You should run in debug mode until you get a correct code, otherwise you may no see some error messages. Also, it is recommended to add error checking after every call to PETSc, otherwise the execution continues and may get blocked. See
the CHKERRA macro in <a href="https://slepc.upv.es/documentation/current/src/eps/tutorials/ex1f90.F90.html">
https://slepc.upv.es/documentation/current/src/eps/tutorials/ex1f90.F90.html</a><br>
<br>
The problem you are probably having is that you are running with several MPI processes, so you need a parallel LU solver. See FAQ #10
<a href="https://slepc.upv.es/documentation/faq.htm">https://slepc.upv.es/documentation/faq.htm</a><br>
<br>
Jose<br>
<br>
<br>
> El 28 mar 2021, a las 11:25, dazza simplythebest <sayosale@hotmail.com> escribió:<br>
> <br>
> Dear All,<br>
> I am seeking to use slepc/petsc to solve a generalised eigenvalue problem<br>
> that arises in a hydrodynamic stability problem. The code is a parallelisation of an existing
<br>
> serial Fortran code. Before I get to grips with this target problem, I would of course like to get<br>
> some relevant examples working. I have installed petsc/ slepc seemingly without any problem,<br>
> and the provided slepc example fortran program ex1f.F, which solves a regular eigenvalue<br>
> problem Ax = lambda x, seemed to compile and run correctly.<br>
> I have now written a short program to instead solve the complex generalised <br>
> problem Ax = lambda B x (see below) . This code compiles and runs w/out errors but<br>
> for some reason hangs when calling EPSSolve - we enter EPSSolve but never leave.<br>
> The matrices appear to be correctly assembled -all the values are correct in the Matview<br>
> printout, so I am not quite sure where I have gone wrong, can anyone spot my mistake?<br>
> ( Note that for the actual problem I wish to solve I have already written the code to construct the matrix,<br>
> which distributes the rows across the processes and it is fully tested and working. Hence I want to specify<br>
> the distribution of rows and not leave it up to a PETS_DECIDE .)<br>
> I would be very grateful if someone can point out what is wrong with this small example code (see below),<br>
> Many thanks,<br>
> Dan.<br>
> <br>
> ! this program MUST be run with NGLOBAL = 10, MY_NO_ROWS = 5<br>
> ! and two MPI processes since row distribution is hard-baked into code <br>
> !<br>
> program main<br>
> #include <slepc/finclude/slepceps.h><br>
> use slepceps<br>
> implicit none<br>
> <br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! Declarations<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> !<br>
> ! Variables:<br>
> ! A , B double complex operator matrices<br>
> ! eps eigenproblem solver context<br>
> <br>
> Mat A,B<br>
> EPS eps<br>
> EPSType tname<br>
> PetscReal tol, error<br>
> PetscScalar kr, ki<br>
> Vec xr, xi<br>
> PetscInt NGLOBAL , MY_NO_ROWS, NL3, owner<br>
> PetscInt nev, maxit, its, nconv<br>
> PetscInt i,j,ic,jc<br>
> PetscReal reala, imaga, realb, imagb, di, dj<br>
> PetscScalar a_entry, b_entry<br>
> PetscMPIInt rank<br>
> PetscErrorCode ierr<br>
> PetscInt,parameter :: zero = 0, one = 1, two = 2, three = 3 <br>
> <br>
> PetscInt M1, N1, mm1, nn1<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! Beginning of program<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> <br>
> call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)<br>
> call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)<br>
> ! make sure you set NGLOBAL = 10, MY_NO_ROWS = 5 and run with two processes<br>
> NGLOBAL = 10<br>
> MY_NO_ROWS = 5<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! Compute the operator matrices that define the eigensystem, Ax=kBx<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> !!!!!!!!! Setup A matrix<br>
> <br>
> call MatCreate(PETSC_COMM_WORLD,A,ierr) <br>
> call MatSetsizes(A,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) <br>
> call MatSetFromOptions(A,ierr)<br>
> call MatGetSize(A,M1,N1,ierr)<br>
> write(*,*)'Rank [',rank,']: global size of A is ',M1, N1<br>
> call MatGetLocalSize(A,mm1,nn1,ierr)<br>
> write(*,*)'Rank [',rank,']: my local size of A is ',mm1, nn1<br>
> call MatMPIAIJSetPreallocation(A,three, PETSC_NULL_INTEGER,one, &<br>
> & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation<br>
> <br>
> !!!!!!!!! Setup B matrix<br>
> call MatCreate(PETSC_COMM_WORLD,B,ierr) <br>
> call MatSetsizes(B,MY_NO_ROWS, MY_NO_ROWS ,NGLOBAL,NGLOBAL,ierr) <br>
> call MatSetFromOptions(B,ierr)<br>
> call MatGetSize(B,M1,N1,ierr)<br>
> write(*,*)'Rank [',rank,']: global size of B is ',M1, N1<br>
> call MatGetLocalSize(B,mm1,nn1,ierr)<br>
> write(*,*)'Rank [',rank,']: my local size of B is ',mm1, nn1<br>
> <br>
> call MatMPIAIJSetPreallocation(B,three, PETSC_NULL_INTEGER,one, &<br>
> & PETSC_NULL_INTEGER,ierr) !parallel (MPI) allocation<br>
> <br>
> ! initalise<br>
> call MatZeroEntries(A,ierr)<br>
> call MatZeroEntries(B,ierr)<br>
> <br>
> ! Fill in values of A, B and assemble matrices<br>
> ! Both matrices are tridiagonal with<br>
> ! Aij = cmplx( (i-j)**2, (i+j)**2)<br>
> ! Bij = cmplx( ij/i + j, (i/j)**2)<br>
> ! (numbering from 1 )<br>
> <br>
> do i = 1, NGLOBAL<br>
> ! a rather crude way to distribute rows <br>
> if (i < 6) owner = 0<br>
> if (i >= 6) owner = 1<br>
> if (rank /= owner) cycle<br>
> do j = 1, NGLOBAL<br>
> if ( abs(i-j) < 2 ) then<br>
> write(*,*)rank,' : Setting ',i,j<br>
> di = dble(i) ; dj = dble(j)<br>
> <br>
> reala = (di - dj)**2 ; imaga = (di + dj)**2<br>
> a_entry = dcmplx(reala, imaga)<br>
> realb = (di*dj)/(di + dj) ; imagb = di**2/dj**2<br>
> b_entry = dcmplx(realb, imagb) <br>
> <br>
> ic = i -1 ; jc = j-1 ! convert to C indexing<br>
> call MatSetValue(A, ic, jc, a_entry, ADD_VALUES,ierr)<br>
> call MatSetValue(B, ic, jc, b_entry, ADD_VALUES,ierr)<br>
> endif<br>
> enddo<br>
> enddo<br>
> <br>
> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) <br>
> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)<br>
> call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) <br>
> call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)<br>
> <br>
> ! Check matrices <br>
> write(*,*)'A matrix ... '<br>
> call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
> write(*,*)'B matrix ... '<br>
> call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)<br>
> <br>
> call MatCreateVecs(A,PETSC_NULL_VEC,xr)<br>
> call MatCreateVecs(A, PETSC_NULL_VEC,xi)<br>
> <br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! Create the eigensolver and display info<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! ** Create eigensolver context<br>
> call EPSCreate(PETSC_COMM_WORLD,eps,ierr)<br>
> <br>
> ! ** Set operators.for general problem Ax = lambda B x<br>
> call EPSSetOperators(eps,A, B, ierr) <br>
> call EPSSetProblemType(eps,EPS_GNHEP,ierr)<br>
> <br>
> ! ** Set solver parameters at runtime<br>
> call EPSSetFromOptions(eps,ierr)<br>
> <br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> ! Solve the eigensystem<br>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> <br>
> write(*,*)'rank',rank, 'entering solver ...'<br>
> call EPSSolve(eps,ierr)<br>
> <br>
> ! ** Free work space<br>
> call EPSDestroy(eps,ierr)<br>
> call MatDestroy(A,ierr)<br>
> call MatDestroy(B,ierr)<br>
> <br>
> call VecDestroy(xr,ierr)<br>
> call VecDestroy(xi,ierr)<br>
> <br>
> call SlepcFinalize(ierr)<br>
> end<br>
<br>
</div>
</span></font></div>
</div>
</body>
</html>