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 <[email protected]> > 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" <[email protected]> >> À: "Jose E. Roman" <[email protected]> >> Cc: "For users of the development version of PETSc" <[email protected]> >> 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" <[email protected]> >>> À: "Franck Houssen" <[email protected]> >>> Cc: "For users of the development version of PETSc" <[email protected]> >>> 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 <[email protected]> >>>> 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" <[email protected]> >>>>> À: "Franck Houssen" <[email protected]> >>>>> Cc: "For users of the development version of PETSc" >>>>> <[email protected]> >>>>> 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 <[email protected]> >>>>>> 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" <[email protected]> >>>>>>> À: "Franck Houssen" <[email protected]> >>>>>>> Cc: "For users of the development version of PETSc" >>>>>>> <[email protected]> >>>>>>> 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 <[email protected]> >>>>>>>> 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)); >>>>>>> >>>>>>> >>>>> >>>>> >>> >>> >>
