[petsc-dev] QR factorization of dense matrix

Fri Dec 15 07:55:18 CST 2017

OK. I have a MPI sparse matrix. I'll find a way to convert it to MPI dense. Thanks Jose.

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 15 Décembre 2017 11:58:41
> Objet: Re: [petsc-dev] QR factorization of dense matrix
> 
> It also works with MPIDense matrices. Use PETSC_COMM_WORLD only for bv and Z.
> R should still be sequential.
> Jose
> 
> > El 15 dic 2017, a las 11:47, Franck Houssen <franck.houssen at inria.fr>
> > escribió:
> > 
> > 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));
> >>>>>>> 
> >>>>>>> 
> >>>>> 
> >>>>> 
> >>> 
> >>> 
> >> 
> 
> 