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