Re: [petsc-dev] QR factorization of dense matrix

2017-10-27 Thread Jose E. Roman 
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  
> 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_(,,A,,tau,work,,));
> src/dm/dt/interface/dtfv.c:  
> LAPACKgeqrf_(,,A,,tau,work,,);
> src/ksp/ksp/impls/gmres/agmres/agmres.c:  
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(, , 
> agmres->hh_origin, , agmres->tau, agmres->work, , ));
> src/ksp/pc/impls/bddc/bddcprivate.c:  
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(_M,_N,qr_basis,_LDA,qr_tau,_work_t,_work,));
> src/ksp/pc/impls/bddc/bddcprivate.c:  
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(_M,_N,qr_basis,_LDA,qr_tau,qr_work,_work,));
> src/ksp/pc/impls/gamg/agg.c:  
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(, , qqc, , TAU, 
> WORK, , ));
> src/tao/leastsquares/impls/pounders/pounders.c:
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(,,mfqP->Q_tmp,,mfqP->tau_tmp,mfqP->mwork,,));
> src/tao/leastsquares/impls/pounders/pounders.c:
> PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(,,mfqP->Q,,mfqP->tau,mfqP->mwork,,));

Re: [petsc-dev] SLEPc failure

2017-10-27 Thread Jose E. Roman 
I cannot load the files you sent. Please send the matrices in binary format. 
The easiest way is to run your program with -eps_view_mat0 binary:Atau.bin 
-eps_view_mat1 binary:Btau.bin

However, the files are written at the end of EPSSolve() so if the solve fails 
then it will not create the files. You can try running with -eps_max_it 1  or 
add code in your main program to write the matrices.

Jose


> El 27 oct 2017, a las 12:28, Franck Houssen  
> escribió:
> 
> Maybe could be convenient for the users to have an option (or an EPSSetXXX) 
> to relax that check ?
> Data are attached.
> 
> Franck 
> 
> - Mail original -
>> De: "Jose E. Roman" 
>> À: "Franck Houssen" 
>> Cc: "For users of the development version of PETSc" 
>> Envoyé: Vendredi 27 Octobre 2017 10:15:44
>> Objet: Re: [petsc-dev] SLEPc failure
>> 
>> There is no new option. What I mean is that from 3.7 to 3.8 we changed the
>> line that produces this error. But it seems that it is still failing in your
>> problem. Maybe your B matrix is indefinite or not exactly symmetric. Can you
>> send me the matrices?
>> Jose
>> 
>>> El 27 oct 2017, a las 9:57, Franck Houssen 
>>> escribió:
>>> 
>>> I use the development version (bitbucket clone). How to relax the check ?
>>> At command line option ?
>>> 
>>> Franck
>>> 
>>> - Mail original -
 De: "Jose E. Roman" 
 À: "Franck Houssen" 
 Cc: "For users of the development version of PETSc"
 
 Envoyé: Jeudi 26 Octobre 2017 18:49:22
 Objet: Re: [petsc-dev] SLEPc failure
 
 
> El 26 oct 2017, a las 18:36, Franck Houssen 
> escribió:
> 
> Here is a stack I end up with when trying to solve an eigen problem
> (real,
> sym, generalized) with SLEPc. My understanding is that, during the Gram
> Schmidt orthogonalisation, the projection of one basis vector turns out
> to
> be null.
> First, is this correct ? Second, in such cases, are there some
> recommended
> "recipe" to test/try (options) to get a clue on the problem ? (I would
> unfortunately perfectly understand the answer could be no !... As this
> totally depends on A/B).
> 
> With arpack, the eigen problem is solved (so the matrix A and B I use
> seems
> to be relevant). But, when I change from arpack to
> krylovschur/ciss/arnoldi, I get the stack below.
> 
> Franck
> 
> [0]PETSC ERROR: #1 BV_SafeSqrt()
> [0]PETSC ERROR: #2 BVNorm_Private()
> [0]PETSC ERROR: #3 BVNormColumn()
> [0]PETSC ERROR: #4 BV_NormVecOrColumn()
> [0]PETSC ERROR: #5 BVOrthogonalizeCGS1()
> [0]PETSC ERROR: #6 BVOrthogonalizeGS()
> [0]PETSC ERROR: #7 BVOrthonormalizeColumn()
> [0]PETSC ERROR: #8 EPSFullLanczos()
> [0]PETSC ERROR: #9 EPSSolve_KrylovSchur_Symm()
> [0]PETSC ERROR: #10 EPSSolve()
 
 Is this with SLEPc 3.8? In SLEPc 3.8 we relaxed this check so I would
 suggest
 trying with it.
 Jose
 
 
>> 
>> 
>

[petsc-dev] QR factorization of dense matrix

2017-10-27 Thread Franck Houssen 
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_(,,A,,tau,work,,));
 
src/dm/dt/interface/dtfv.c: LAPACKgeqrf_(,,A,,tau,work,,); 
src/ksp/ksp/impls/gmres/agmres/agmres.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(, , agmres->hh_origin, 
, agmres->tau, agmres->work, , )); 
src/ksp/pc/impls/bddc/bddcprivate.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(_M,_N,qr_basis,_LDA,qr_tau,_work_t,_work,));
 
