<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">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_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>