[petsc-users] Singlar values of the GMRES Hessenberg matrix
Dave Lee
davelee2804 at gmail.com
Mon May 27 02:55:33 CDT 2019
Hi Matt and PETSc.
Thanks again for the advice.
So I think I know what my problem might be. Looking at the comments above
the function
KSPInitialResidual()
in
src/ksp/ksp/interface/itres.c
I see that the initial residual, as passed into VEC_VV(0) is the residual
of the *preconditioned* system (and that the original residual goes
temporarily to gmres->vecs[1]).
So I'm wondering, is the Hessenberg, as derived via the *HES(row,col) macro
the Hessenberg for the original Krylov subspace, or the preconditioned
subspace?
Secondly, do the vecs within the KSP_GMRES structure, as accessed via
VEC_VV() correspond to the (preconditioned) Krylov subspace or the
orthonormalized vectors that make up the matrix Q_k in the Arnoldi
iteration? This isn't clear to me, and I need to access the vectors in Q_k in
order to expand the corrected hookstep solution.
Thanks again, Dave.
On Sat, May 25, 2019 at 6:18 PM Dave Lee <davelee2804 at gmail.com> wrote:
> Thanks Matt, this is where I'm adding in my hookstep code.
>
> Cheers, Dave.
>
> On Fri, May 24, 2019 at 10:49 PM Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Fri, May 24, 2019 at 8:38 AM Dave Lee <davelee2804 at gmail.com> wrote:
>>
>>> Thanks Matt, great suggestion.
>>>
>>> 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
>>> *HES(row,col) macro, and all the singular values are non-zero. However
>>> the solution to example src/ksp/ksp/examples/tutorials/ex1.c 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 VEC_VV(), which I assume are the orthonormal vectors
>>> of the Q_k matrix in the Arnoldi iteration....
>>>
>>
>> Here we are building the solution:
>>
>>
>> https://bitbucket.org/petsc/petsc/src/7c23e6aa64ffbff85a2457e1aa154ec3d7f238e3/src/ksp/ksp/impls/gmres/gmres.c#lines-331
>>
>> There are some subtleties if you have a nonzero initial guess or a
>> preconditioner.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> Thanks again for your advice, I'll keep digging.
>>>
>>> Cheers, Dave.
>>>
>>> On Thu, May 23, 2019 at 8:20 PM Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Thu, May 23, 2019 at 5:09 AM Dave Lee via petsc-users <
>>>> petsc-users at mcs.anl.gov> wrote:
>>>>
>>>>> 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.
>>>>>
>>>>
>>>> First, write out all the SVD matrices you get and make sure that they
>>>> reconstruct the input matrix (that
>>>> you do not have something transposed somewhere).
>>>>
>>>> Matt
>>>>
>>>>
>>>>> 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.
>>>>>
>>>>
>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which their
>>>> experiments lead.
>>>> -- Norbert Wiener
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/
>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>
>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190527/63ded96d/attachment.html>
More information about the petsc-users
mailing list