[petsc-users] Singlar values of the GMRES Hessenberg matrix

[Dave Lee]

Hi PETSc,

I'm trying to add a "hook step" to the SNES trust region solver (at the end
of the function: KSPGMRESBuildSoln())

I'm testing this using the (linear) example:
src/ksp/ksp/examples/tutorials/ex1.c
as
gdb  --args ./test -snes_mf -snes_type newtontr -ksp_rtol 1.0e-12
-snes_stol 1.0e-12 -ksp_converged_reason -snes_converged_reason
-ksp_monitor -snes_monitor
(Ignore the SNES stuff, this is for when I test nonlinear examples).

When I call the LAPACK SVD routine via PETSc as
PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_(...))
I get the following singular values:

  0 KSP Residual norm 7.071067811865e-01
  1 KSP Residual norm 3.162277660168e-01
  2 KSP Residual norm 1.889822365046e-01
  3 KSP Residual norm 1.290994448736e-01
  4 KSP Residual norm 9.534625892456e-02
  5 KSP Residual norm 8.082545620881e-16

1 0.5 -7.85046e-16 1.17757e-15
0.5 1 0.5 1.7271e-15
0 0.5 1 0.5
0 0 0.5 1
0 0 0 0.5

singular values: 2.36264 0.409816 1.97794e-15 6.67632e-16

Linear solve converged due to CONVERGED_RTOL iterations 5

Where the lines above the singular values are the Hessenberg matrix that
I'm doing the SVD on.

When I build the solution in terms of the leading two right singular
vectors (and subsequently the first two orthonormal basis vectors in VECS_VV
I get an error norm as:
Norm of error 3.16228, Iterations 5

My suspicion is that I'm creating the Hessenberg incorrectly, as I would
have thought that this problem should have more than two non-zero leading
singular values.

Within my modified version of the GMRES build solution function (attached)
I'm creating this (and passing it to LAPACK as):

    nRows = gmres->it+1;
    nCols = nRows-1;

    ierr = PetscBLASIntCast(nRows,&nRows_blas);CHKERRQ(ierr);
    ierr = PetscBLASIntCast(nCols,&nCols_blas);CHKERRQ(ierr);
    ierr = PetscBLASIntCast(5*nRows,&lwork);CHKERRQ(ierr);
    ierr = PetscMalloc1(5*nRows,&work);CHKERRQ(ierr);
    ierr = PetscMalloc1(nRows*nCols,&R);CHKERRQ(ierr);
    ierr = PetscMalloc1(nRows*nCols,&H);CHKERRQ(ierr);
    for (jj = 0; jj < nRows; jj++) {
      for (ii = 0; ii < nCols; ii++) {
        R[jj*nCols+ii] = *HES(jj,ii);
      }
    }
    // Duplicate the Hessenberg matrix as the one passed to the SVD solver
is destroyed
    for (ii=0; ii<nRows*nCols; ii++) H[ii] = R[ii];

    ierr = PetscMalloc1(nRows*nRows,&U);CHKERRQ(ierr);
    ierr = PetscMalloc1(nCols*nCols,&VT);CHKERRQ(ierr);
    ierr = PetscMalloc1(nRows*nRows,&UT);CHKERRQ(ierr);
    ierr = PetscMalloc1(nCols*nCols,&V);CHKERRQ(ierr);
    ierr = PetscMalloc1(nRows,&p);CHKERRQ(ierr);
    ierr = PetscMalloc1(nCols,&q);CHKERRQ(ierr);
    ierr = PetscMalloc1(nCols,&y);CHKERRQ(ierr);

    // Perform an SVD on the Hessenberg matrix - Note: this call destroys
the input Hessenberg
    ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nRows_blas,&nCols_blas,R,&nRows_blas,S,UT,&nRows_blas,V,&nCols_blas,work,&lwork,&lierr));
    if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack
routine %d",(int)lierr);
    ierr = PetscFPTrapPop();CHKERRQ(ierr);

    // Find the number of non-zero singular values
    for(nnz=0; nnz<nCols; nnz++) {
      if(fabs(S[nnz]) < 1.0e-8) break;
    }
    printf("number of nonzero singular values: %d\n",nnz);

    trans(nRows,nRows,UT,U);
    trans(nCols,nCols,V,VT);

    // Compute p = ||r_0|| U^T e_1
    beta = gmres->res_beta;
    for (ii=0; ii<nCols; ii++) {
      p[ii] = beta*UT[ii*nRows];
    }
    p[nCols] = 0.0;

    // Original GMRES solution (\mu = 0)
    for (ii=0; ii<nnz; ii++) {
      q[ii] = p[ii]/S[ii];
    }

    // Expand y in terms of the right singular vectors as y = V q
    for (jj=0; jj<nnz; jj++) {
      y[jj] = 0.0;
      for (ii=0; ii<nCols; ii++) {
        y[jj] += V[jj*nCols+ii]*q[ii]; // transpose of the transpose
      }
    }

    // Pass the orthnomalized Krylov vector weights back out
    for (ii=0; ii<nnz; ii++) {
      nrs[ii] = y[ii];
    }

I just wanted to check that this is the correct way to extract the
Hessenberg from the KSP_GMRES structure, and to pass it to LAPACK, and if
so, should I really be expecting only two non-zero singular values in
return for this problem?

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