src/ksp/pc/impls/bddc/bddcprivate.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(_M,_N,qr_basis,_LDA,qr_tau,qr_work,_work,));
 
src/ksp/pc/impls/gamg/agg.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(, , qqc, , TAU, WORK, 
, )); 
src/tao/leastsquares/impls/pounders/pounders.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(,,mfqP->Q_tmp,,mfqP->tau_tmp,mfqP->mwork,,));
 
src/tao/leastsquares/impls/pounders/pounders.c: 
PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(,,mfqP->Q,,mfqP->tau,mfqP->mwork,,));

Re: [petsc-dev] SLEPc failure

2017-10-27 Thread Jose E. Roman 
There is no new option. What I mean is that from 3.7 to 3.8 we changed the line 
that produces this error. But it seems that it is still failing in your 
problem. Maybe your B matrix is indefinite or not exactly symmetric. Can you 
send me the matrices?
Jose

> El 27 oct 2017, a las 9:57, Franck Houssen  escribió:
> 
> I use the development version (bitbucket clone). How to relax the check ? At 
> command line option ?
> 
> Franck
> 
> - Mail original -
>> De: "Jose E. Roman" 
>> À: "Franck Houssen" 
>> Cc: "For users of the development version of PETSc" 
>> Envoyé: Jeudi 26 Octobre 2017 18:49:22
>> Objet: Re: [petsc-dev] SLEPc failure
>> 
>> 
>>> El 26 oct 2017, a las 18:36, Franck Houssen 
>>> escribió:
>>> 
>>> Here is a stack I end up with when trying to solve an eigen problem (real,
>>> sym, generalized) with SLEPc. My understanding is that, during the Gram
>>> Schmidt orthogonalisation, the projection of one basis vector turns out to
>>> be null.
>>> First, is this correct ? Second, in such cases, are there some recommended
>>> "recipe" to test/try (options) to get a clue on the problem ? (I would
>>> unfortunately perfectly understand the answer could be no !... As this
>>> totally depends on A/B).
>>> 
>>> With arpack, the eigen problem is solved (so the matrix A and B I use seems
>>> to be relevant). But, when I change from arpack to
>>> krylovschur/ciss/arnoldi, I get the stack below.
>>> 
>>> Franck
>>> 
>>> [0]PETSC ERROR: #1 BV_SafeSqrt()
>>> [0]PETSC ERROR: #2 BVNorm_Private()
>>> [0]PETSC ERROR: #3 BVNormColumn()
>>> [0]PETSC ERROR: #4 BV_NormVecOrColumn()
>>> [0]PETSC ERROR: #5 BVOrthogonalizeCGS1()
>>> [0]PETSC ERROR: #6 BVOrthogonalizeGS()
>>> [0]PETSC ERROR: #7 BVOrthonormalizeColumn()
>>> [0]PETSC ERROR: #8 EPSFullLanczos()
>>> [0]PETSC ERROR: #9 EPSSolve_KrylovSchur_Symm()
>>> [0]PETSC ERROR: #10 EPSSolve()
>> 
>> Is this with SLEPc 3.8? In SLEPc 3.8 we relaxed this check so I would suggest
>> trying with it.
>> Jose
>> 
>>

Re: [petsc-dev] SLEPc failure

2017-10-27 Thread Franck Houssen 
I use the development version (bitbucket clone). How to relax the check ? At 
command line option ?

Franck

- Mail original -
> De: "Jose E. Roman" 
> À: "Franck Houssen" 
> Cc: "For users of the development version of PETSc" 
> Envoyé: Jeudi 26 Octobre 2017 18:49:22
> Objet: Re: [petsc-dev] SLEPc failure
> 
> 
> > El 26 oct 2017, a las 18:36, Franck Houssen 
> > escribió:
> > 
> > Here is a stack I end up with when trying to solve an eigen problem (real,
> > sym, generalized) with SLEPc. My understanding is that, during the Gram
> > Schmidt orthogonalisation, the projection of one basis vector turns out to
> > be null.
> > First, is this correct ? Second, in such cases, are there some recommended
> > "recipe" to test/try (options) to get a clue on the problem ? (I would
> > unfortunately perfectly understand the answer could be no !... As this
> > totally depends on A/B).
> > 
> > With arpack, the eigen problem is solved (so the matrix A and B I use seems
> > to be relevant). But, when I change from arpack to
> > krylovschur/ciss/arnoldi, I get the stack below.
> > 
> > Franck
> > 
> > [0]PETSC ERROR: #1 BV_SafeSqrt()
> > [0]PETSC ERROR: #2 BVNorm_Private()
> > [0]PETSC ERROR: #3 BVNormColumn()
> > [0]PETSC ERROR: #4 BV_NormVecOrColumn()
> > [0]PETSC ERROR: #5 BVOrthogonalizeCGS1()
> > [0]PETSC ERROR: #6 BVOrthogonalizeGS()
> > [0]PETSC ERROR: #7 BVOrthonormalizeColumn()
> > [0]PETSC ERROR: #8 EPSFullLanczos()
> > [0]PETSC ERROR: #9 EPSSolve_KrylovSchur_Symm()
> > [0]PETSC ERROR: #10 EPSSolve()
> 
> Is this with SLEPc 3.8? In SLEPc 3.8 we relaxed this check so I would suggest
> trying with it.
> Jose
> 
>