[petsc-dev] SLEPc General Eigenvalue Shift-Invert Crash
Miguel Arriaga
mta2122 at columbia.edu
Mon Dec 17 13:39:27 CST 2012
Thanks for the reply,
I had tried changing the target before but the problem is that the
eigenvalues will change depending on my target. I'm using Fortran and
my code looks like this:
call EPSCreate(PETSC_COMM_WORLD,eps,ierr)
(MatAssembly stuff here)
call EPSSetOperators(eps,Kmat,Mmat,ierr)
call EPSSetProblemType(eps,EPS_GNHEP,ierr)
call EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL,ierr)
call EPSSetTolerances(eps,tol,maxit,ierr)
call EPSSetDimensions(eps,nev,PETSC_DECIDE,PETSC_DECIDE,ierr)
call EPSSetFromOptions(eps,ierr)
call EPSSolve(eps,ierr)
Since I don't know where my eigenvalues will be (because they change
as the material yields) I'm a little concerned with using target.
Maybe I'm doing something wrong? I noticed that although I set
EPS_LARGEST_REAL in my code, if I set it again in the command line it
will change the output. Another thing is that I'm using Shift-Invert
to look at a Largest Real problem, which although it seemed odd to me,
it showed a much faster convergence.
Below I show some test results I obtained for just the first time step.
_Results:
-I'm asking for 6 eigenvalues and I have a tol of 1E-5
-Notice that the eigenvalues are different from the purely imaginary I
got from Matlab and SLEPc without SInvert.
BASE = -st_type sinvert -st_ksp_type preonly -st_pc_type lu
-st_pc_factor_mat_solver_package umfpack (has the mentioned error)
BASE -eps_target 0.2
Global eig 0 Real 8.08057606271894713E-007 Im
0.0000000000000000
Global eig 1 Real 5.88946613699903310E-010 Im
1.11478341210332509E-005
Global eig 2 Real 5.88946613699903310E-010 Im
-1.11478341210332509E-005
Global eig 3 Real 1.00727620688800812E-010 Im
1.60248025360939258E-010
Global eig 4 Real 1.00727620688800812E-010 Im
-1.60248025360939258E-010
Global eig 5 Real -2.82380008265192828E-010 Im
1.27051563295447506E-010
Global eig 6 Real -2.82380008265192828E-010 Im
-1.27051563295447506E-010
Global eig 7 Real -45.307348017077182 Im
0.0000000000000000
BASE -eps_target 1000
Global eig 0 Real 1.33590310724684969E-006 Im
5.69992625689008225E-003
Global eig 1 Real 1.33590310724684969E-006 Im
-5.69992625689008225E-003
Global eig 2 Real -45.307261982302634 Im
0.0000000000000000
Global eig 3 Real -90.614531192221420 Im
0.0000000000000000
Global eig 4 Real -182.34960882354494 Im
0.0000000000000000
Global eig 5 Real -227.65688316119918 Im 0.0000000000000000
BASE -eps_target 10000
Global eig 0 Real 0.43857832608046010 Im
0.0000000000000000
Global eig 1 Real -0.44680019931001880 Im
0.0000000000000000
Global eig 2 Real -45.307619734525360 Im
0.0000000000000000
Global eig 3 Real -90.609456729185695 Im
0.0000000000000000
Global eig 4 Real -1.23910373804392293E-003 Im
1842.6494652654314
Global eig 5 Real -1.23910373804392293E-003 Im -1842.6494652654314
BASE -eps_target 0.2 -eps_largest_real
Global eig 0 Real 5.88925547218011047E-010 Im
1.11479604597693868E-005
Global eig 1 Real 5.88925547218011047E-010 Im
-1.11479604597693868E-005
Global eig 2 Real 3.41273675985576119E-011 Im
2.48751366232282650E-011
Global eig 3 Real 3.41273675985576119E-011 Im
-2.48751366232282650E-011
Global eig 4 Real 1.94852467494399662E-012 Im
5.58171840043041153E-011
Global eig 5 Real 1.94852467494399662E-012 Im
-5.58171840043041153E-011
Global eig 6 Real -3.06560610230377506E-011 Im
2.38634932296064250E-011
Global eig 7 Real -3.06560610230377506E-011 Im
-2.38634932296064250E-011
Global eig 8 Real -3.14280823587864688E-010 Im 0.0000000000000000
BASE -eps_target 1000 -eps_largest_real
Global eig 0 Real 3.41060513164848089E-012 Im
1.32716666241605970E-005
Global eig 1 Real 3.41060513164848089E-012 Im
-1.32716666241605970E-005
Global eig 2 Real -3.36171979142818600E-010 Im
1842.6504379364101
Global eig 3 Real -3.36171979142818600E-010 Im
-1842.6504379364101
Global eig 4 Real -2.93482571578351781E-009 Im
2008.6622800422490
Global eig 5 Real -2.93482571578351781E-009 Im -2008.6622800422490
BASE -eps_target 10000 -eps_largest_real
Global eig 0 Real 8.44766719455947168E-003 Im
0.23198091192374859
Global eig 1 Real 8.44766719455947168E-003 Im
-0.23198091192374859
Global eig 2 Real 2.18278728425502777E-011 Im
1842.6504379364712
Global eig 3 Real 2.18278728425502777E-011 Im
-1842.6504379364712
Global eig 4 Real -2.00088834390044212E-011 Im
2008.6622800484372
Global eig 5 Real -2.00088834390044212E-011 Im -2008.6622800484372
Thank you in advance for your help,
Miguel
On Mon, Dec 17, 2012 at 5:08 AM, Jose E. Roman <jroman at dsic.upv.es> wrote:
> This may be caused by the KSP coefficient matrix being singular. With shift-and-invert the KSP matrix is A-sigma*B where the shift sigma is equal to the target. Are you specifying a target? By default the target in shift-and-invert runs is 0 so in that case the KSP matrix is A. Probably your A matrix is singular or has tiny eigenvalues. I would suggest trying with a nonzero target, e.g. -eps_target 0.2 (or a larger value if you know that your unstable eigenvalues are further away from the real axis).
>
> If this does not solve the problem, send configure and make logs to slepc-maint, together with the output of -eps_view
>
> Jose
>
>
> El 17/12/2012, a las 06:52, Barry Smith escribió:
>
>>
>> A google of DHSEQR 25 led to this page: http://www.netlib.org/lapack/explore-html/d8/d66/dhseqr_8f.html
>>
>> INFO is INTEGER
>> = 0: successful exit
>> .LT. 0: if INFO = -i, the i-th argument had an illegal
>> value
>> .GT. 0: if INFO = i, DHSEQR failed to compute all of
>> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
>> and WI contain those eigenvalues which have been
>> successfully computed. (Failures are rare.)
>>
>> If INFO .GT. 0 and JOB = 'E', then on exit, the
>> remaining unconverged eigenvalues are the eigen-
>> values of the upper Hessenberg matrix rows and
>> columns ILO through INFO of the final, output
>> value of H.
>>
>>
>> Being at the end of my knowledge, I'll suggest perhaps you can run everything in __float128 precision? You need a non-ancient version of GCC and
>> configure PETSc using -with-precision=__float128 --download-f2cblaslapack (I hope you don't have Fortran code cause I don't know if that will build properly with __float28). Also unlikely umfpack or any other external direct solver package supports __float128?
>>
>> Hopefully the eigenvalue experts will have a better solution.
>>
>> Barry
>>
>>
>>
>> On Dec 16, 2012, at 11:16 PM, Miguel Arriaga <mta2122 at columbia.edu> wrote:
>>
>>> Hello.
>>>
>>> I'm using slepc and I'm trying to do a stability analysis on a solid
>>> mechanics problem which in my case is associated with looking at the
>>> highest real part of the eigenvalues of the general problem Ax=λBx. My
>>> problem has a real non-symmetric A and a positive definite B. When I
>>> try to ask for the largest real eigenvalues it takes a great number of
>>> iterations to solve the problem, especially when I reach the yielding
>>> phase of my problem, in which case I can no longer obtain the required
>>> eigenvalues.
>>>
>>> The eigenvalues on the first steps are similar to what I obtained in
>>> Matlab for largest magnitude (Matlab doesn't converge for the largest real
>>> option). These eigenvalues are purely imaginary. When I was doing a
>>> standard eigenvalue analysis on this problem I found that I could
>>> greatly increase my performance with a Shift-Invert transformation.
>>> However, when I try the SI in the general case I get the following
>>> error:
>>>
>>> [0]PETSC ERROR: Error in external library!
>>> [0]PETSC ERROR: Error in Lapack xHSEQR 25!
>>>
>>> I'm using umfpack to do a preonly LU solve. Preconditioned
>>> Eigensolvers seem not to converge at all so for now I'm stuck with
>>> these. When I removed the umfpack I got the following error:
>>>
>>> [0]PETSC ERROR: Detected zero pivot in LU factorization:
>>> see http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot!
>>> [0]PETSC ERROR: Zero pivot row 0 value 0 tolerance 2.22045e-14!
>>>
>>> Based on the FAQ I tried using -st_pc_factor_shift_type
>>> POSITIVE_DEFINITE. While this made sinvert not crash and converge fast
>>> it also messed with my eigenvalues.
>>>
>>> Any ideas on how I could solve this?
>>>
>>> Thank you so much,
>>> Miguel Arriaga
>>
>
More information about the petsc-dev
mailing list