--- a/liboctave/numeric/eigs-base.cc	Fri Mar 10 13:15:23 2017 +0100
+++ b/liboctave/numeric/eigs-base.cc	Thu Mar 16 09:44:44 2017 +0100
@@ -78,12 +78,33 @@
 static M
 ltsolve (const SM& L, const ColumnVector& Q, const M& m)
 {
+  // Solve (Q_mat * L) * x = m, that is L * x = Q_mat' * m = m(Q)
   octave_idx_type n = L.cols ();
   octave_idx_type b_nc = m.cols ();
   octave_idx_type err = 0;
   double rcond;
   MatrixType ltyp (MatrixType::Lower);
-  M tmp = L.solve (ltyp, m, err, rcond, 0);
+  M retval (n, b_nc);
+  const double* qv = Q.fortran_vec ();
+  for (octave_idx_type j = 0; j < b_nc; j++)
+    {
+      for (octave_idx_type i = 0; i < n; i++)
+        retval.elem (i,j) = m.elem (static_cast<octave_idx_type>(qv[i]), j);
+    }
+  return L.solve (ltyp, retval, err, rcond, 0);
+}
+
+template <typename SM, typename M>
+static M
+utsolve (const SM& U, const ColumnVector& Q, const M& m)
+{
+  // Solve (U * Q_mat') * x = m by U * tmp = m, x(Q) = tmp (Q_mat * tmp = x)
+  octave_idx_type n = U.cols ();
+  octave_idx_type b_nc = m.cols ();
+  octave_idx_type err = 0;
+  double rcond;
+  MatrixType utyp (MatrixType::Upper);
+  M tmp = U.solve (utyp, m, err, rcond, 0);
   M retval;
   const double* qv = Q.fortran_vec ();
 
@@ -101,26 +122,6 @@
   return retval;
 }
 
-template <typename SM, typename M>
-static M
-utsolve (const SM& U, const ColumnVector& Q, const M& m)
-{
-  octave_idx_type n = U.cols ();
-  octave_idx_type b_nc = m.cols ();
-  octave_idx_type err = 0;
-  double rcond;
-  MatrixType utyp (MatrixType::Upper);
-
-  M retval (n, b_nc);
-  const double* qv = Q.fortran_vec ();
-  for (octave_idx_type j = 0; j < b_nc; j++)
-    {
-      for (octave_idx_type i = 0; i < n; i++)
-        retval.elem (i,j) = m.elem (static_cast<octave_idx_type>(qv[i]), j);
-    }
-  return U.solve (utyp, retval, err, rcond, 0);
-}
-
 static bool
 vector_product (const SparseMatrix& m, const double* x, double* y)
 {
@@ -199,7 +200,7 @@
     return false;
   else
     {
-      bt = fact.chol_matrix ();
+      bt = fact.chol_matrix (); // upper triangular
       b = bt.transpose ();
       permB = ColumnVector (n);
       for (octave_idx_type i = 0; i < n; i++)
@@ -218,7 +219,7 @@
     return false;
   else
     {
-      b = fact.L ();
+      b = fact.L (); // lower triangular
       bt = b.transpose ();
       permB = fact.perm () - 1.0;
       return true;
@@ -236,7 +237,7 @@
     return false;
   else
     {
-      bt = fact.chol_matrix ();
+      bt = fact.chol_matrix (); // upper triangular
       b = bt.hermitian ();
       permB = ColumnVector (n);
       for (octave_idx_type i = 0; i < n; i++)
@@ -256,7 +257,7 @@
     return false;
   else
     {
-      b = fact.L ();
+      b = fact.L (); // lower triangular
       bt = b.hermitian ();
       permB = fact.perm () - 1.0;
       return true;
@@ -766,7 +767,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);
@@ -854,7 +855,7 @@
             }
 
           if (note3)
-            eig_vec = ltsolve (b, permB, eig_vec);
+            eig_vec = utsolve (bt, permB, eig_vec);
         }
     }
   else

--- a/scripts/sparse/eigs.m	Fri Mar 10 13:15:23 2017 +0100
+++ b/scripts/sparse/eigs.m	Thu Mar 16 09:44:44 2017 +0100
@@ -143,8 +143,8 @@
 ##
 ## @item permB
 ## The permutation vector of the Cholesky@tie{}factorization of @var{B} if
-## @code{cholB} is true.  That is @code{chol (@var{B}(permB, permB))}.  The
-## default is @code{1:@var{n}}.
+## @code{cholB} is true.  It is obtained by @code{[R, ~, permB] =
+## chol (@var{B}, @qcode{"vector"})}. The default is @code{1:@var{n}}.
 ##
 ## @end table
 ##
@@ -743,6 +743,42 @@
 %! for i=1:k
 %!   assert (max (abs ((A - d1(i)*speye (n))*v1(:,i))), 0, 1e-11);
 %! endfor
+%!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, "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));
+%! R = chol (B);
+%! opts.cholB = R;
+%! [v, d] = eigs (A, R, 5, "lm", opts);
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+%!testif HAVE_ARPACK, HAVE_UMFPACK
+%! A = toeplitz (sparse (1:10));
+%! B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
+%! [R, ~, permB] = chol (B, "vector");
+%! opts.cholB = R;
+%! opts.permB = permB;
+%! [v, d] = eigs (A, R, 5, "lm", opts);
+%! for i = 1:5
+%!   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
+%! endfor
+
 
 #### FULL MATRIX VERSIONS ####