# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1489749940 -3600
#      Fri Mar 17 12:25:40 2017 +0100
# Branch stable
# Node ID dcf30fdacb273d190669b78cca232b73386c2ead
# Parent  3b2dbed2676282ac9e5e0007ec8b747c11ddf606
Fix eigs for generalized nonsymmetric and shift-invert problems (bug #39573)

* liboctave/numeric/eigs-base.cc: fix the operator which computes mtmp in
EigsRealNonSymmetricMatrix and EigsComplexNonSymmetricMatrix. Fix
vector_product for ido=-1.

* scripts/sparse/eigs.m: add tests for the generalized nonsymmetric
problem (nonsymmetric and complex) and for the generalized shift-invert
problems (symmetric, nonsymmetric and complex).

diff -r 3b2dbed26762 -r dcf30fdacb27 liboctave/numeric/eigs-base.cc
--- a/liboctave/numeric/eigs-base.cc	Thu Mar 16 09:44:44 2017 +0100
+++ b/liboctave/numeric/eigs-base.cc	Fri Mar 17 12:25:40 2017 +0100
@@ -1022,7 +1022,7 @@ EigsRealSymmetricMatrixShift (const M& m
                 {
                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
 
-                  vector_product (m, workd+iptr(0)-1, dtmp);
+                  vector_product (b, workd+iptr(0)-1, dtmp);
 
                   Matrix tmp(n, 1);
 
@@ -1566,7 +1566,7 @@ EigsRealNonSymmetricMatrix (const M& m, 
               for (octave_idx_type i = 0; i < n; i++)
                 mtmp(i,0) = workd[i + iptr(0) - 1];
 
-              mtmp = utsolve (bt, permB, m * ltsolve (b, permB, mtmp));
+              mtmp = ltsolve (b, permB, m * utsolve (bt, permB, mtmp));
 
               for (octave_idx_type i = 0; i < n; i++)
                 workd[i+iptr(1)-1] = mtmp(i,0);
@@ -1701,7 +1701,7 @@ EigsRealNonSymmetricMatrix (const M& m, 
             }
 
           if (note3)
-            eig_vec = ltsolve (M(b), permB, eig_vec);
+            eig_vec = utsolve (bt, permB, eig_vec);
         }
     }
   else
@@ -1870,7 +1870,7 @@ EigsRealNonSymmetricMatrixShift (const M
                 {
                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
 
-                  vector_product (m, workd+iptr(0)-1, dtmp);
+                  vector_product (b, workd+iptr(0)-1, dtmp);
 
                   Matrix tmp(n, 1);
 
@@ -2521,7 +2521,7 @@ EigsComplexNonSymmetricMatrix (const M& 
               ComplexMatrix mtmp (n,1);
               for (octave_idx_type i = 0; i < n; i++)
                 mtmp(i,0) = workd[i + iptr(0) - 1];
-              mtmp = utsolve (bt, permB, m * ltsolve (b, permB, mtmp));
+              mtmp = ltsolve (b, permB, m * utsolve (bt, permB, mtmp));
               for (octave_idx_type i = 0; i < n; i++)
                 workd[i+iptr(1)-1] = mtmp(i,0);
 
@@ -2610,7 +2610,7 @@ EigsComplexNonSymmetricMatrix (const M& 
             }
 
           if (note3)
-            eig_vec = ltsolve (b, permB, eig_vec);
+            eig_vec = utsolve (bt, permB, eig_vec);
         }
     }
   else
@@ -2785,7 +2785,7 @@ EigsComplexNonSymmetricMatrixShift (cons
                 {
                   OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
 
-                  vector_product (m, workd+iptr(0)-1, ctmp);
+                  vector_product (b, workd+iptr(0)-1, ctmp);
 
                   ComplexMatrix tmp(n, 1);
 
diff -r 3b2dbed26762 -r dcf30fdacb27 scripts/sparse/eigs.m
--- a/scripts/sparse/eigs.m	Thu Mar 16 09:44:44 2017 +0100
+++ b/scripts/sparse/eigs.m	Fri Mar 17 12:25:40 2017 +0100
@@ -760,6 +760,88 @@ endfunction
 %! v1 = v1(:, idx);
 %! assert (abs (v), abs (R \ v1), 1e-12);
 %!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [v, d] = eigs (A, B, 5, "lm");
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%! ddiag = diag (d);
+%! [ddiag, idx] = sort (ddiag);
+%! v = v(:, idx);
+%! R = chol (B);
+%! [v1, d1] = eigs (R' \ A / R, 5, "lm");
+%! d1diag = diag (d1);
+%! [d1diag, idx] = sort (d1diag);
+%! v1 = v1(:, idx);
+%! assert (abs (v), abs (R \ v1), 1e-12);
+%!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10) -...
+%! 1i * spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [v, d] = eigs (A, B, 5, "lm");
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%! ddiag = diag (d);
+%! [ddiag, idx] = sort (ddiag);
+%! v = v(:, idx);
+%! R = chol (B);
+%! [v1, d1] = eigs (R' \ A / R, 5, "lm");
+%! d1diag = diag (d1);
+%! [d1diag, idx] = sort (d1diag);
+%! v1 = v1(:, idx);
+%! assert (abs (v), abs (R \ v1), 1e-12);
+%!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = toeplitz (sparse (1:10));
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [v, d] = eigs (A, B, 5, 1);
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%! ddiag = diag (d);
+%! [ddiag, idx] = sort (ddiag);
+%! v = v(:, idx);
+%! R = chol (B);
+%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! d1diag = diag (d1);
+%! [d1diag, idx] = sort (d1diag);
+%! v1 = v1(:, idx);
+%! assert (abs (v), abs (R \ v1), 1e-12);
+%!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [v, d] = eigs (A, B, 5, 1);
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%! ddiag = diag (d);
+%! [ddiag, idx] = sort (ddiag);
+%! v = v(:, idx);
+%! R = chol (B);
+%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! d1diag = diag (d1);
+%! [d1diag, idx] = sort (d1diag);
+%! v1 = v1(:, idx);
+%! assert (abs (v), abs (R \ v1), 1e-12);
+%!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10) -...
+%! 1i * spdiags ([[1./(2:11)]',[-5:-2:-23]',[1:10]'],-1:1,10,10);
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [v, d] = eigs (A, B, 5, 1);
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%! ddiag = diag (d);
+%! [ddiag, idx] = sort (ddiag);
+%! v = v(:, idx);
+%! R = chol (B);
+%! [v1, d1] = eigs (R' \ A / R, 5, 1);
+%! d1diag = diag (d1);
+%! [d1diag, idx] = sort (d1diag);
+%! v1 = v1(:, idx);
+%! assert (abs (v), abs (R \ v1), 1e-12);
+%!testif HAVE_ARPACK, HAVE_UMFPACK
 %! A = toeplitz (sparse (1:10));
 %! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
 %! R = chol (B);
