<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.<br></div><div><br></div><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>