[petsc-dev] QR factorization of dense matrix

Fri Dec 15 04:47:40 CST 2017

Coming back on that question: does BVOrthogonalize works also with distributed matrices ? Or only sequential ones ?

This works fine with sequential matrices :

Mat Z; MatCreateSeqDense(PETSC_COMM_SELF, n, m, NULL, &Z);
... // MatSetValues(Z, ...)
BVCreate(PETSC_COMM_SELF, &bv);
BVCreateFromMat(Z, &bv); // Z is tall-skinny
Mat R; MatCreateSeqDense(PETSC_COMM_SELF, m, m, NULL, &R);
BVOrthogonalize(bv, R);

If it could be possible to got this to work with distributed matrices, do I need to replace PETSC_COMM_SELF with PETSC_COMM_WORLD all over the place, or, only for Z and bv ? (the doc says R must be sequential dense)

Franck

----- Mail original -----
> De: "Franck Houssen" <franck.houssen at inria.fr>
> À: "Jose E. Roman" <jroman at dsic.upv.es>
> Cc: "For users of the development version of PETSc" <petsc-dev at mcs.anl.gov>
> Envoyé: Mardi 31 Octobre 2017 16:17:28
> Objet: Re: [petsc-dev] QR factorization of dense matrix
> 
> OK, got it. Thanks.
> 
> Franck
> 
> ----- Mail original -----
> > De: "Jose E. Roman" <jroman at dsic.upv.es>
> > À: "Franck Houssen" <franck.houssen at inria.fr>
> > Cc: "For users of the development version of PETSc" <petsc-dev at mcs.anl.gov>
> > Envoyé: Lundi 30 Octobre 2017 17:37:52
> > Objet: Re: [petsc-dev] QR factorization of dense matrix
> > 
> > Any BV type will do. The default BVSVEC is generally best.
> > Jose
> > 
> > 
> > > El 30 oct 2017, a las 17:18, Franck Houssen <franck.houssen at inria.fr>
> > > escribió:
> > > 
> > > It was not clear to me when I read the doc. That's OK now: got it to
> > > work,
> > > thanks Jose !
> > > Just to make sure, to make it work, I had to set a BV type: I chose BVMAT
> > > as I use BVCreateFromMat. Is that the good type ? (BVVECS works too)
> > > 
> > > Franck
> > > 
> > > ----- Mail original -----
> > >> De: "Jose E. Roman" <jroman at dsic.upv.es>
> > >> À: "Franck Houssen" <franck.houssen at inria.fr>
> > >> Cc: "For users of the development version of PETSc"
> > >> <petsc-dev at mcs.anl.gov>
> > >> Envoyé: Samedi 28 Octobre 2017 16:56:22
> > >> Objet: Re: [petsc-dev] QR factorization of dense matrix
> > >> 
> > >> Matrix R must be mxm.
> > >> BVOrthogonalize computes Z=Q*R, where Q overwrites Z.
> > >> Jose
> > >> 
> > >>> El 28 oct 2017, a las 13:11, Franck Houssen <franck.houssen at inria.fr>
> > >>> escribió:
> > >>> 
> > >>> I've seen that !... But can't get BVOrthogonalize to work.
> > >>> 
> > >>> I tried:
> > >>> Mat Z; MatCreateSeqDense(PETSC_COMM_SELF, n, m, NULL, &Z);
> > >>> ...; // MatSetValues(Z, ...)
> > >>> BVCreate(PETSC_COMM_SELF, &bv);
> > >>> BVCreateFromMat(Z, &bv); // Z is tall-skinny
> > >>> Mat R; MatCreateSeqDense(PETSC_COMM_SELF, n, m, NULL, &R); // Same n, m
> > >>> than Z.
> > >>> BVOrthogonalize(bv, R);
> > >>> 
> > >>> But BVOrthogonalize fails with :
> > >>>> [0]PETSC ERROR: Nonconforming object sizes
> > >>>> [0]PETSC ERROR: Mat argument is not square, it has 1 rows and 3
> > >>>> columns
> > >>> 
> > >>> So, as I didn't get what's wrong, I was looking for another way to do
> > >>> this.
> > >>> 
> > >>> Franck
> > >>> 
> > >>> ----- Mail original -----
> > >>>> De: "Jose E. Roman" <jroman at dsic.upv.es>
> > >>>> À: "Franck Houssen" <franck.houssen at inria.fr>
> > >>>> Cc: "For users of the development version of PETSc"
> > >>>> <petsc-dev at mcs.anl.gov>
> > >>>> Envoyé: Vendredi 27 Octobre 2017 19:03:37
> > >>>> Objet: Re: [petsc-dev] QR factorization of dense matrix
> > >>>> 
> > >>>> Franck,
> > >>>> 
> > >>>> SLEPc has some support for this, but it is intended only for
> > >>>> tall-skinny
> > >>>> matrices, that is, when the number of columns is much smaller than
> > >>>> rows.
> > >>>> For
> > >>>> an almost square matrix you should not use it.
> > >>>> 
> > >>>> Have a look at this
> > >>>> http://slepc.upv.es/documentation/current/docs/manualpages/BV/BVOrthogonalize.html
> > >>>> http://slepc.upv.es/documentation/current/docs/manualpages/BV/BVOrthogBlockType.html
> > >>>> 
> > >>>> You can see there are three methods. All of them have drawbacks:
> > >>>> GS: This is a Gram-Schmidt QR, computed column by column, so it is
> > >>>> slower
> > >>>> than the other two. However, it is robust.
> > >>>> CHOL: Cholesky QR, it is not numerically stable. In the future we will
> > >>>> add
> > >>>> Cholesky QR2.
> > >>>> TSQR: Unfortunately this is not implemented in parallel. I wanted to
> > >>>> add
> > >>>> the
> > >>>> parallel version for 3.8, but didn't have time. It will be added soon.
> > >>>> 
> > >>>> You can use BVCreateFromMat() to create a BV object from a Mat.
> > >>>> 
> > >>>> Jose
> > >>>> 
> > >>>> 
> > >>>>> El 27 oct 2017, a las 18:39, Franck Houssen <franck.houssen at inria.fr>
> > >>>>> escribió:
> > >>>>> 
> > >>>>> I am looking for QR factorization of (sequential) dense matrix: is
> > >>>>> this
> > >>>>> available in PETSc ? I "just" need the diagonal of R (I do not need
> > >>>>> neither the full content of R, nor Q)
> > >>>>> 
> > >>>>> I found that (old !) thread
> > >>>>> https://lists.mcs.anl.gov/pipermail/petsc-users/2013-November/019577.html
> > >>>>> that says it could be implemented: has it been done ?
> > >>>>> As for a direct solve, the way to go is "KSPSetType(ksp, KSPPREONLY);
> > >>>>> PCSetType(pc, PCLU);", I was expecting something like
> > >>>>> "KSPSetType(ksp,
> > >>>>> KSPPREONLY); PCSetType(pc, PCQR);"... But it seems there is no PCQR
> > >>>>> available. Or is it possible to do that using "an iterative way" with
> > >>>>> a
> > >>>>> specific kind of KSP that triggers a Gram Schmidt orthogonalization
> > >>>>> in
> > >>>>> back-end ? (I have seen a KSPLSQR but could I get Q and R back ? As I
> > >>>>> understand this, I would say no: I would say the user can only get
> > >>>>> the
> > >>>>> solution)
> > >>>>> 
> > >>>>> Is it possible to QR a (sequential) dense matrix in PETSc ? If yes,
> > >>>>> what
> > >>>>> are the steps to follow ?
> > >>>>> 
> > >>>>> Franck
> > >>>>> 
> > >>>>> My understanding is that DGEQRF from lapack can do "more" than what I
> > >>>>> need,
> > >>>>> but, no sure to get if I can use it from PETSc through a KSP:
> > >>>>>>> git grep DGEQRF
> > >>>>> include/petscblaslapack_stdcall.h:#  define LAPACKgeqrf_ DGEQRF
> > >>>>>>> git grep LAPACKgeqrf_
> > >>>>> include/petscblaslapack.h:PETSC_EXTERN void
> > >>>>> LAPACKgeqrf_(PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
> > >>>>> include/petscblaslapack_mangle.h:#define LAPACKgeqrf_
> > >>>>> PETSCBLAS(geqrf,GEQRF)
> > >>>>> include/petscblaslapack_stdcall.h:#  define LAPACKgeqrf_ SGEQRF
> > >>>>> include/petscblaslapack_stdcall.h:#  define LAPACKgeqrf_ DGEQRF
> > >>>>> include/petscblaslapack_stdcall.h:#  define LAPACKgeqrf_ CGEQRF
> > >>>>> include/petscblaslapack_stdcall.h:#  define LAPACKgeqrf_ ZGEQRF
> > >>>>> include/petscblaslapack_stdcall.h:PETSC_EXTERN void PETSC_STDCALL
> > >>>>> LAPACKgeqrf_(PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
> > >>>>> src/dm/dt/interface/dt.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
> > >>>>> src/dm/dt/interface/dtfv.c:
> > >>>>> LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
> > >>>>> src/ksp/ksp/impls/gmres/agmres/agmres.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&lC, &KspSize,
> > >>>>> agmres->hh_origin, &ldH, agmres->tau, agmres->work, &lwork, &info));
> > >>>>> src/ksp/pc/impls/bddc/bddcprivate.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,&lqr_work_t,&lqr_work,&lierr));
> > >>>>> src/ksp/pc/impls/bddc/bddcprivate.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,qr_work,&lqr_work,&lierr));
> > >>>>> src/ksp/pc/impls/gamg/agg.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA,
> > >>>>> TAU, WORK, &LWORK, &INFO));
> > >>>>> src/tao/leastsquares/impls/pounders/pounders.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&blasnp,&blasnplus1,mfqP->Q_tmp,&blasnpmax,mfqP->tau_tmp,mfqP->mwork,&blasmaxmn,&info));
> > >>>>> src/tao/leastsquares/impls/pounders/pounders.c:
> > >>>>> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->mwork,&blasmaxmn,&info));
> > >>>> 
> > >>>> 
> > >> 
> > >> 
> > 
> > 
> 