[petsc-dev] QR factorization of dense matrix

Fri Dec 15 04:58:41 CST 2017

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));
>>>>>>> 
>>>>>>> 
>>>>> 
>>>>> 
>>> 
>>> 
>> 