Mercurial > octave

changeset 23297:f093b9d4ef30 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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).
author Marco Caliari <marco.caliari@univr.it>
date Fri, 17 Mar 2017 12:25:40 +0100
parents 3b2dbed26762
children 85ffe1bdd3a0
files liboctave/numeric/eigs-base.cc scripts/sparse/eigs.m
diffstat 2 files changed, 89 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- 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 @@
                 {
                   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 @@
               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 @@
             }
 
           if (note3)
-            eig_vec = ltsolve (M(b), permB, eig_vec);
+            eig_vec = utsolve (bt, permB, eig_vec);
         }
     }
   else
@@ -1870,7 +1870,7 @@
                 {
                   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 @@
               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 @@
             }
 
           if (note3)
-            eig_vec = ltsolve (b, permB, eig_vec);
+            eig_vec = utsolve (bt, permB, eig_vec);
         }
     }
   else
@@ -2785,7 +2785,7 @@
                 {
                   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);
 
--- 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 @@
 %! 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);