[petsc-users] Calling LAPACK routines from PETSc

[Dave Lee]

Hi Petsc,

I'm attempting to implement a "hookstep" for the SNES trust region solver.
Essentially what I'm trying to do is replace the solution of the least
squares problem at the end of each GMRES solve with a modified solution
with a norm that is constrained to be within the size of the trust region.

In order to do this I need to perform an SVD on the Hessenberg matrix,
which copying the function KSPComputeExtremeSingularValues(), I'm trying to
do by accessing the LAPACK function dgesvd() via the PetscStackCallBLAS()
machinery. One thing I'm confused about however is the ordering of the 2D
arrays into and out of this function, given that that C and FORTRAN arrays
use reverse indexing, ie: C[j+1][i+1] = F[i,j].

Given that the Hessenberg matrix has k+1 rows and k columns, should I be
still be initializing this as H[row][col] and passing this into
PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))
or should I be transposing this before passing it in?

Also for the left and right singular vector matrices that are returned by
this function, should I be transposing these before I interpret them as C
arrays?

I've attached my modified version of gmres.c in case this is helpful. If
you grep for DRL (my initials) then you'll see my changes to the code.

Cheers, Dave.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190520/e7e7e233/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gmres.c
Type: application/octet-stream
Size: 38311 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190520/e7e7e233/attachment-0001.obj>