[petsc-dev] SLEPc failure

Jose E. Roman jroman at dsic.upv.es
Thu Nov 2 04:03:54 CDT 2017 

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