Re: [petsc-users] Calling LAPACK routines from PETSc
Thanks Barry,
I found some helpful examples on the intel lapack site - moral of the
story: using C ordering for input matrix, but transposed output matrices
leads to a consistent solution.
Cheers, Dave.
On Mon, May 20, 2019 at 6:07 PM Smith, Barry F. wrote:
>
>
> > On May 20, 2019, at 2:28 AM, Dave Lee wrote:
> >
> > Thanks Jed and Barry,
> >
> > So, just to confirm,
> >
> > -- From the KSP_GMRES structure, if I call *HH(a,b), that will return
> the row a, column b entry of the Hessenberg matrix (while the back end
> array *hh_origin array is ordering using the Fortran convention)
> >
> > -- Matrices are passed into and returned from
> PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_() using Fortran indexing, and
> need to be transposed to get back to C ordering
>
> In general, I guess depending on what you want to do with them you don'
> need to transpose them. Why would you want to? Just leave them as little
> column oriented blogs and with them what you need directly.
>
>Just do stuff and you'll find it works out.
>
> >
> > Are both of these statements correct?
> >
> > Cheers, Dave.
> >
> > On Mon, May 20, 2019 at 4:34 PM Smith, Barry F.
> wrote:
> >
> >The little work arrays in GMRES tend to be stored in Fortran
> ordering; there is no C style p[][] indexing into such arrays. Thus the
> arrays can safely be sent to LAPACK. The only trick is knowing the two
> dimensions and as Jed say the "leading dimension parameter. He gave you a
> place to look
> >
> > > On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
> > >
> > > Dave Lee via petsc-users writes:
> > >
> > >> 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?
> > >
> > > LAPACK terminology is with respect to Fortran ordering. There is a
> > > "leading dimension" parameter so that you can operate on non-contiguous
> > > blocks. See KSPComputeExtremeSingularValues_GMRES for an example.
> > >
> > >> 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.
> > >>
> > >> /*
> > >>This file implements GMRES (a Generalized Minimal Residual) method.
> > >>Reference: Saad and Schultz, 1986.
> > >>
> > >>
> > >>Some comments on left vs. right preconditioning, and restarts.
> > >>Left and right preconditioning.
> > >>If right preconditioning is chosen, then the problem being solved
> > >>by gmres is actually
> > >> My = AB^-1 y = f
> > >>so the initial residual is
> > >> r = f - Mx
> > >>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> > >>residual is
> > >> r = f - A x
> > >>The final solution is then
> > >> x = B^-1 y
> > >>
> > >>If left preconditioning is chosen, then the problem being solved is
> > >> My = B^-1 A x = B^-1 f,
> > >>and the initial residual is
> > >> r = B^-1(f - Ax)
> > >>
> > >>Restarts: Restarts are basically solves with x0 not equal to zero.
> > >>Note that we can eliminate an extra application of B^-1 between
> > >>restarts as long as we don't require that the solution at the end
> > >>of an unsuccessful gmres iteration always be the solution x.
> > >> */
> > >>
> > >> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I
> "petscksp.h" I*/
> > >> #include // DRL
> > >> #define GMRES_DELTA_DIRECTIONS 10
> > >> #define GMRES_DEFAULT_MAXK 30
> > >> static PetscErrorCode
> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> > >> static PetscErrorCode
> KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
> > >>
> > >> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> > >> {
> > >> PetscInt
Re: [petsc-users] Calling LAPACK routines from PETSc
> On May 20, 2019, at 2:28 AM, Dave Lee wrote:
>
> Thanks Jed and Barry,
>
> So, just to confirm,
>
> -- From the KSP_GMRES structure, if I call *HH(a,b), that will return the row
> a, column b entry of the Hessenberg matrix (while the back end array
> *hh_origin array is ordering using the Fortran convention)
>
> -- Matrices are passed into and returned from
> PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_() using Fortran indexing, and
> need to be transposed to get back to C ordering
In general, I guess depending on what you want to do with them you don' need
to transpose them. Why would you want to? Just leave them as little column
oriented blogs and with them what you need directly.
Just do stuff and you'll find it works out.
>
> Are both of these statements correct?
>
> Cheers, Dave.
>
> On Mon, May 20, 2019 at 4:34 PM Smith, Barry F. wrote:
>
>The little work arrays in GMRES tend to be stored in Fortran ordering;
> there is no C style p[][] indexing into such arrays. Thus the arrays can
> safely be sent to LAPACK. The only trick is knowing the two dimensions and
> as Jed say the "leading dimension parameter. He gave you a place to look
>
> > On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users
> > wrote:
> >
> > Dave Lee via petsc-users writes:
> >
> >> 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?
> >
> > LAPACK terminology is with respect to Fortran ordering. There is a
> > "leading dimension" parameter so that you can operate on non-contiguous
> > blocks. See KSPComputeExtremeSingularValues_GMRES for an example.
> >
> >> 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.
> >>
> >> /*
> >>This file implements GMRES (a Generalized Minimal Residual) method.
> >>Reference: Saad and Schultz, 1986.
> >>
> >>
> >>Some comments on left vs. right preconditioning, and restarts.
> >>Left and right preconditioning.
> >>If right preconditioning is chosen, then the problem being solved
> >>by gmres is actually
> >> My = AB^-1 y = f
> >>so the initial residual is
> >> r = f - Mx
> >>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> >>residual is
> >> r = f - A x
> >>The final solution is then
> >> x = B^-1 y
> >>
> >>If left preconditioning is chosen, then the problem being solved is
> >> My = B^-1 A x = B^-1 f,
> >>and the initial residual is
> >> r = B^-1(f - Ax)
> >>
> >>Restarts: Restarts are basically solves with x0 not equal to zero.
> >>Note that we can eliminate an extra application of B^-1 between
> >>restarts as long as we don't require that the solution at the end
> >>of an unsuccessful gmres iteration always be the solution x.
> >> */
> >>
> >> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I "petscksp.h"
> >> I*/
> >> #include // DRL
> >> #define GMRES_DELTA_DIRECTIONS 10
> >> #define GMRES_DEFAULT_MAXK 30
> >> static PetscErrorCode
> >> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> >> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
> >>
> >> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> >> {
> >> PetscInt hh,hes,rs,cc;
> >> PetscErrorCode ierr;
> >> PetscInt max_k,k;
> >> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
> >>
> >> PetscFunctionBegin;
> >> max_k = gmres->max_k; /* restart size */
> >> hh= (max_k + 2) * (max_k + 1);
> >> hes = (max_k + 1) * (max_k + 1);
> >> rs= (max_k + 2);
> >> cc= (max_k + 1);
> >>
> >> ierr =
> >>
Re: [petsc-users] Calling LAPACK routines from PETSc
The little work arrays in GMRES tend to be stored in Fortran ordering; there
is no C style p[][] indexing into such arrays. Thus the arrays can safely be
sent to LAPACK. The only trick is knowing the two dimensions and as Jed say
the "leading dimension parameter. He gave you a place to look
> On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users
> wrote:
>
> Dave Lee via petsc-users writes:
>
>> 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?
>
> LAPACK terminology is with respect to Fortran ordering. There is a
> "leading dimension" parameter so that you can operate on non-contiguous
> blocks. See KSPComputeExtremeSingularValues_GMRES for an example.
>
>> 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.
>>
>> /*
>>This file implements GMRES (a Generalized Minimal Residual) method.
>>Reference: Saad and Schultz, 1986.
>>
>>
>>Some comments on left vs. right preconditioning, and restarts.
>>Left and right preconditioning.
>>If right preconditioning is chosen, then the problem being solved
>>by gmres is actually
>> My = AB^-1 y = f
>>so the initial residual is
>> r = f - Mx
>>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
>>residual is
>> r = f - A x
>>The final solution is then
>> x = B^-1 y
>>
>>If left preconditioning is chosen, then the problem being solved is
>> My = B^-1 A x = B^-1 f,
>>and the initial residual is
>> r = B^-1(f - Ax)
>>
>>Restarts: Restarts are basically solves with x0 not equal to zero.
>>Note that we can eliminate an extra application of B^-1 between
>>restarts as long as we don't require that the solution at the end
>>of an unsuccessful gmres iteration always be the solution x.
>> */
>>
>> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I "petscksp.h"
>> I*/
>> #include // DRL
>> #define GMRES_DELTA_DIRECTIONS 10
>> #define GMRES_DEFAULT_MAXK 30
>> static PetscErrorCode
>> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
>> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
>>
>> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
>> {
>> PetscInt hh,hes,rs,cc;
>> PetscErrorCode ierr;
>> PetscInt max_k,k;
>> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
>>
>> PetscFunctionBegin;
>> max_k = gmres->max_k; /* restart size */
>> hh= (max_k + 2) * (max_k + 1);
>> hes = (max_k + 1) * (max_k + 1);
>> rs= (max_k + 2);
>> cc= (max_k + 1);
>>
>> ierr =
>> PetscCalloc5(hh,>hh_origin,hes,>hes_origin,rs,>rs_origin,cc,>cc_origin,cc,>ss_origin);CHKERRQ(ierr);
>> ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs +
>> 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);
>>
>> if (ksp->calc_sings) {
>>/* Allocate workspace to hold Hessenberg matrix needed by lapack */
>>ierr = PetscMalloc1((max_k + 3)*(max_k + 9),>Rsvd);CHKERRQ(ierr);
>>ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k +
>> 9)*sizeof(PetscScalar));CHKERRQ(ierr);
>>ierr = PetscMalloc1(6*(max_k+2),>Dsvd);CHKERRQ(ierr);
>>ierr =
>> PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);
>> }
>>
>> /* Allocate array to hold pointers to user vectors. Note that we need
>> 4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
>> gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
>>
>> ierr = PetscMalloc1(gmres->vecs_allocated,>vecs);CHKERRQ(ierr);
>> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>user_work);CHKERRQ(ierr);
>> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>mwork_alloc);CHKERRQ(ierr);
>> ierr =
>>
Re: [petsc-users] Calling LAPACK routines from PETSc
Dave Lee via petsc-users writes:
> 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?
LAPACK terminology is with respect to Fortran ordering. There is a
"leading dimension" parameter so that you can operate on non-contiguous
blocks. See KSPComputeExtremeSingularValues_GMRES for an example.
> 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.
>
> /*
> This file implements GMRES (a Generalized Minimal Residual) method.
> Reference: Saad and Schultz, 1986.
>
>
> Some comments on left vs. right preconditioning, and restarts.
> Left and right preconditioning.
> If right preconditioning is chosen, then the problem being solved
> by gmres is actually
>My = AB^-1 y = f
> so the initial residual is
> r = f - Mx
> Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> residual is
> r = f - A x
> The final solution is then
> x = B^-1 y
>
> If left preconditioning is chosen, then the problem being solved is
>My = B^-1 A x = B^-1 f,
> and the initial residual is
>r = B^-1(f - Ax)
>
> Restarts: Restarts are basically solves with x0 not equal to zero.
> Note that we can eliminate an extra application of B^-1 between
> restarts as long as we don't require that the solution at the end
> of an unsuccessful gmres iteration always be the solution x.
> */
>
> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I "petscksp.h" I*/
> #include // DRL
> #define GMRES_DELTA_DIRECTIONS 10
> #define GMRES_DEFAULT_MAXK 30
> static PetscErrorCode
> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
>
> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> {
> PetscInt hh,hes,rs,cc;
> PetscErrorCode ierr;
> PetscInt max_k,k;
> KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
>
> PetscFunctionBegin;
> max_k = gmres->max_k; /* restart size */
> hh= (max_k + 2) * (max_k + 1);
> hes = (max_k + 1) * (max_k + 1);
> rs= (max_k + 2);
> cc= (max_k + 1);
>
> ierr =
> PetscCalloc5(hh,>hh_origin,hes,>hes_origin,rs,>rs_origin,cc,>cc_origin,cc,>ss_origin);CHKERRQ(ierr);
> ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs +
> 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);
>
> if (ksp->calc_sings) {
> /* Allocate workspace to hold Hessenberg matrix needed by lapack */
> ierr = PetscMalloc1((max_k + 3)*(max_k + 9),>Rsvd);CHKERRQ(ierr);
> ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k +
> 9)*sizeof(PetscScalar));CHKERRQ(ierr);
> ierr = PetscMalloc1(6*(max_k+2),>Dsvd);CHKERRQ(ierr);
> ierr =
> PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);
> }
>
> /* Allocate array to hold pointers to user vectors. Note that we need
>4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
> gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
>
> ierr = PetscMalloc1(gmres->vecs_allocated,>vecs);CHKERRQ(ierr);
> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>user_work);CHKERRQ(ierr);
> ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>mwork_alloc);CHKERRQ(ierr);
> ierr =
> PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt))
> + gmres->vecs_allocated*sizeof(Vec));CHKERRQ(ierr);
>
> if (gmres->q_preallocate) {
> gmres->vv_allocated = VEC_OFFSET + 2 + max_k;
>
> ierr =
> KSPCreateVecs(ksp,gmres->vv_allocated,>user_work[0],0,NULL);CHKERRQ(ierr);
> ierr =
> PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);CHKERRQ(ierr);
>
> gmres->mwork_alloc[0] = gmres->vv_allocated;
> gmres->nwork_alloc= 1;
>
