bugGNU Octave - Bugs: bug #39573, Wrong result from eigs in...

 
 

bug #39573: Wrong result from eigs in generalized shift-invert mode

Submitted by:  None
Submitted on:  Fri 26 Jul 2013 01:25:53 PM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  3.6.4
Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Fri 17 Mar 2017 04:56:49 PM UTC, comment #8: 

I pushed the patch in bug #50546.  Nice work Marco, closing out a 4-year old bug.

Rik <rik5>
Project Administrator
Fri 17 Mar 2017 12:28:52 PM UTC, comment #7: 

Fixed in #50546.

Marco Caliari <caliari>
Project Member
Thu 16 Mar 2017 09:52:07 AM UTC, comment #6: 

Taking the input matrices from #50546, you see that Octave gives

error: __eigs__: eigs: error -9999 in dsaupd

for the fourth case (shift-invert mode). With the modifications described in comment #1, I get

error = 7.6926e-13

If there are no additional comments, I will make a patch.

(file #40014)

Marco Caliari <caliari>
Project Member
Wed 15 Mar 2017 12:12:23 PM UTC, comment #5: 

Octave now works simply because for rows (A) < 9, eig is called instead.

Marco Caliari <caliari>
Project Member
Wed 16 Nov 2016 10:44:07 PM UTC, comment #4: 

The original example here gives a different result with the current version of Octave

>> A = [1 0 1; 0 1 0; 1 0 0];
>> B = [1 0 0 ; 0 1 0 ; 0 0 0];
>> [V, D] = eigs (A, B, 1, 0);
>> residual_norm = norm (A*V - B*V*D)
residual_norm = 0

So at least something has changed with respect to this example since 3.6.4. If eigs is still not correct for some cases, please feel free to comment on this bug report or open a new bug report with details about what still is not working correctly.

Mike Miller <mtmiller>
Project Administrator
Tue 17 Sep 2013 07:59:47 AM UTC, comment #3: 

You are right, documentation requires positive semi-definite matrices but ARPACK performs quite badly.
For instance, matlab gives "error -9999 in dsaupd" with your example and the same would give Octave with the modification you suggested.

Marco

Marco Caliari <caliari>
Project Member
Mon 16 Sep 2013 09:45:51 AM UTC, comment #2: 

ARPACK documentation in DSAUPD.f states that for `mode = 3`,
the requirement is
"""
K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
"""
Positive definiteness is not required here.

However, what I'm not sure is whether this is actually correct, as ARPACK appears to behave quite badly for this class of problems.

Anonymous
Mon 05 Aug 2013 12:15:06 PM UTC, comment #1: 

Hi,

your example is not correct, since B is not positive definite. On the other hand, I think you are right and line 1275 (and 2238 and 3259) of eigs-base.cc should be changed to

vector_product (b, workd+iptr(0)-1, dtmp);

as for the case ido == 2.

Marco

Marco Caliari <caliari>
Project Member
Fri 26 Jul 2013 01:25:53 PM UTC, original submission:  

Consider the following problem:

A = [1 0 1; 0 1 0; 1 0 0 ]
B = [1 0 0 ; 0 1 0 ; 0 0 0 ]
[V,D]=eigs(A, B, 1, 0);
residual_norm=norm(A*V - B*V*D)

This gives large residual norms of order 1. The generalized eigenvalue problem is apparently solved incorrectly.

The reason appears to be that in `eigs-base.cc` the routine EigsRealSymmetricMatrixShift does the following for IDO=-1

1275 vector_product (m, workd+iptr(0)-1, dtmp);
1276
1277 Matrix tmp(n, 1);
1278
1279 for (octave_idx_type i = 0; i < n; i++)
1280 tmp(i,0) = dtmp[P[i]];
1281
1282 lusolve (L, U, tmp);

The ARPACK mode used is mode==3, for which the ARPACK DSAUPD documentation states

c  Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
c           ===> OP = (inv[K - sigma*M])*M  and  B = M.
c           ===> Shift-and-Invert mode
...
c          IDO = -1: compute  Y = OP * X  where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y.
c                    This is for the initialization phase to force the
c                    starting vector into the range of OP.

However, Octave is apparently computing here Y = inv[K - sigma*M]*K, as in the Octave routine,
`m` is apparently the matrix K, and `b` the matrix M. Is this correct?

The issue probably doesn't occur in typical usage as the IDO=-1 seems to be emitted by ARPACK only in the initial steps.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #40014:  eigenvalues.m added by caliari (11KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-03-17 rik5 StatusPatch Submitted => Fixed
        Open/ClosedOpen => Closed
    2017-03-17 caliari StatusIn Progress => Patch Submitted
        Open/ClosedClosed => Open
    2017-03-17 caliari StatusFixed => In Progress
    2017-03-16 caliari Attached File- => Added eigenvalues.m, #40014
    2016-11-16 mtmiller StatusNone => Fixed
        Open/ClosedOpen => Closed

    Back to the top


    Powered by Savane 3.5