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.hous...@inria.fr> > À: "Jose E. Roman" <jro...@dsic.upv.es> > Cc: "For users of the development version of PETSc" <petsc-dev@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" <jro...@dsic.upv.es> > > À: "Franck Houssen" <franck.hous...@inria.fr> > > Cc: "For users of the development version of PETSc" <petsc-dev@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.hous...@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" <jro...@dsic.upv.es> > > >> À: "Franck Houssen" <franck.hous...@inria.fr> > > >> Cc: "For users of the development version of PETSc" > > >> <petsc-dev@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.hous...@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" <jro...@dsic.upv.es> > > >>>> À: "Franck Houssen" <franck.hous...@inria.fr> > > >>>> Cc: "For users of the development version of PETSc" > > >>>> <petsc-dev@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.hous...@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)); > > >>>> > > >>>> > > >> > > >> > > > > >