[petsc-dev] QR factorization of dense matrix

[Jose E. Roman]

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