OK, thanks. Unfortunately, this does not help. So I need to find a another way to extract B from the data set in order to get a "better" (= "less" singular) B. Thanks for tips and help.
Franck ----- Mail original ----- > De: "Jose E. Roman" <jro...@dsic.upv.es> > À: "Franck Houssen" <franck.hous...@inria.fr> > Cc: "For users of the development version of PETSc" <petsc-dev@mcs.anl.gov> > Envoyé: Jeudi 2 Novembre 2017 10:03:54 > Objet: Re: [petsc-dev] SLEPc failure > > Could you please try the following modified function? It should replace the > one in $SLEPC_DIR/include/slepc/private/bvimpl.h > Thanks. > > PETSC_STATIC_INLINE PetscErrorCode BV_SafeSqrt(BV bv,PetscScalar > alpha,PetscReal *res) > { > PetscErrorCode ierr; > PetscReal absal,realp; > > PetscFunctionBegin; > absal = PetscAbsScalar(alpha); > realp = PetscRealPart(alpha); > if (absal<PETSC_MACHINE_EPSILON) { > ierr = PetscInfo(bv,"Zero norm, either the vector is zero or a semi-inner > product is being used\n");CHKERRQ(ierr); > } > #if defined(PETSC_USE_COMPLEX) > if (PetscAbsReal(PetscImaginaryPart(alpha))>PETSC_MACHINE_EPSILON && > PetscAbsReal(PetscImaginaryPart(alpha))/absal>100*PETSC_MACHINE_EPSILON) > SETERRQ1(PetscObjectComm((PetscObject)bv),1,"The inner product is not well > defined: nonzero imaginary part %g",PetscImaginaryPart(alpha)); > #endif > if (bv->indef) { > *res = (realp<0.0)? -PetscSqrtReal(-realp): PetscSqrtReal(realp); > } else { > if (realp<-10*PETSC_MACHINE_EPSILON) > SETERRQ(PetscObjectComm((PetscObject)bv),1,"The inner product is not > well defined: indefinite matrix"); > *res = (realp<0.0)? 0.0: PetscSqrtReal(realp); > } > PetscFunctionReturn(0); > } > > > > > El 31 oct 2017, a las 16:17, Franck Houssen <franck.hous...@inria.fr> > > escribió: > > > > Thanks, this is helpful ! At least, I have some clues to dig deeper into. > > > > The data set I have seems to be very sensitive: I knew B would likely be > > close to singular (even arpack, which seems to be the most stable, is not > > always robust enough depending on use cases and/or number of domains). To > > get things stable with arpack, I had to add "-mat_mumps_cntl_1 0.01 > > -mat_mumps_cntl_3 -0.0001 -mat_mumps_cntl_4 0.0001": so I let all that > > also when testing with krylovschur. I've just tried to suppress these > > extra options: I get the MUMPS error in numerical factorization. At least, > > it make sense... > > > > Now, I added -eps_gen_non_hermitian and failed with EPS_DIVERGED_ITS. > > Increasing -eps_max_it does not help. So, I guess this is the end of the > > road. > > > > Would you consider to expose (in future release) the tolerance of this > > check ? Or is this something you really want to keep private ? (whatever B > > is or not singular - I guess in my case, this would have not helped > > anyway) > > > > Franck > > > > ----- Mail original ----- > >> De: "Jose E. Roman" <jro...@dsic.upv.es> > >> À: "Franck Houssen" <franck.hous...@inria.fr> > >> Cc: "For users of the development version of PETSc" > >> <petsc-dev@mcs.anl.gov> > >> Envoyé: Lundi 30 Octobre 2017 18:11:48 > >> Objet: Re: [petsc-dev] SLEPc failure > >> > >> I am getting a MUMPS error in numerical factorization... > >> Anyway, your B matrix is singular, with a high-dimensional nullspace. > >> Maybe > >> this is producing small negative values when computing v'*B*v. > >> There is no way to relax the check. You should solve the problem as > >> non-symmetric. Or use Arpack if it works for you. > >> > >> Jose > >> > >> > >> > >>> El 30 oct 2017, a las 17:16, Franck Houssen <franck.hous...@inria.fr> > >>> escribió: > >>> > >>> I deal with domain decomposition. It was faster/easier to generate files > >>> with 1 proc: I guess the "dirichlet" and "neumann" matrices are the same > >>> in this case, so one get the same files in the end... I didn't realized > >>> that when I sent the files. My mistake. > >>> > >>> At my side, when I use krylovschur, SLEPc fails using 1, 2, 4, 8, ... > >>> procs > >>> (each proc performs a SLEPc solve that may fail or not - difficult to > >>> catch one). For instance, I attached data of one failed domain out of 8 > >>> (8 > >>> MPI procs): matrices are very close but different. Moreover, I added EPS > >>> logs of the same run and domain but replacing krylovschur with arpack > >>> when > >>> SLEPc does not fail (regarding your remark that B could be indefinite, I > >>> added -mat_mumps_icntl_33 1 to get the determinant). > >>> > >>> Anyway, I don't expect you spend too much time on this. My understanding > >>> is > >>> that there is no way to relax this check ? Correct ? > >>> > >>> Franck > >>> > >>> ----- Mail original ----- > >>>> De: "Jose E. Roman" <jro...@dsic.upv.es> > >>>> À: "Franck Houssen" <franck.hous...@inria.fr> > >>>> Cc: "For users of the development version of PETSc" > >>>> <petsc-dev@mcs.anl.gov> > >>>> Envoyé: Samedi 28 Octobre 2017 17:40:41 > >>>> Objet: Re: [petsc-dev] SLEPc failure > >>>> > >>>> The two matrices are the same ... > >>>> > >>>>> El 28 oct 2017, a las 13:11, Franck Houssen <franck.hous...@inria.fr> > >>>>> escribió: > >>>>> > >>>>> I just added that before EPSSetOperators: > >>>>> PetscViewer viewerA; > >>>>> PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Atau.out",FILE_MODE_WRITE,&viewerA); > >>>>> MatView(A,viewerA); > >>>>> PetscViewer viewerB; > >>>>> PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Btau.out",FILE_MODE_WRITE,&viewerB); > >>>>> MatView(B,viewerB); > >>>>> > >>>>> At first, I avoided binary as I didn't know if the format handles > >>>>> big/little endianness... So I prefered ASCII. > >>>>> > >>>>> Binary data are attached. > >>>>> > >>>>> Franck > >>>>> > >>>>> PS : running debian with little endian. > >>>>>>> python -c "import sys;print(0 if sys.byteorder=='big' else 1)" > >>>>> 1 > >>>>> > >>>>> > >>>>> ----- Mail original ----- > >>>>>> De: "Jose E. Roman" <jro...@dsic.upv.es> > >>>>>> À: "Franck Houssen" <franck.hous...@inria.fr> > >>>>>> Cc: "For users of the development version of PETSc" > >>>>>> <petsc-dev@mcs.anl.gov> > >>>>>> Envoyé: Vendredi 27 Octobre 2017 18:52:56 > >>>>>> Objet: Re: [petsc-dev] SLEPc failure > >>>>>> > >>>>>> 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 <franck.hous...@inria.fr> > >>>>>>> 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" <jro...@dsic.upv.es> > >>>>>>>> À: "Franck Houssen" <franck.hous...@inria.fr> > >>>>>>>> Cc: "For users of the development version of PETSc" > >>>>>>>> <petsc-dev@mcs.anl.gov> > >>>>>>>> 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 > >>>>>>>>> <franck.hous...@inria.fr> > >>>>>>>>> escribió: > >>>>>>>>> > >>>>>>>>> I use the development version (bitbucket clone). How to relax the > >>>>>>>>> check > >>>>>>>>> ? > >>>>>>>>> At command line option ? > >>>>>>>>> > >>>>>>>>> Franck > >>>>>>>>> > >>>>>>>>> ----- Mail original ----- > >>>>>>>>>> De: "Jose E. Roman" <jro...@dsic.upv.es> > >>>>>>>>>> À: "Franck Houssen" <franck.hous...@inria.fr> > >>>>>>>>>> Cc: "For users of the development version of PETSc" > >>>>>>>>>> <petsc-dev@mcs.anl.gov> > >>>>>>>>>> Envoyé: Jeudi 26 Octobre 2017 18:49:22 > >>>>>>>>>> Objet: Re: [petsc-dev] SLEPc failure > >>>>>>>>>> > >>>>>>>>>> > >>>>>>>>>>> El 26 oct 2017, a las 18:36, Franck Houssen > >>>>>>>>>>> <franck.hous...@inria.fr> > >>>>>>>>>>> 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 > >>>>>>>>>> > >>>>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>> <ABtau.tar.gz> > >>>>>> > >>>>>> > >>>>> <ABtau.tar.gz> > >>>> > >>>> > >>> <slepc.tar.gz> > >> > >> > >