[Octave-bug-tracker] [bug #39573] Wrong result from eigs in generalized shift-invert mode

[anonymous]

URL:
  <http://savannah.gnu.org/bugs/?39573>

                 Summary: Wrong result from eigs in generalized shift-invert
mode
                 Project: GNU Octave
            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: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.6.4
        Operating System: Any

    _______________________________________________________

Details:

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.





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?39573>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/