<div dir="ltr">Hi Matt and PETSc.<div><br></div><div>Thanks again for the advice.</div><div><br></div><div>So I think I know what my problem might be. Looking at the comments above the function</div><div><font face="courier new, monospace">KSPInitialResidual()</font><br></div><div><font face="arial, sans-serif">in</font></div><div><font face="courier new, monospace">src/ksp/ksp/interface/itres.c<br></font></div><div>I see that the initial residual, as passed into <font face="courier new, monospace">VEC_VV(0)</font> is the residual of the <i>preconditioned</i> system (and that the original residual goes temporarily to <font face="courier new, monospace">gmres->vecs[1]</font>). </div><div><br></div><div>So I'm wondering, is the Hessenberg, as derived via the <font face="courier new, monospace">*HES(row,col)</font> macro the Hessenberg for the original Krylov subspace, or the preconditioned subspace?</div><div><br></div><div>Secondly, do the vecs within the <font face="courier new, monospace">KSP_GMRES</font> structure, as accessed via <font face="courier new, monospace">VEC_VV()</font> correspond to the (preconditioned) Krylov subspace or the orthonormalized vectors that make up the matrix <font face="courier new, monospace">Q_k</font> in the Arnoldi iteration? This isn't clear to me, and I need to access the vectors in <font face="courier new, monospace">Q_k</font> in order to expand the corrected hookstep solution.</div><div><br></div><div>Thanks again, Dave.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, May 25, 2019 at 6:18 PM Dave Lee <<a href="mailto:davelee2804@gmail.com">davelee2804@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">Thanks Matt, this is where I'm adding in my hookstep code.<div><br></div><div>Cheers, Dave.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, May 24, 2019 at 10:49 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 Fri, May 24, 2019 at 8:38 AM Dave Lee <<a href="mailto:davelee2804@gmail.com" target="_blank">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" target="_blank">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_-7729034548834137429gmail-m_-548281515559377497gmail-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-m_-7729034548834137429gmail-m_-548281515559377497gmail_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>