<div dir="ltr"><div dir="ltr">On Fri, May 24, 2019 at 8:38 AM Dave Lee <<a href="mailto:davelee2804@gmail.com">davelee2804@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Thanks Matt, great suggestion.<div><br></div><div>I did indeed find a transpose error this way. The SVD as reconstructed via U S V^T now matches the input Hessenberg matrix as derived via the <font face="courier new, monospace">*HES(row,col)</font> macro, and all the singular values are non-zero. However the solution to example <span style="font-family:"courier new",monospace">src/ksp/ksp/examples/</span><span style="font-family:"courier new",monospace">tutorials/ex1.c </span><font face="arial, sans-serif">as determined via the expansion over the singular vectors is still not correct. I suspect I'm doing something wrong with regards to the expansion over the vec array </font><font face="courier new, monospace">VEC_VV()</font><font face="arial, sans-serif">, which I assume are the orthonormal vectors of the </font><font face="courier new, monospace">Q_k</font><font face="arial, sans-serif"> matrix in the Arnoldi iteration....</font></div></div></blockquote><div><br></div><div>Here we are building the solution:</div><div><br></div><div>  <a href="https://bitbucket.org/petsc/petsc/src/7c23e6aa64ffbff85a2457e1aa154ec3d7f238e3/src/ksp/ksp/impls/gmres/gmres.c#lines-331">https://bitbucket.org/petsc/petsc/src/7c23e6aa64ffbff85a2457e1aa154ec3d7f238e3/src/ksp/ksp/impls/gmres/gmres.c#lines-331</a></div><div><br></div><div>There are some subtleties if you have a  nonzero initial guess or a preconditioner.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><font face="arial, sans-serif">Thanks again for your advice, I'll keep digging.</font></div><div><font face="arial, sans-serif"><br></font></div><div><font face="arial, sans-serif">Cheers, Dave.</font></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, May 23, 2019 at 8:20 PM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Thu, May 23, 2019 at 5:09 AM Dave Lee via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hi PETSc,<div><br></div><div>I'm trying to add a "hook step" to the SNES trust region solver (at the end of the function: <font face="courier new, monospace">KSPGMRESBuildSoln()</font>)</div><div><br></div><div>I'm testing this using the (linear) example: </div><div><font face="courier new, monospace">src/ksp/ksp/examples/tutorials/ex1.c</font><br></div><div>as </div><div><font face="courier new, monospace">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</font><br></div><div>(Ignore the SNES stuff, this is for when I test nonlinear examples).</div><div><br></div><div>When I call the LAPACK SVD routine via PETSc as</div><div><font face="courier new, monospace">PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_(...))</font><br></div><div>I get the following singular values:</div><div><br></div><div><font face="courier new, monospace">  0 KSP Residual norm 7.071067811865e-01 <br>  1 KSP Residual norm 3.162277660168e-01 <br>  2 KSP Residual norm 1.889822365046e-01 <br>  3 KSP Residual norm 1.290994448736e-01 <br>  4 KSP Residual norm 9.534625892456e-02 <br>  5 KSP Residual norm 8.082545620881e-16 <br><br>  1       0.5     -7.85046e-16    1.17757e-15<br>   0.5     1       0.5     1.7271e-15<br>    0       0.5     1       0.5<br>   0       0       0.5     1<br>     0       0       0       0.5<br><br>singular values:         2.36264 0.409816        1.97794e-15     6.67632e-16<br><br>Linear solve converged due to CONVERGED_RTOL iterations 5</font><br></div><div><br></div><div>Where the lines above the singular values are the Hessenberg matrix that I'm doing the SVD on.</div></div></blockquote><div><br></div><div>First, write out all the SVD matrices you get and make sure that they reconstruct the input matrix (that</div><div>you do not have something transposed somewhere).</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>When I build the solution in terms of the leading two right singular vectors (and subsequently the first two orthonormal basis vectors in <font face="courier new, monospace">VECS_VV</font> I get an error norm as:</div><div><font face="courier new, monospace">Norm of error 3.16228, Iterations 5</font><br></div><div><br></div><div>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.</div><div><br></div><div>Within my modified version of the GMRES build solution function (attached) I'm creating this (and passing it to LAPACK as):</div><div><br></div><div><font face="courier new, monospace">    nRows = gmres->it+1;<br>    nCols = nRows-1;<br><br>    ierr = PetscBLASIntCast(nRows,&nRows_blas);CHKERRQ(ierr);<br>    ierr = PetscBLASIntCast(nCols,&nCols_blas);CHKERRQ(ierr);<br>    ierr = PetscBLASIntCast(5*nRows,&lwork);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(5*nRows,&work);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nRows*nCols,&R);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nRows*nCols,&H);CHKERRQ(ierr);<br>    for (jj = 0; jj < nRows; jj++) {<br>      for (ii = 0; ii < nCols; ii++) {<br>        R[jj*nCols+ii] = *HES(jj,ii);<br>      }<br>    }<br>    // Duplicate the Hessenberg matrix as the one passed to the SVD solver is destroyed<br>    for (ii=0; ii<nRows*nCols; ii++) H[ii] = R[ii];<br></font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">    ierr = PetscMalloc1(nRows*nRows,&U);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nCols*nCols,&VT);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nRows*nRows,&UT);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nCols*nCols,&V);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nRows,&p);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nCols,&q);CHKERRQ(ierr);<br>    ierr = PetscMalloc1(nCols,&y);CHKERRQ(ierr);<br><br>    // Perform an SVD on the Hessenberg matrix - Note: this call destroys the input Hessenberg<br>    ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nRows_blas,&nCols_blas,R,&nRows_blas,S,UT,&nRows_blas,V,&nCols_blas,work,&lwork,&lierr));<br>    if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr);<br>    ierr = PetscFPTrapPop();CHKERRQ(ierr);</font><br></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">    // Find the number of non-zero singular values<br>    for(nnz=0; nnz<nCols; nnz++) {<br>      if(fabs(S[nnz]) < 1.0e-8) break;<br>    }<br>    printf("number of nonzero singular values: %d\n",nnz);<br><br>    trans(nRows,nRows,UT,U);<br>    trans(nCols,nCols,V,VT);<br><br>    // Compute p = ||r_0|| U^T e_1<br>    beta = gmres->res_beta;<br>    for (ii=0; ii<nCols; ii++) {<br>      p[ii] = beta*UT[ii*nRows];<br>    }<br>    p[nCols] = 0.0;<br></font></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">    // Original GMRES solution (\mu = 0)<br>    for (ii=0; ii<nnz; ii++) {<br>      q[ii] = p[ii]/S[ii];<br>    }<br></font></div><div><br></div><div><font face="courier new, monospace">    // Expand y in terms of the right singular vectors as y = V q<br>    for (jj=0; jj<nnz; jj++) {<br>      y[jj] = 0.0;<br>      for (ii=0; ii<nCols; ii++) {<br>        y[jj] += V[jj*nCols+ii]*q[ii]; // transpose of the transpose<br>      }<br>    }</font><br></div><div><font face="courier new, monospace"><br></font></div><div><font face="courier new, monospace">    // Pass the orthnomalized Krylov vector weights back out<br>    for (ii=0; ii<nnz; ii++) {<br>      nrs[ii] = y[ii];<br>    }<br></font></div><div><br></div><div>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?</div><div><br></div><div>Cheers, Dave.</div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail-m_8895125370612354309gmail-m_-4001304069175436629gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>