# HG changeset patch
# User Carnë Draug <carandraug@octave.org>
# Date 1435936713 -3600
# Node ID 0b9d23557506be8cc175fb29f46dd3dc0e391f42
# Parent  557979395ca9f75e3ccf2fae361368aaba6701eb
gallery: fix randsvd by adding missing dependency qmult().

* scripts/special-matrix/gallery.m (randsvd) was copied from the
  Test Matrix toolbox by Nicholas J. Higham.  It made use of qmult()
  which was also part of that toolbox but was left behind.  This qmult()
  implementation is also recovered from the Test Matrix toolbox.  Note
  that Octave itself used to have qmult() which is part of the legacy
  quaternion package (now also removed from the new quaternion package).
  See cset 21904fe299c8 for when qmult() was removed from Octave.
  Also fix the default value for KL and KU so that a 2 element vector
  can be used as N.

diff -r 557979395ca9 -r 0b9d23557506 scripts/special-matrix/gallery.m
--- a/scripts/special-matrix/gallery.m	Fri Jul 03 15:01:38 2015 +0100
+++ b/scripts/special-matrix/gallery.m	Fri Jul 03 16:18:33 2015 +0100
@@ -2150,7 +2150,8 @@
 
 endfunction
 
-function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = n-1, ku = kl)
+function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = max (n) -1,
+                      ku = kl)
   ## RANDSVD  Random matrix with pre-assigned singular values.
   ##   RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N
   ##   with COND(A) = KAPPA and singular values from the distribution MODE.
@@ -2851,6 +2852,53 @@
   endif
 endfunction
 
+## NOTE: qmult is part of the Test Matrix Toolbox and is used by randsvd()
+function B = qmult (A)
+  ## QMULT  Pre-multiply by random orthogonal matrix.
+  ##   QMULT(A) is Q*A where Q is a random real orthogonal matrix from
+  ##   the Haar distribution, of dimension the number of rows in A.
+  ##   Special case: if A is a scalar then QMULT(A) is the same as
+  ##   QMULT(EYE(A)).
+  ##
+  ##   Called by RANDSVD.
+  ##
+  ## Reference:
+  ##   G.W. Stewart, The efficient generation of random
+  ##   orthogonal matrices with an application to condition estimators,
+  ##   SIAM J. Numer. Anal., 17 (1980), 403-409.
+
+  [n, m] = size (A);
+
+  ##  Handle scalar A
+  if (isscalar (A))
+    n = A;
+    A = eye (n);
+  endif
+
+  d = zeros (n);
+
+  for k = n-1:-1:1
+    ## Generate random Householder transformation.
+    x = randn (n-k+1, 1);
+    s = norm (x);
+    sgn = sign (x(1)) + (x(1) == 0); # Modification for sign(1)=1.
+    s = sgn*s;
+    d(k) = -sgn;
+    x(1) = x(1) + s;
+    beta = s*x(1);
+
+    ## Apply the transformation to A.
+    y = x'*A(k:n,:);
+    A(k:n,:) = A(k:n,:) - x*(y/beta);
+  endfor
+
+  ## Tidy up signs
+  for i = 1:n-1
+    A(i,:) = d(i)*A(i,:);
+  endfor
+  A(n,:) = A(n,:) * sign (randn);
+  B = A;
+endfunction
 
 ## BIST testing for just a few functions to verify that the main gallery
 ## dispatch function works.