[petsc-dev] SLEPc failure

Franck Houssen franck.houssen at inria.fr
Fri Nov 3 06:41:25 CDT 2017 

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" <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é: 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.houssen at 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" <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 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.houssen at 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" <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 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.houssen at 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" <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 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.houssen at 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" <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 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.houssen at 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" <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é: Jeudi 26 Octobre 2017 18:49:22
> >>>>>>>>>> Objet: Re: [petsc-dev] SLEPc failure
> >>>>>>>>>> 
> >>>>>>>>>> 
> >>>>>>>>>>> El 26 oct 2017, a las 18:36, Franck Houssen
> >>>>>>>>>>> <franck.houssen at 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>
> >> 
> >> 
> 
>

More information about the petsc-dev mailing list