[petsc-users] Singlar values of the GMRES Hessenberg matrix
Matthew Knepley
knepley at gmail.com
Mon May 27 06:52:31 CDT 2019
On Mon, May 27, 2019 at 7:46 AM Dave Lee <davelee2804 at gmail.com> wrote:
> Thanks Matt, so just to clarify:
>
> -- VEC_VV() contains the Krylov subspace vectors (left preconditioned if
> PC_LEFT) and not the orthonomalized vectors that make up Q_k?
> - if so, is it possible to obtain Q_k?
>
> -- HES(row,col) contains the entries of the Hessenberg matrix
> corresponding to the Arnoldi iteration for the preconditioned Krylov
> vectors (if PC_LEFT)?
>
Everything in GMRES refers to the preconditioned operator. That is how
preconditioned GMRES works.
If you need the unpreconditioned space, you would have to use FGMRES, since
it pays the storage overhead.
Thanks,
Matt
> Cheers, Dave.
>
> On Mon, May 27, 2019 at 8:50 PM Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Mon, May 27, 2019 at 3:55 AM Dave Lee <davelee2804 at gmail.com> wrote:
>>
>>> 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?
>>>
>>
>> Left-preconditioning changes the operator, so you get he Arnoldi subspace
>> for the transforned operator, starting with a transformed rhs.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> 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/>
>>>>>
>>>>
>>
>> --
>> 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/6f062155/attachment-0001.html>
More information about the petsc-users
mailing list