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

 
 

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

Submitter:  None
Submitted:  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 Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
Group administrator
Fri 17 Mar 2017 12:28:52 PM UTC, comment #7: 

Fixed in #50546.

Marco Caliari <caliari>
Group 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>
Group 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>
Group 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>
Group Member
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>
Group 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>
Group 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)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group 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.13-02a9.
    Corresponding source code