CD Concept 6/CD Concept 06.iso/mac/UTILITAIRE/RLaB/test

Text File | 1995-01-26 | 69.5 KB | 3,119 lines | [TEXT/ttxt]

//
// Test Suite for RLaB.
// This test suite is supposed to test as much as possible, and still
// be runable on most platforms. This means we cannot do graphics or
// other platform specific stuff like piping. However, that is OK,
// since we are mostly interested in assuring the builder that RLaB
// will produce correct numerical answers.

// We use randsvd, which in turn uses RLaB's random number
// generator. We try to be carefull and seed the random number
// generator so that each user will get similar results (or at least
// similar inputs to the tests).

// Test the other style comments
1 + 2; // A simple statement with trailing comment.
#  Optional RLaB comment style.
1 + 2; #  A simple statement with trailing comment.
%  Optional Matlab comment style.
1 + 2; %  A simple statement with trailing comment.

srand (SEED = 10);   // Seems to produce reasonable results
rand("default");

tic();    // Start timing the tests...

//
// Test Parameters and some functions we will need later
//

pi = 4.0*atan(1.0);
X = 3;    // should be 3 (heuristic).

//
// Compute machine epsilon
//

epsilon = function()
{
  eps = 1.0;
  while((eps + 1.0) != 1.0)
  {
    eps = eps/2.0;
  }
  return 2*eps;
};

eps = epsilon();

eye = function( m , n )
{
  if (!exist (n))
  {
    if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
    new = zeros (m[1], m[2]);
    N = min ([m[1], m[2]]);
  else
    if (class (m) == "string" || class (n) == "string") {
      error ("eye(), string arguments not allowed");
    }
    if (max (size (m)) == 1 && max (size (n)) == 1)
    {
      new = zeros (m[1], n[1]);
      N = min ([m[1], n[1]]);
    else
      error ("matrix arguments to eye() must be 1x1");
    }
  }
  for(i in 1:N)
  {
    new[i;i] = 1.0;
  }
  return new;
};

symm = function( A )
{
  return (A + A')./2;
};

//-------------------------------------------------------------------//

// Synopsis:    Pascal matrix.

// Syntax:    P = pascal ( N )

// Description:

//    The Pascal matrix of order N: a symmetric positive definite
//    matrix with integer entries taken from Pascal's triangle. The
//    Pascal matrix is totally positive and its inverse has integer
//    entries.  Its eigenvalues occur in reciprocal pairs. COND(P)
//    is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
//    lower triangular Cholesky factor (up to signs of columns) of
//    the Pascal matrix.   It is involutary (is its own
//    inverse). PASCAL(N,2) is a transposed and permuted version of
//    PASCAL(N,1) which is a cube root of the identity.

//      References:
//      R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
//           Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
//           gives a factorization of L = PASCAL(N,1) and a formula for the
//           elements of L^k).
//      S. Karlin, Total Positivity, Volume 1, Stanford University Press,
//           1968.  (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
//                   PASCAL(N) is a submatrix of this matrix.)
//      M. Newman and J. Todd, The evaluation of matrix inversion programs,
//           J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
//      H. Rutishauser, On test matrices, Programmation en Mathematiques
//           Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
//           1966, pp. 349-365.  (Gives an integral formula for the
//           elements of PASCAL(N).)
//      J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
//           Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
//      H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
//           Blackie, London and Glasgow, 1929.  (PASCAL(N,2) on page 332.)

//    This file is a translation of pascal.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

//-------------------------------------------------------------------//

pascal = function ( n , k )
{
  local(n, k)

  if (!exist (k)) { k = 0; }

  P = diag( (-1).^[0:n-1] );
  P[; 1] = ones(n,1);

  //  Generate the Pascal Cholesky factor (up to signs).

  for (j in 2:n-1)
  {
    for (i in j+1:n)
    {
      P[i;j] = P[i-1;j] - P[i-1;j-1];
    }
  }

  if (k == 0)
  {
    P = P*P';
  else if (k == 2) {
    P = rot90(P,3);
    if (n/2 == round(n/2)) { P = -P; }
  }}

  return P;
};

//-------------------------------------------------------------------//

// Synopsis:    Random, orthogonal upper Hessenberg matrix.

// Syntax:      H = ohess ( N )

// Description:

//      H is an N-by-N real, random, orthogonal upper Hessenberg
//      matrix. Alternatively, H = OHESS(X), where X is an arbitrary
//      real N-vector (N > 1) constructs H non-randomly using the
//      elements of X as parameters. In both cases H is constructed
//      via a product of N-1 Givens rotations.

//      Note: See Gragg (1986) for how to represent an N-by-N
//      (complex) unitary Hessenberg matrix with positive subdiagonal
//      elements in terms of 2N-1 real parameters (the Schur
//      parametrization). This M-file handles the real case only and
//      is intended simply as a convenient way to generate random or
//      non-random orthogonal Hessenberg matrices.

//      Reference:
//      W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
//      J. Comp. Appl. Math., 16 (1986), pp. 1-8.

//    This file is a translation of ohess.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

//-------------------------------------------------------------------//

ohess = function ( x )
{
  local (x)
  global (pi)

  if (any (imag (x))) { error("Parameter must be real."); }

  n = max(size(x));

  if (n == 1)
  {
    //  Handle scalar x.
    n = x;
    x = rand(n-1, 1)*2*pi;
    H = eye(n,n);
    H[n;n] = sign(rand());
  else
    H = eye(n,n);
    H[n;n] = sign(x[n]) + (x[n]==0);   // Second term ensures H[n;n] nonzero.
  }

  for (i in n:2:-1)
  {
    // Apply Givens rotation through angle x(i-1).
    theta = x[i-1];
    c = cos(theta);
    s = sin(theta);
    H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
                       -s*H[i-1;]+c*H[i;] ];
  }

  return H;
};

//-------------------------------------------------------------------//

// Synopsis:    Random matrix with pre-assigned singular values.

// Syntax:    R = randsvd (N, KAPPA, MODE, KL, KU)

// Description:

//    R is a (banded) random matrix of order N with COND(A) = KAPPA
//    and singular values from the distribution MODE.

//      N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
//      Available types:
//             MODE = 1:   one large singular value,
//             MODE = 2:   one small singular value,
//             MODE = 3:   geometrically distributed singular values,
//             MODE = 4:   arithmetically distributed singular values,
//             MODE = 5:   random singular values with unif. dist. logarithm.

//      If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
//      If MODE < 0 then the effect is as for ABS(MODE) except that in the
//      original matrix of singular values the order of the diagonal entries
//      is reversed: small to large instead of large to small.
//      KL and KU are the lower and upper bandwidths respectively; if they
//      are omitted a full matrix is produced.
//      If only KL is present, KU defaults to KL.
//      Special case: if KAPPA < 0 then a random full symmetric positive
//                    definite matrix is produced with COND(A) = -KAPPA and
//                    eigenvalues distributed according to MODE.
//                    KL and KU, if present, are ignored.

//      This routine is similar to the more comprehensive Fortran routine xLATMS
//      in the following reference:
//      J.W. Demmel and A. McKenney, A test matrix generation suite,
//      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
//      New York, 1989.

//    This file is a translation of randsvd.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

// Dependencies
#   rfile bandred
#   rfile qmult

//-------------------------------------------------------------------//

randsvd = function (n, kappa, mode, kl, ku)
{
  local (n, kappa, mode, kl, ku)

  if (!exist (kappa)) { kappa = sqrt(1/eps); }
  if (!exist (mode))  { mode = 3; }
  if (!exist (kl))    { kl = n-1; }    // Full matrix.
  if (!exist (ku))    { ku = kl; }     // Same upper and lower bandwidths.

  if (abs(kappa) < 1) {
    error("Condition number must be at least 1!");
  }

  posdef = 0;
  if (kappa < 0)            // Special case.
  {
    posdef = 1;
    kappa = -kappa;
  }

  p = min(n);
  m = n[1];        // Parameter n specifies dimension: m-by-n.
  n = n[max(size(n))];

  if (p == 1)        // Handle case where A is a vector.
  {
    rand("normal", -10, 10);
    A = rand(m, n);
    A = A/norm(A, "2");
    return A;
  }

  j = abs(mode);

  // Set up vector sigma of singular values.
  if (j == 3)
  {
    factor = kappa^(-1/(p-1));
    sigma = factor.^[0:p-1];

  else if (j == 4) {
    sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);

  else if (j == 5) {    // In this case cond(A) <= kappa.
    rand("uniform", 0, 1)
    sigma = exp( -rand(p,1)*log(kappa) );

  else if (j == 2) {
    sigma = ones(p,1);
    sigma[p] = 1/kappa;

  else if (j == 1) {
    sigma = ones(p,1)./kappa;
    sigma[1] = 1;

  }}}}}


  // Convert to diagonal matrix of singular values.
  if (mode < 0) {
    sigma = sigma[p:1:-1];
  }

  sigma = diag(sigma);

  if (posdef)        // Handle special case.
  {
    Q = qmult(p);
    A = Q'*sigma*Q;
    A = (A + A')/2;    // Ensure matrix is symmetric.
    return A;
  }

  if (m != n)
  {
    sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  }

  if (kl == 0 && ku == 0) // Diagonal matrix requested - nothing more to do.
  {
    A = sigma;
    return A;
  }

  // A = U*sigma*V, where U, V are random orthogonal matrices from the
  // Haar distribution.
  A = qmult(sigma');
  A = qmult(A');

  if (kl < n-1 || ku < n-1)    // Bandwidth reduction.
  {
   A = bandred(A, kl, ku);
  }

  rand("default");
  return A;
};

//-------------------------------------------------------------------//

// Synopsis:    Band reduction by two-sided unitary transformations.

// Syntax:    bandred ( A , KL, KU )

// Description:

//    bandred(A, KL, KU) is a matrix unitarily equivalent to A with
//    lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
//    i > j+KL or j > i+KU).  The reduction is performed using
//    Householder transformations. If KU is omitted it defaults to
//    KL.

//    Called by randsvd.
//    This is a `standard' reduction.  Cf. reduction to bidiagonal
//    form prior to computing the SVD.  This code is a little
//    wasteful in that it computes certain elements which are
//    immediately set to zero!
//
//      Reference:
//      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
//      Johns Hopkins University Press, Baltimore, Maryland, 1989.
//      Section 5.4.3.

//    This file is a translation of bandred.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

// Dependencies
#   rfile house

//-------------------------------------------------------------------//

bandred = function ( A , kl , ku )
{
  local (A, kl, ku)

  if (!exist (ku)) { ku = kl; else ku = ku; }

  if (kl == 0 && ku == 0) {
    error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  }

  // Check for special case where order of left/right transformations matters.
  // Easiest approach is to work on the transpose, flipping back at the end.

  flip = 0;
  if (ku == 0)
  {
    A = A';
    temp = kl; kl = ku; ku = temp; flip = 1;
  }

  m = A.nr; n = A.nc;

  for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  {
    if (j+kl+1 <= m)
    {
       ltmp = house(A[j+kl:m;j]);
       beta = ltmp.beta; v = ltmp.v;
       temp = A[j+kl:m;j:n];
       A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
       A[j+kl+1:m;j] = zeros(m-j-kl,1);
    }

    if (j+ku+1 <= n)
    {
       ltmp = house(A[j;j+ku:n]');
       beta = ltmp.beta; v = ltmp.v;
       temp = A[j:m;j+ku:n];
       A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
       A[j;j+ku+1:n] = zeros(1,n-j-ku);
    }
  }

  if (flip) {
    A = A';
  }

  return A;
};

//-------------------------------------------------------------------//

// Synopsis:    Lehmer matrix - symmetric positive definite.

// Syntax:      A = lehmer ( N )

// Description:

//      A is the symmetric positive definite N-by-N matrix with
//                     A(i,j) = i/j for j >= i.
//      A is totally nonnegative.  INV(A) is tridiagonal, and explicit
//      formulas are known for its entries.

//      N <= COND(A) <= 4*N*N.

//      References:
//        M. Newman and J. Todd, The evaluation of matrix inversion
//           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
//        Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
//           of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
//        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
//           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.

//    This file is a translation of lehmer.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

//-------------------------------------------------------------------//

lehmer = function ( n )
{
  local (n)
  global (tril)

  A = ones(n,1)*(1:n);
  A = A./A';
  A = tril(A) + tril(A,-1)';

  return A;
};

//-------------------------------------------------------------------//

// Synopsis:    Pre-multiply by random orthogonal matrix.

// Syntax:    B = qmult ( A )

// Description:

//    B 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.

//    This file is a translation of qmult.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

//-------------------------------------------------------------------//

qmult = function ( A )
{
  local (A)

  n = A.nr; m = A.nc;

  //  Handle scalar A.
  if (max(n,m) == 1)
  {
    n = A;
    A = eye(n,n);
  }

  d = zeros(n,n);

  for (k in n-1:1:-1)
  {
    // Generate random Householder transformation.
    rand("normal", 0, 1);
    x = rand(n-k+1,1);
    s = norm(x, "2");
    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);
  }

  // Tidy up signs.
  for (i in 1:n-1)
  {
    A[i;] = d[i]*A[i;];
  }

  A[n;] = A[n;]*sign(rand());
  B = A;

  rand("default");
  return B;
};

//-------------------------------------------------------------------//

// Synopsis:    Create a Householder matrix.

// Syntax:    house ( X )

//    If HOUSE(x), which returns a list containing elements `v' and
//    `beta',  then H = EYE - beta*v*v' is a Householder matrix such
//    that Hx = -sign(x(1))*norm(x)*e_1.
//    NB: If x = 0 then v = 0, beta = 1 is returned.
//            x can be real or complex.
//            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).

//    Theory: (textbook references Golub & Van Loan 1989, 38-43;
//             Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
//      Hx = y: (I - beta*v*v')x = -s*e_1.
//      Must have |s| = norm(x), v = x+s*e_1, and
//      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
//      So take s = sign(x(1))*norm(x) (which avoids cancellation).
//      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
//          = 2*norm(x)*(norm(x) + |x(1)|).

//    References:
//        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
//           Johns Hopkins University Press, Baltimore, Maryland, 1989.
//        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
//           New York, 1973,
//        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
//           Press, 1965.

//    This file is a translation of house.m from version 2.0 of
//    "The Test Matrix Toolbox for Matlab", described in Numerical
//    Analysis Report No. 237, December 1993, by N. J. Higham.

//-------------------------------------------------------------------//

house = function ( x )
{
  local (x)

  m = x.nr; n = x.nc;
  if (n > 1) { error ("Argument must be a column vector."); }

  s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  v = x;
  if (s == 0)    // Quit if x is the zero vector.
  {
    beta = 1;
    return << beta = beta ; v = v >>;
  }

  v[1] = v[1] + s;
  beta = 1/(s'*v[1]);        // NB the conjugated s.

  // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  // But beta as above can be non-real (due to rounding) only
  // when x is complex.

  return << beta = beta ; v = v >>;
};

//-------------------------------------------------------------------//

//  Syntax:    tril ( A )
//        tril ( A , K )

//  Description:

//  tril(x) returns the lower triangular part of A.

//  tril(A,K) returns the elements on and below the K-th diagonal of
//  A.

//  K = 0: main diagonal
//  K > 0: above the main diag.
//  K < 0: below the main diag.

//  See Also: triu
//-------------------------------------------------------------------//

tril = function(x, k)
{
  local(x, k)

  if (!exist (k)) { k = 0; }
  nr = x.nr; nc = x.nc;
  if(k > 0)
  {
    if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  else
    if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  }

  y = zeros(nr, nc);

  for(i in max( [1,1-k] ):nr) {
    j = 1:min( [nc, i+k] );
    y[i;j] = x[i;j];
  }

  return y;
};

//-------------------------------------------------------------------//

//  Syntax:    triu ( A )
//        triu ( A , K )

//  Description:

//  triu(x) returns the upper triangular part of A.

//  tril(x; k) returns the elements on and above the k-th diagonal of
//  A.

//  K = 0: main diagonal
//  K > 0: above the main diag.
//  K < 0: below the main diag.

//  See Also: tril
//-------------------------------------------------------------------//

triu = function(x, k)
{
  local(x, k)

  if (!exist (k)) { k = 0; }
  nr = x.nr; nc = x.nc;

  if(k > 0)
  {
    if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  else
    if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  }

  y = zeros(nr, nc);

  for(j in max( [1,1+k] ):nc) {
    i = 1:min( [nr, j-k] );
    y[i;j] = x[i;j];
  }

  return y;
};

//
//-------------------- Test Relational Expressions -------------------
//
printf("\tstart scalar tests...\n");
printf("\tstart relational tests...\n");

//    SCALAR CONSTANTS (REAL)
if( !(1<2) ) { error(); }
if( !(1<=2) ) { error(); }
if( 1>2 ) { error(); }
if( 1>=2 ) { error(); }
if( 1==2 ) { error(); }
if( !(1!=2) ) { error(); }
if( !1 ) { error(); }
if( !!!1) { error(); }

if( !([1]<[2]) ) { error(); }
if( !([1]<=[2]) ) { error(); }
if( [1]>[2] ) { error(); }
if( [1]>=[2] ) { error(); }
if( [1]==[2] ) { error(); }
if( !([1]!=[2]) ) { error(); }
if( ![1] ) { error(); }
if( !!![1]) { error(); }

//    SCALAR CONSTANTS (COMPLEX)
if( !(1+2i<2+3i) ) { error(); }
if( !(1+2i<=2+3i) ) { error(); }
if( 1+2i>2+3i ) { error(); }
if( 1+2i>=2+3i ) { error(); }
if( 1+2i==2+3i ) { error(); }
if( !(1+2i!=2+3i) ) { error(); }
if( !1+2i ) { error(); }
if( !!!1+2i) { error(); }

if( !([1+2i]<[2+3i]) ) { error(); }
if( !([1+2i]<=[2+3i]) ) { error(); }
if( [1+2i]>[2+3i] ) { error(); }
if( [1+2i]>=[2+3i] ) { error(); }
if( [1+2i]==[2+3i] ) { error(); }
if( !([1+2i]!=[2+3i]) ) { error(); }
if( ![1+2i] ) { error(); }
if( !!![1+2i]) { error(); }

//    SCALAR ENTITIES (REAL)
a=1;b=2;
if( !(a<b) ) { error(); }
if( !(a<=b) ) { error(); }
if( a>b ) { error(); }
if( a>=b ) { error(); }
if( a==b ) { error(); }
if( !(a!=b) ) { error(); }
if( !a ) { error(); }
if( !!!a) { error(); }

if( !([a]<[b]) ) { error(); }
if( !([a]<=[b]) ) { error(); }
if( [a]>[b] ) { error(); }
if( [a]>=[b] ) { error(); }
if( [a]==[b] ) { error(); }
if( !([a]!=[b]) ) { error(); }
if( ![a] ) { error(); }
if( !!![a]) { error(); }

//    SCALAR ENTITIES (COMPLEX)
a=1+2i;b=2+3i;
if( !(a<b) ) { error(); }
if( !(a<=b) ) { error(); }
if( a>b ) { error(); }
if( a>=b ) { error(); }
if( a==b ) { error(); }
if( !(a!=b) ) { error(); }
if( !a ) { error(); }
if( !!!a) { error(); }

if( !([a]<[b]) ) { error(); }
if( !([a]<=[b]) ) { error(); }
if( [a]>[b] ) { error(); }
if( [a]>=[b] ) { error(); }
if( [a]==[b] ) { error(); }
if( !([a]!=[b]) ) { error(); }
if( ![a] ) { error(); }
if( !!![a]) { error(); }

if (! all (all (-2 <  rand(4,4)))) { error (); }
if (! all (all (-2 <= rand(4,4)))) { error (); }
if (! all (all ( 2 >  rand(4,4)))) { error (); }
if (! all (all ( 2 >= rand(4,4)))) { error (); }

if (! all (all (rand(4,4) >  -2))) { error (); }
if (! all (all (rand(4,4) >= -2))) { error (); }
if (! all (all (rand(4,4) <   2))) { error (); }
if (! all (all (rand(4,4) <=  2))) { error (); }

if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }

//------------------------- Test REAL SCALARS ------------------------
//
//    CONSTANTS
//      Addition
if(1+2 != 3) { error(); }
//      Subtraction
if(1-2 != -1) { error(); }
//      Multiply
if(1*2 != 2) { error(); }
//      Divide
if(1/2 != 0.5) { error(); }
//      Power
if(2^2 != 4) { error(); }
if(4^0 != 1) { error(); }
//      Unary Minus
if(2-3 != -1) { error(); }
//
//    ENTITIES
//
a = 1; b = 2; c = 3; d = 0.5;
//      Addition
if(a+b != c) { error(); }
//      Subtraction
if(a-b != -a) { error(); }
//      Multiply
if(a*b != b) { error(); }
//      Divide
if(a/b != d) { error(); }
//      Power
if(b^b != b*b) { error(); }
if(b^0 != 1) { error (); }
//      Unary Minus
if(b-c != -a) { error(); }
//
//  ENTITIES & CONSTANTS
//
if(a+2 != c) { error(); }
//      Subtraction
if(1-b != -a) { error(); }
//      Multiply
if(a*2 != 2) { error(); }
//      Divide
if(1/b != d) { error(); }
//      Power
if(2^b != b*b) { error(); }
//      Unary Minus
if(b-3 != -a) { error(); }
//
//------------------------Test COMPLEX SCALARS -------------------------
//
//    CONSTANTS
if(sqrt(-1) != 1i) { error(); }
//      Addition
if((1+2i)+(2+3i) != (3+5i)) { error(); }
//      Subtraction
if((1+2i)-(3+4i) != (-2-2i)) { error(); }
//      Multiply
if((1+2i)*(3+4i) != (-5+10i)) { error(); }
//      Divide
if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
//      Power
//      Precision problems prevent us from testing these. Have to
//      be checked by hand.
//  (1+2i)^2 = -3 + 4i
//  (1+2i)^.5 = 1.272 + 7.862e-1i
//  if((1+2i)^2 != (-3+4i)) { error(); }
//  if((1+2i)^10 != (237+3116i)) { error(); }
//      Unary Minus
if(-(1+2i) != -1-2i) { error(); }
//
//    ENTITIES
//
a = 1+2i; b = 3+4i; c = 4+6i;
if(a+b != c) { error(); }
//      Subtraction
if(a-b != -2-2i) { error(); }
//      Multiply
if(a*b != -5+10i) { error(); }
//      Divide
//if(a/(3-4i) != -.2+.4i) { error(); }
//      Power
//  if(b^b != b*b) { error(); }
//      Unary Minus
if(-a != -1-2i) { error(); }
//
//    ENTITIES & CONSTANTS
//
if(a+(3+4i) != c) { error(); }
//      Subtraction
if((1+2i)-b != -2-2i) { error(); }
//      Multiply
if(a*(3+4i) != -5+10i) { error(); }
//      Divide
//if(a/(3-4i) != -.2+.4i) { error(); }
//      Power
//if(b^b != b*b) { error(); }
//      Unary Minus
if(-(1+2i) != -1-2i) { error(); }
//
// String - Numerical Equalities
//
if ((1 == "1")) { error(); }
if (([1] == "1")) { error(); }
if (("1" == 1)) { error(); }
if (("1" == [1])) { error(); }

if (rand(3,3) == "str") { error(); }
if ("str" == rand(3,3)) { error(); }
if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }

//
//----------------------- Test REAL MATRICES ---------------------------
//

printf("\tstart matrix tests...\n");
printf("\t\treal-matrices\n");

//  Read in test matrices
//

read("test.input");

//
//  Matrix construction
//

if(any([1;2;3] != [1,2,3]')) {
  error();
}

if(any (any (m0 != zeros(2,2)))) { error(); }
if(any (any (m1 != 1+zeros(2,2)))) { error(); }
if(any (any (m2 != [1,2;3,4]))) { error(); }
if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
if(any (any ([m2] != m2))) { error(); }

//
//  Matrix indexing
//

p = pascal(6);
if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }

//
//  Sub-Matrix promotion
//

if(m2[2;2] != 4) { error(); }
if(any(m2[2;] != [3,4])) { error(); }
if(any(m2[;2] != [2,4]')) { error(); }
i=2;j=1;
if(m2[i;j] != 3) { error(); }
i=1;j=2;
if(m2[i;j] != 2) { error(); }
m = [1,2,3;4,5,6;7,8,9];

if(any(m[1;1,2] != [1,2]))
{
  error();
}

if(any(m[1,2;1] != [1;4]))
{
  error();
}

if(any (any (m[1,2;1,2] != [1,2;4,5])))
{
  error();
}

//

if(m3[2;2] != (5+6i)) { error(); }
if(any(m3[2;] != [3+4i,5+6i])) { error(); }
if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }

//
//  Automatic creation, extension
//

if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
a=[1,2,3;4,5,6];
if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
a=[1,2;3,4];
if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8])))
{
  error();
}

mmm[2;] = a[1;];

//
//  Matrix binary operations
//

a = m2; b = [5,6;7,8];
if(any (any (a+a != [2,4;6,8]))) { error(); }
if(any (any (a-a != zeros(2,2)))) { error(); }
if(any (any (2+a != [3,4;5,6]))) { error(); }
if(any (any (2-a != [1,0;-1,-2]))) { error(); }
if(any (any (a-2 != [-1,0;1,2]))) { error(); }
if(any (any (2*a != [2,4;6,8]))) { error(); }
if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
if(any (any (a*a != [7,10;15,22]))) { error(); }
if(any (any (a*a*a != [37,54;81,118]))) { error(); }
if(any (any (a .* a != [1,4;9,16]))) { error(); }
if(any (any (a./a != [1,1;1,1]))) { error(); }
if(any (any (a' != [1,3;2,4]))) { error(); }

if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }

if(any ([1;2;3]+[4;5;6] != [5;7;9]))
{
  error();
}

if(any ([1;2;3]-[4;5;6] != [-3;-3;-3]))
{
  error();
}

if(any ([2;2;2] ./ [1;1;1] != [2;2;2]))
{
  error();
}

if(any ([1;2;3] .* [4;5;6] != [4;10;18]))
{
  error();
}

if (type (1^.33333) != "real") { error (); }
if (type (1^(1/3)) != "real") { error (); }
if (type ([1]^.33333) != "real") { error (); }
if (type (1^[.33333]) != "real") { error (); }
if (type ([1]^[.33333]) != "real") { error (); }

//
// Test row-wise matrix addition
//

a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }

a = [1,2,3] + [1,2,3]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
if (!all (all (a .+ b == c))) { error (); }
if (!all (all (b .+ a == c))) { error (); }

printf("\t\tpassed matrix row-wise add test...\n");

//
// Test row-wise matrix subtraction
//

a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }

a = [1,1,1] + [1,1,1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
if (!all (all (a .- b == -c))) { error (); }
if (!all (all (b .- a ==  c))) { error (); }

printf("\t\tpassed matrix row-wise subtraction test...\n");

//
// Test col-wise matrix addition
//

a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }

a = [1;1;1;1] + [1;1;1;1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
if (!all (all (a .+ b == c))) { error (); }
if (!all (all (b .+ a == c))) { error (); }

printf("\t\tpassed matrix col-wise add test...\n");

//
// Test col-wise matrix subtraction
//

a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }

a = [1;1;1;1] + [1;1;1;1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
if (!all (all (a .- b == -c))) { error (); }
if (!all (all (b .- a ==  c))) { error (); }

printf("\t\tpassed matrix col-wise subtraction test...\n");

a = [1,2,3];
b = [1,2,3;4,5,6;7,8,9];
c = [1,4,9;4,10,18;7,16,27];

if (!all (all (a .* b == c))) { error (); }
if (!all (all (b .* a == c))) { error (); }

za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;

if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
if (!all (all (b .* za == b .* [za;za;za]))) { error (); }

printf("\t\tpassed matrix row-wise multiplication test...\n");

a = [1,2,3];
b = [1,2,3;4,6,6;7,8,9];
c = [1,1,1;4,3,2;7,4,3];

if (!all (all (b ./ a == c))) { error (); }
if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }

za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;

if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }

printf("\t\tpassed matrix row-wise division test...\n");

a = [1;2;3];
b = [1,2,3;4,5,6;7,8,9];

if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
if (!all (all (b .* a == b .* [a,a,a]))) { error (); }

za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;

if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }

printf("\t\tpassed matrix column-wise multiplication test...\n");

a = [1;2;3];
b = [1,2,3;4,6,6;7,8,9];

if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }

za = a + rand (size(a))*1j;
zb = b + rand (size(b))*1j;

if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }

printf("\t\tpassed matrix column-wise division test...\n");


//
//--------------------- Test COMPLEX MATRICES --------------------------
//
//  Automatic creation, extension
//
printf("\t\tcomplex-matrices\n");
if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
a=[1,2,3;4,5,6];
if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
//
a = m3;
if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
if(any (any (a-a != zeros(2,2)))) { error(); }
if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
//
// The following test may not work on some machines
//
if(any (any (a./a != [1,1;1,1]))) {
  printf("\t\t***complex division inaccuracy, check manually***\n");
}

if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
//
//--------------------- Test NULL MATRICES -------------------------
//
printf("\t\tnull-matrices\n");
// Create a NULL matrix
mnull = [];
// Test it with SCALARS
if( any([1,mnull] != 1)) {
  error();
}
if( any([mnull,1] != 1)) {
  error();
}
// Test with MATRIX construction
m = [1,2;3,4;5,6];
if( any([mnull;1] != [1])) {
  error();
}
if( any([1;mnull] != [1])) {
  error();
}
if( any([mnull;1,2,3] != [1,2,3])) {
  error();
}
if( any([1,2,3;mnull] != [1,2,3])) {
  error();
}
if(any (any ([mnull,m] != m))) {
  error();
}
if(any (any ([m,mnull] != m))) {
  error();
}
if(any (any ([mnull;m] != m))) {
  error();
}
if(any (any ([m;mnull] != m))) {
  error();

mnull = matrix();
mnull[1] = [1];
}

//--------------------- Test Matrix Multiply --------------------------

i = sqrt(-1);
a = [1,2,3;4,5,6;7,8,9];
b = [4,5,6;7,8,9;10,11,12];
c = [ 48,  54,  60 ;
     111, 126, 141 ;
     174, 198, 222 ];

if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }

az = a + b*i;
bz = b + a*i;

cz = [-18+141*i , -27+162*i , -36+183*i ;
        9+240*i ,   0+279*i ,  -9+318*i ;
       36+339*i ,  27+396*i ,  18+453*i ];

czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
       111+66*i , 126+81*i  , 141+96*i ;
       174+102*i, 198+126*i , 222+150*i ];

czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
        111+174*i , 126+198*i , 141+222*i ;
        174+237*i , 198+270*i , 222+303*i ];

if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }

a = [a,a];
b = [b;b];
c = [  96 , 108 , 120 ;
      222 , 252 , 282 ;
      348 , 396 , 444 ];

if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }

az = [az,az];
bz = [bz;bz];

cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
         18+480*i ,    0+558*i ,  -18+636*i ;
         72+678*i ,   54+792*i ,   36+906*i ];

czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
        222+132*i , 252+162*i , 282+192*i ;
        348+204*i , 396+252*i , 444+300*i ];

czzz = [  96+222*i , 108+252*i , 120+282*i ;
         222+348*i , 252+396*i , 282+444*i ;
         348+474*i , 396+540*i , 444+606*i ];

if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }

printf("\t\tpassed matrix multiply test...\n");

//--------------------- Test STRING MATRICES --------------------------
//
printf("\t\tstring-matrices\n");
sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
if(sm[1] != "s1") { error(); }
if( sm[1;3] != "sm3" ) { error(); }
if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
{
  error();
}
if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }

//  Test string-matrix add functionality

sm = sm[1,2;1,2];
if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}

if ("c"+["1"] != "c1") { error (); }
if (["c"]+"1" != "c1") { error (); }

printf("\tpassed matrix test...\n");

//
//---------------------------- List Tests --------------------------
//
//  List creation
listest = << << 11; 12 >>; << 21; 22>> >>;
if( listest.[1].[2] != 12 ) { error(); }
if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
mlist.[0] = m;
if(any(any(mlist.[0] != m))) { error(); }

// Test list functions...

listtest.fun = function ( a ) { return sin(a); };
if (!(listtest.fun(0.5) == sin(0.5))) { error() ; }

// Test open-list assignment...
</v;d/> = eig(rand(3,3));

printf("\tpassed list test...\n");

//
// Reset random number generator seed...
//

rand("default");

//
//-------------------------- Test printf () --------------------------
//

sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
printf("\tpassed sprintf test...\n");

//
//------------------------- Test strtod()  ----------------------------
//

if (123.456 != strtod ("123.456")) { error (); }
if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  { error (); }
printf("\tpassed strtod test...\n");

//
//------------------------- Test getline()  ---------------------------
//
//

close( "test.getline" );

x = getline( "test.getline" );
if ( x.[1] !=  123.456 ) { error(); }
if ( x.[2] != -123.456 ) { error(); }
if ( x.[3] !=  123.456 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  .123 ) { error(); }
if ( x.[2] != -.123 ) { error(); }
if ( x.[3] !=  .123 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  123 ) { error(); }
if ( x.[2] != -123 ) { error(); }
if ( x.[3] !=  123 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  1e6 ) { error(); }
if ( x.[2] != -1e6 ) { error(); }
if ( x.[3] !=  1e6 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  1e5 ) { error(); }
if ( x.[2] != -1e5 ) { error(); }
if ( x.[3] !=  1e5 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  123.456e3 ) { error(); }
if ( x.[2] != -123.456e3 ) { error(); }
if ( x.[3] !=  123.456e3 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  123.456e3 ) { error(); }
if ( x.[2] != -123.456e3 ) { error(); }
if ( x.[3] !=  123.456e3 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  123.456e-3 ) { error(); }
if ( x.[2] != -123.456e-3 ) { error(); }
if ( x.[3] !=  123.456e-3 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  .123e3 ) { error(); }
if ( x.[2] != -.123e3 ) { error(); }
if ( x.[3] !=  .123e3 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] !=  123e3 ) { error(); }
if ( x.[2] != -123e3 ) { error(); }
if ( x.[3] !=  123e3 ) { error(); }

x = getline( "test.getline" );
if ( x.[1] != "123abc" ) { error(); }
if ( x.[2] != "abc123" ) { error(); }
if ( x.[3] != "123.abc" ) { error(); }

x = getline( "test.getline" );
if ( x.[1] != "quoted string" ) { error(); }
if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }

x = getline( "test.getline" );
if ( x.[1] != "quoted string" ) { error(); }
if ( x.[2] !=  1.23e3 ) { error(); }
if ( x.[3] !=  100 ) { error(); }
if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
if ( x.[5] !=  200 ) { error(); }

printf("\tpassed getline() test...\n");

//
//---------------------- Test readb()/writeb()  --------------------
//

a = rand (5,5);
z = rand(3,3) + rand(3,3)*1j;
strm = what()[1:5;1:5];
pi = 4*atan(1.0);
sc = 2*pi;
str = "this is a sample string\ttab";
l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;

writeb ("jnk.rb", a, z, strm, sc, str, l);

# Set aside matrices for later tests

check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;

clear (a, z, strm, sc, str, l);

close ("jnk.rb");

readb ("jnk.rb");

#
# Now do checks
#

if (!all (all (a == check.a))) { error (); }
if (!all (all (z == check.z))) { error (); }
if (!all (all (strm == check.strm))) { error (); }
if (sc != check.sc) { error (); }
if (str != check.str) { error (); }

if (length (l) != 5) { error (); }
if (!all (all (l.a == check.a))) { error (); }
if (!all (all (l.z == check.z))) { error (); }
if (!all (all (l.strm == check.strm))) { error (); }
if (l.sc != check.sc) { error (); }
if (l.str != check.str) { error (); }


printf("\tpassed binary I/O test...\n");

//
//---------------------- Test read()/write()  --------------------
//

a = rand (5,5);
z = rand(3,3) + rand(3,3)*1j;
strm = what()[1:5;1:5];
pi = 4*atan(1.0);
sc = 2*pi;
str = "this is a sample string\ttab";
l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;

write ("jnk.ra", a, z, strm, sc, str, l);

# Set aside matrices for later tests

check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;

clear (a, z, strm, sc, str, l);

close ("jnk.ra");

read ("jnk.ra");

#
# Now do checks
#

if (a.nr != 5 || a.nc != 5) { error (); }
if (z.nr != 3 || z.nc != 3) { error (); }
if (strm.nr != 5 || strm.nc != 5) { error (); }
if (str != check.str) { error (); }

if (length (l) != 5) { error (); }
if (l.a.nr != 5 || l.a.nc != 5) { error (); }
if (l.z.nr != 3 || l.z.nc != 3) { error (); }
if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
if (l.str != check.str) { error (); }

printf("\tpassed ascii I/O test...\n");

//
//-------------------------- Test eval () ------------------------------
//

if (1 + 2 != eval("1 + 2")) { error ("eval() error"); }
x = function (s, a) { return eval(s); };
str = "yy = 2 + x(\"2*a\",3)";
z = eval(str);
if (z != yy) { error ("eval() error"); }
printf("\tpassed eval test...\n");


//-------------------------- ------------------------------

printf ("\ttest builtin function correctness...\n");

//-------------------------- ------------------------------

//
//-------------------------- Test abs () ------------------------------
//

A = rand(5,5);
T = ( A == abs (A));
if (!all (all (A))) { error ("abs() incorrect"); }
printf("\tabs()");


//
//-------------------------- Test max () ------------------------------
//

A = [1,10,100;2,20,200;1,2,3];
B = A./2;
ZA = A + rand (3,3)*A*1i;
ZB = B + rand (3,3)*B*1i;
if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
if (max (max(A)) != 200) { error ("max() incorrect"); }
if (any (any (max (A, B) != A))) { error (); }
if (any (any (max (B, A) != max (A, B)))) { error (); }
if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
printf("\tmax()");

//
//-------------------------- Test min () ------------------------------
//

if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
if (min (min(A)) != 1) { error ("min() incorrect"); }
if (any (any (min (A, B) != B))) { error (); }
if (any (any (min (B, A) != min (A, B)))) { error (); }
if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
printf("\tmin()");

//
//-------------------------- Test maxi () -----------------------------
//

if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
printf("\tmaxi()");

//
//-------------------------- Test mini () -----------------------------
//

if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
printf("\tmini()");

//
//-------------------------- Test floor () ----------------------------
//

if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  error ("floor output incorrect");
}
printf("\tfloor()");

//
//-------------------------- Test ceil () 0----------------------------
//

if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  error ("ceil output incorrect");
}
printf("\tceil()\n");

//
//-------------------------- Test round () ----------------------------
//

if (round (1.8) != 2) { error ("round() output incorrect"); }
if (round (1.4) != 1) { error ("round() output incorrect"); }
if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  error ("round output incorrect");
}
printf("\tround()");

//
//-------------------------- Test sum () ------------------------------
//

S = [1:4; 4:7; 8:11];
if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
printf("\tsum()");

//
//-------------------------- Test int () ------------------------------
//

if (int (1.9999) != 1) { error ("int() output incorrect"); }
if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  error ("int() output incorrect");
}
printf("\tint()");

//
//-------------------------- Test mod () ------------------------------
//

if (mod (1,1) != 0) { error ("mod() output incorrect"); }
if (mod (4,2) != 0) { error ("mod() output incorrect"); }
if (mod (3,2) != 1) { error ("mod() output incorrect"); }
if (mod (5,3) != 2) { error ("mod() output incorrect"); }
printf("\tmod()");

//
//-------------------------- Test find () ------------------------------
//

if (find ([0,1]) != 2) { error ("find() output incorrect"); }
if (find ([1,0]) != 1) { error ("find() output incorrect"); }
if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }

printf("\tfind()\n");


//-------------------------- ------------------------------

printf ("\ttest builtin function operation...\n");

//-------------------------- ------------------------------

//
// At least use most of the builtins....
//

A = rand(5,5);
Z = rand(5,5) + rand(5,5)*1j;
S = what()[3:6;2:4];

abs (A);
abs ([A]);
abs (3.14);

abs (Z);
abs ([Z]);
abs (3.14j);

printf ("\tabs");

acos (A);
acos ([A]);
acos (3.14);

acos (Z);
acos ([Z]);
acos (3.14j);

printf ("\tacos");

all (A);
all ([A]);
all (3.14);

all (Z);
all ([Z]);
all (3.14j);

printf ("\tall");

any (A);
any ([A]);
any (3.14);

any (Z);
any ([Z]);
any (3.14j);

printf ("\tany");

asin (A);
asin ([A]);
asin (3.14);

asin (Z);
asin ([Z]);
asin (3.14j);

printf ("\tasin");

atan (A);
atan ([A]);
atan (3.14);

atan (Z);
atan ([Z]);
atan (3.14j);

printf ("\tatan");

atan2 (A,A);
atan2 ([A],[A]);
atan2 (3.14,3.14);

#atan2 (Z,Z);
#atan2 ([Z],[Z]);
#atan2 (3.14j,3.14j);

printf ("\tatan2");
printf ("\n");

cd (":");
printf ("\tcd");

ceil (A);
ceil ([A]);
ceil (3.14);

ceil (Z);
ceil ([Z]);
ceil (3.14j);

printf ("\tceil");

class (A);
class ([A]);
class (3.14);
class ("string");
class (S);
class ([S]);
class (<< rand(3,3); rand(4,4) >>);

class (Z);
class ([Z]);
class (3.14j);

printf ("\tclass");

clear (rand(3,3));
clear ([rand(3,3)]);

printf ("\tclear");

conj (A);
conj ([A]);
conj (3.14);

conj (Z);
conj ([Z]);
conj (3.14j);

printf ("\tconj");

cos (A);
cos ([A]);
cos (3.14);

cos (Z);
cos ([Z]);
cos (3.14j);

printf ("\tcos");

diag (A);
diag ([A]);
diag (3.14);

diag (Z);
diag ([Z]);
diag (3.14j);

printf ("\tdiag");
printf ("\n");

exist (A);
exist (Z);

printf ("\texist");

exp (A);
exp ([A]);
exp (3.14);

exp (Z);
exp ([Z]);
exp (3.14j);

printf ("\texp");

find (A);
find ([A]);
find (3.14);

find (Z);
find ([Z]);
find (3.14j);

printf ("\tfind");

floor (A);
floor ([A]);
floor (3.14);

floor (Z);
floor ([Z]);
floor (3.14j);

printf ("\tfloor");

format (1,1);
format ([2],[3]);
format ();
format ([3], 3);
format ();

printf ("\tformat");

imag (A);
imag ([A]);
imag (3.14);

imag (Z);
imag ([Z]);
imag (3.14j);

printf ("\timag");

inf ();
printf ("\tinf");
printf ("\n");

int (A);
int ([A]);
int (3.14);

int (Z);
int ([Z]);
int (3.14j);

printf ("\tint");

issymm (A);
issymm ([A]);
issymm (3.14);

issymm (Z);
issymm ([Z]);
issymm (3.14j);

printf ("\tissymm");

length (A);
length ([A]);
length (3.14);

length (Z);
length ([Z]);
length (3.14j);

length (3.14);
length (S);
length (<< rand(2,2); rand(10,10); 3.14 >>);

printf ("\tlength");

log (A);
log ([A]);
log (3.14);

log (Z);
log ([Z]);
log (3.14j);

printf ("\tlog");

log10 (A);
log10 ([A]);
log10 (3.14);

log10 (Z);
log10 ([Z]);
log10 (3.14j);

printf ("\tlog10");

max (A);
max ([A]);
max (3.14);

max (Z);
max ([Z]);
max (3.14j);

printf ("\tmax");

maxi (A);
maxi ([A]);
maxi (3.14);

maxi (Z);
maxi ([Z]);
maxi (3.14j);

printf ("\tmaxi");
printf ("\n");

members (<< rand(2,2); rand(); rand(3,3) >>);

printf ("\tmembers");

min (A);
min ([A]);
min (3.14);

min (Z);
min ([Z]);
min (3.14j);

printf ("\tmin");

mini (A);
mini ([A]);
mini (3.14);

mini (Z);
mini ([Z]);
mini (3.14j);

printf ("\tmini");

mod (A,A);
mod ([A],A);
mod (3.14,2);

mod (Z,Z);
mod ([Z],Z);
mod (3.14j,2);

printf ("\tmod");

nan ();

ones(3,3);
ones([4,4]);

printf ("\tones");

prod (A);
prod ([A]);
prod (3.14);

prod (Z);
prod ([Z]);
prod (3.14j);

printf ("\tprod");

real (A);
real ([A]);
real (3.14);

real (Z);
real ([Z]);
real (3.14j);

printf ("\treal");
printf ("\n");

reshape (A, A.nr*A.nc, 1);
reshape (A, [A.nr*A.nc], [1]);
reshape (Z, Z.nr*Z.nc, 1);
reshape (Z, [Z.nr*Z.nc], [1]);

printf ("\treshape");

round (A);
round ([A]);
round (3.14);

round (Z);
round ([Z]);
round (3.14j);

printf ("\tround");

show(A);
show(Z);
show([A]);
show([Z]);

printf ("\tshow");

sign (A);
sign ([A]);
sign (3.14);

sign (Z);
sign ([Z]);
sign (3.14j);

printf ("\tsign");

sin (A);
sin ([A]);
sin (3.14);

sin (Z);
sin ([Z]);
sin (3.14j);

printf ("\tsin");

size (A);
size ([A]);
size (3.14);

size (Z);
size ([Z]);
size (3.14j);

printf ("\tsize");

sizeof (A);
sizeof ([A]);
sizeof (3.14);

sizeof (Z);
sizeof ([Z]);
sizeof (3.14j);

printf ("\tsizeof");
printf ("\n");

sort (A);
sort ([A]);
sort (3.14);

sort (Z);
sort ([Z]);
sort (3.14j);

printf ("\tsort");

sqrt (A);
sqrt ([A]);
sqrt (3.14);

sqrt (Z);
sqrt ([Z]);
sqrt (3.14j);

printf ("\tsqrt");

srand ();
srand ("clock");
srand (SEED);

printf ("\tsrand");

strsplt (S);
printf ("\tstrsplt");

strtod ("string");
strtod (S);
strtod ("1.2");
strtod (["1.2", "1e3"]);

printf ("\tstrtod");

sum (A);
sum ([A]);
sum (3.14);

sum (Z);
sum ([Z]);
sum (3.14j);

printf ("\tsum");

tan (A);
tan ([A]);
tan (3.14);

tan (Z);
tan ([Z]);
tan (3.14j);

printf ("\ttan");
printf ("\n");

tmpnam ();
printf ("\ttmpnam");

type (A);
type ([A]);
type (Z);
type ([Z]);
type (3.14);
type (S);
type ("string");

printf ("\ttype");

zeros (3,3);
zeros ([3], [3]);
zeros ([3,3]);

printf ("\tzeros");

printf ("\n");

srand(SEED);
rand("default");

//
//-------------------------- Test rand () -----------------------------
//
rand ("normal", 5, 1);
xrand = rand(4000,1);

mean = function(x)
{
  local(m);

  m = size (x)[1];
  if( m == 1 )
  {
    m = size (x)[2];
  }

  return sum( x ) / m;
};

std = function(x)
{
  local(i, m, s);

  if(class(x) != "num") { error("std() requires NUMERICAL input"); }

  m = x.nr;
  if( m == 1 )
  {
    return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  else
    for( i in 1:x.nc) {
      s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
    }
    return s;
  }
};

if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1))
  { error ("error in random"); }
if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  { error ("error in random"); }
printf("\tpassed rand test...\n");

rand("default");

//
//-------------------------- Test norm () -----------------------------
//

tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
printf("\tpassed norm test...\n");

//
//-------------------------- Test qr () ------------------------------
//

a = ohess(4);
qa = qr (a);
if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  { error ("possible qr() problems"); }

z = ohess (4) + ohess(4)*1i;
qz = qr (z);
if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  { error ("possible qr() problems"); }

printf("\tpassed qr test...\n");

//
//-------------------------- Test schur () ----------------------------
//

a = randsvd (10, 10);
sa = schur (a);
if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  { error ("possible schur() problems"); }

z = rand (4,4) + rand(4,4)*1i;
sz = schur (z);
if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  { error ("possible schur() problems"); }

a = randsvd (10, -10);
sa = schur (a);
if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  { error ("possible schur() problems"); }

z = rand (4,4) + rand(4,4)*1i;
sz = schur (z);
if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  { error ("possible schur() problems"); }

printf("\tpassed schur test...\n");

//
//-------------------------- Test schord () ----------------------------
//

s2a = schord (sa, 2, 4);
if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
{
  error ("possible schord() problems");
}

s2z = schord (sz, 3, 1);
if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
{
  error ("possible schord() problems");
}

printf("\tpassed schord test...\n");

//
//-------------------------- Test chol () -----------------------------
//

c = lehmer(10);
u = chol (c);
if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
{
  error ("possible chol() problems");
}

cz = lehmer(10) + lehmer(10)*1j;
cz = symm(cz);
uz = chol (cz);
if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
{
  error ("possible chol() problems");
}

printf("\tpassed chol test...\n");

//
//--------------------------- Test inv () --------------------------------
//

a = randsvd(10,10);
b = ones(10,1);
x = inv(a) * b;
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
  printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  printf ("\tA*X - B:\n");
  abs (a*x - b)
  error ("possible inv() problems\n");
}

az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
xz = inv (az)*bz;
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
  printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  printf ("\tA*X - B:\n");
  abs (az*xz - bz)
  error ("possible inv() problems\n");
}

printf("\tpassed inv test...\n");

//
//-------------------------- Test solve () -----------------------------
//

//
// Real - General case
//

a = randsvd(10,10);
b = ones(10,1);
x = solve (a,b);
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
  printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  printf ("\tA*X - B:\n");
  abs (a*x - b)
  error ("possible solve() problems\n");
}

//
// Real - Symmetric case
//

s = symm (randsvd(10,10));
b = ones(10,1);
x = solve (s,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
  printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  printf ("\tA*X - B:\n");
  abs (s*x - b)
  error ("possible solve() problems\n");
}

//
// Complex - General  case
//

az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
xz = solve (az,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
  printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  printf ("\tA*X - B:\n");
  abs (az*xz - bz)
  error ("possible solve() problems\n");
}

//
// Complex - Symmetric case
//

sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
bz = rand(10,1) + rand(10,1)*1j;
xz = solve (sz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
  printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  printf ("\tA*X - B:\n");
  abs (sz*xz - bz)
  error ("possible solve() problems\n");
}

printf("\tpassed solve test...\n");

//
//-------------------------- Test factor() / backsub() -----------------------------
//

//
// Real - General case
//

a = randsvd(10,10);
b = ones(10,1);
f = factor (a);
x = backsub (f,b);
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
  printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  printf ("\tA*X - B:\n");
  abs (a*x - b)
  error ("possible factor/backsub problems\n");
}

//
// Real - Symmetric case
//

s = symm (randsvd(10,10));
f = factor (s);
x = backsub (f,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
  printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  printf ("\tA*X - B:\n");
  abs (s*x - b)
  error ("possible factor/backsub problems\n");
}

s = symm (randsvd(10,10));
f = factor (s, "s");
x = backsub (f,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
  printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  printf ("\tA*X - B:\n");
  abs (s*x - b)
  error ("possible factor/backsub problems\n");
}

//
// Complex - General  case
//

az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
fz = factor (az);
xz = backsub (fz,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
  printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  printf ("\tA*X - B:\n");
  abs (az*xz - bz)
  error ("possible factor/backsub problems\n");
}

az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
fz = factor (az, "g");
xz = backsub (fz,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
  printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  printf ("\tA*X - B:\n");
  abs (az*xz - bz)
  error ("possible factor/backsub problems\n");
}

//
// Complex - Symmetric case
//

sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
fz = factor(sz);
xz = backsub (fz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
  printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  printf ("\tA*X - B:\n");
  abs (sz*xz - bz)
  error ("possible factor/backsub problems\n");
}

sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
fz = factor(sz, "s");
xz = backsub (fz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
  printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  printf ("\tA*X - B:\n");
  abs (sz*xz - bz)
  error ("possible factor/backsub problems\n");
}

printf("\tpassed factor/backsub test...\n");

//
//------------------------------ Test lu() ---------------------------------
//

static (swap);

lu = function ( A )
{
  local (i, l, u, pvt, x)

  if (A.nr != A.nc) { error ("lu() requires square A"); }

  x = factor (A, "g");    // Do the factorization

  //
  // Now create l, u, and pvt from lu and pvt.
  //

  l = tril (x.lu, -1) + eye (size (x.lu));
  u = triu (x.lu);
  pvt = eye (size (x.lu));

  //
  // Now re-arange the columns of pvt
  //

  for (i in 1:max (size (x.lu)))
  {
    pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  }
  return << l = l; u = u; pvt = pvt >>;
};

//
//  In vector V, swap elements I, J
//

swap = function ( V, I, J )
{
  local (v, tmp);
  v = V;
  tmp = v[I];
  v[I] = v[J];
  v[J] = tmp;
  return v;
};

a = randsvd(10,10);
lua = lu (a);
if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
{
  printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  printf ("\tA - p*l*u:\n");
  abs (a - lua.pvt*lua.l*lua.u)
  error ("possible lu()/factor() problems\n");
}

//
// Real
az = randsvd(10,10) + randsvd(10,10)*1j;
luz = lu (az);
if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
{
  printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  printf ("\tA - p*l*u:\n");
  abs (az - luz.ovt*luz.l*luz.u)
  error ("possible lu()/factor()() problems\n");
}

printf("\tpassed lu/factor test...\n");

//
//-------------------------- Test svd ()   -----------------------------
//

a = randsvd(10,10);
s = svd (a);
if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
{
  error ("possible svd() problems");
}

z = randsvd(10,10) + rand(10,10)*1j;
sz = svd (z);
if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
{
  error ("possible svd() problems");
}

printf("\tpassed svd test...\n");

//
//-------------------------- Test hess ()  -----------------------------
//

a = randsvd(10,10);
h = hess (a);
if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
{
  error ("possible hess() problems");
}

z = randsvd(10,10) + randsvd(10,10)*1j;
hz = hess (z);
if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
{
  error ("possible hess() problems");
}

printf("\tpassed hess test...\n");

//
//-------------------------- Test lyap () ------------------------------
//

lyap = function ( A, B, C )
{
  local (A, B, C)

  if (!exist (B))
  {
    B = A';    // Solve the special form: A*X + X*A' = -C
  }

  if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
    error ("Dimensions do not agree.");
  }

  //
  // Schur decomposition on A and B
  //

  sa = schur (A);
  sb = schur (B);

  //
  // transform C
  //

  tc = sa.z' * C * sb.z;

  X = sylv (sa.t, sb.t, tc);

  //
  // Undo the transformation
  //

  X = sa.z * X * sb.z';

  return X;
};

a = randsvd (10,10);
b = rand (10,10);
c = rand (10,10);

x = lyap (a, b, c);
if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
{
  error ("possible problems with lyap() or sylv()");
}

printf("\tpassed lyap test...\n");

//
//-------------------------- Test eig () ------------------------------
//

trace = function(m)
{
  local(i, tr);

  if(m.class != "num") {
    error("must provide NUMERICAL input to trace()");
  }

  tr = 0;
  for(i in 1:min( [m.nr, m.nc] )) {
    tr = tr + m[i;i];
  }

  return tr;
};

eye = function( m , n )
{
  local(i, N, new);

  if (!exist (n))
  {
    if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
    new = zeros (m[1], m[2]);
    N = min ([m[1], m[2]]);
  else
    if (class (m) == "string" || class (n) == "string") {
      error ("eye(), string arguments not allowed");
    }
    if (max (size (m)) == 1 && max (size (n)) == 1)
    {
      new = zeros (m[1], n[1]);
      N = min ([m[1], n[1]]);
    else
      error ("matrix arguments to eye() must be 1x1");
    }
  }
  for(i in 1:N)
  {
    new[i;i] = 1.0;
  }
  return new;
};

//
// Standard eigenvalue problem
//

a = randsvd(10,10);
ta = trace (a);
sa = symm (a);
tsa = trace (sa);

z = randsvd(10,10) + randsvd(10,10)*1i;
tz = trace (z);
sz = symm (z);
tsz = trace (sz);

tol = 1.e-6;

if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
{
  error ("error in eig");
}
if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
{
  error ("error in eig");
}
if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
{
  error ("error in eig");
}
if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
{
  error ("error in eig");
}

//
// Generalized eigenvalue problem
//

b = randsvd(10,10);
sb = symm (b) + eye(size(b))*3;
tb = trace (b);
tsb = trace (sb);

zb = randsvd(10,10) + randsvd(10,10)*1i;
szb = symm (zb) + eye(size(zb))*3;
tzb = trace (zb);
tszb = trace (szb);

eig(a,b);    // not sure of a good way to check these yet
eigs(sa,sb);
eigs(sa,sb);

eig(z, zb);
eigs(sz, szb);
eigs(sz, szb);

printf("\tpassed eig test...\n");

//
//-------------------------- Test fft () -----------------------------
//

if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
printf("\tpassed fft test...\n");

//
//------------------------- Fibonacci Test -------------------------
//
//  Calculate Fibonacci numbers
//

i=1;
while ( i < 2 ) {
  i=i+1;
  a=0; b=1;
  while ( b < 10000 ) {
    c = b;
    b = a+b;
    a = c;
  }
}
if ( b != 10946 ) {
  error("failed fibonacci test");
else
  printf("\tpassed fibonacci test...\n");
}

//
//------------------------- Factorial Test -------------------------
//

fac = function(a)
{
  if(a <= 1)
  {
    return 1
  else
    return a*$self(a-1)
  }
};

if(fac(10) != 3628800)
{
  error(); else printf("\tpassed factorial test...\n");
}

//
//--------------------------- ACK Test ----------------------------
//

ack = function(a, b)
{
  if(a == 0) { return b + 1; }
  if(b == 0) { return $self(a-1, 1); }
  return $self(a-1, $self(a, b-1));
};

if(ack(2,2) != 7)
{
  error("error in ack() test");
  else
  printf("\tpassed ACK test...\n");
}

//
//------------------------- Prime Test -----------------------------
//
// An example that finds all primes less than limit
//

primes = function (limit)
{
  local(prime, cnt, i, j, k);

  i = 1; cnt = 0;
  for(k in 2:limit)
  {
    j = 2;
    while(mod(k,j) != 0)
    {
      j++;
    }
    if(j == k)             // Found prime
    {
      cnt++;
      prime[i;1] = k;
      i++;
    }
  }
  return prime;
};

if(max(size(primes(100))) != 25)
{
  error("error in prime test");
  else
  printf("\tpassed prime test...\n");
}

//
//--------------------------- Fib Min Test -----------------------------
//
//    fibmin() will minimize an arbitrary function
//    in 1D using Fibonacci search

f065 = function(x)
{
  return (x - 0.65) * (x - 0.65);
};

fib = function(x)
{
  local (n, a, b);

  a = 1;
  b = 1;
  if (x >= 2)
  {
    n = x - 1;
    for (n in n:1:-1)
    {
      c = b;
      b = a + b;
      a = c;
      n = n - 1;
    }
  }
  return b;
};

//  Minimize a 1D function using Fibonacci search
//  f = function to minimize
//  xlo = lower bound
//  xhi = upper bound
//  n = number of iterations (the bigger the more accurate)

fibmin = function(f, xlo, xhi, n)
{
  local(a, b, x, y, ex, ey, k, lo, hi);

  lo = xlo;
  hi = xhi;
  k = n;
  for (k in k:2:-1)
  {
    a = fib(k - 2) / fib(k);
    b = fib(k - 1) / fib(k);
    x = lo + (hi - lo) * a;
    y = lo + (hi - lo) * b;
    ex = f(x);
    ey = f(y);
    if (ex >= ey)
    {
      lo = x;
      else
      hi = y;
    }
    //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  }
  return (lo + hi) / 2;
};

//
// Simple example using above function to mimize f065. Answer is 0.65
//

x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
if (abs(x - 0.65) > 1e-6)
{
  printf("x = %f\n", x);
  error("failed fibmin test");
}

printf("\tpassed fibmin test...\n");

//
//--------------------- Nasty Function Test ------------------------
//

printf("\tStarting Nasty Function Test...");
printf("\tthis will take awhile\n");
check = function( a, b, c, d, e, f, g, h ) {
  if ( a+b+c+d == e+f+g+h && ...
       a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
       a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
    return 1;
  else
    return 0;
  }
};

for(a in 8:10) {
  for(b in 7:(a-1)) {
    for(c in 6:(b-1)) {
      for(d in 5:(c-1)) {
        for(e in 4:(d-1)) {
          for(f in 3:(e-1)) {
            for(g in 2:(f-1)) {
              for(h in 1:(g-1)) {
              if(check( a, b, c, d,  e, f, g, h ) || ...
                     check( a, e, c, d,  b, f, g, h ) || ...
                     check( a, f, c, d,  e, b, g, h ) || ...
                     check( a, g, c, d,  e, f, b, h ) || ...
                     check( a, h, c, d,  e, f, g, b ) || ...
                     check( a, b, e, d,  c, f, g, h ) || ...
                     check( a, b, f, d,  e, c, g, h ) || ...
                     check( a, b, g, d,  e, f, c, h ) || ...
                     check( a, b, h, d,  e, f, g, c ) || ...
                     check( a, b, c, e,  d, f, g, h ) || ...
                     check( a, b, c, f,  e, d, g, h ) || ...
                     check( a, b, c, g,  e, f, d, h ) || ...
                     check( a, b, c, h,  e, f, g, d ) || ...
                     check( a, e, f, d,  b, c, g, h ) || ...
                     check( a, e, g, d,  b, f, c, h ) || ...
                     check( a, e, h, d,  b, f, g, c ) || ...
                     check( a, f, g, d,  e, b, c, h ) || ...
                     check( a, f, h, d,  e, b, g, c ) || ...
                     check( a, g, h, d,  e, f, b, c ) || ...
                     check( a, b, e, f,  c, d, g, h ) || ...
                     check( a, b, e, g,  c, f, d, h ) || ...
                     check( a, b, e, h,  c, f, g, d ) || ...
                     check( a, b, f, g,  e, c, d, h ) || ...
                     check( a, b, f, h,  e, c, g, d ) || ...
                     check( a, b, g, h,  e, f, c, d ) || ...
                     check( a, e, f, g,  e, f, g, h ) || ...
                     check( a, e, f, h,  e, f, g, h ) || ...
                     check( a, e, g, h,  e, f, g, h ) || ...
                     check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
              }
            }
          }
        }
      }
    }
  }
}

if(1) {  // figure out the value of cnt, and check!
  printf("\tpassed nasty function test...\n");
else
  error();
}

//
//------------------ Test More Advanced Functions --------------------
//

printf( "\tStarting the lqr/ode test..." );
printf( "\tthis will take awhile\n" );

lqr = function( a, b, q, r )
{
  local( k, s,...
         m, n, mb, nb, mq, nq,...
         e, v, d );

  m = size(a)[1]; n = size(a)[2];
  mb = size(b)[1]; nb = size(b)[2];
  mq = size(q)[1]; nq = size(q)[2];

  if ( m != mq || n != nq )
  {
    fprintf( "stderr", "A and Q must be the same size.\n" );
    quit
  }

  mr = size(r)[1]; nr = size(r)[2];
  if ( mr != nr || nb != mr )
  {
    fprintf( "stderr", "B and R must be consistent.\n" );
    quit
  }

  nn = zeros( m, nb );

  // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:

  e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  v = e.vec; d = e.val;

  index = sort( real( d ) ).ind;
  d = real( d[ index ] );

  if ( !( d[n] < 0 && d[n+1] > 0 ) )
  {
    fprintf( "stderr", "Can't order eigenvalues.\n" );
    quit
  }

  chi = v[ 1:n; index[1:n] ];
  lambda = v[ (n+1):(2*n); index[1:n] ];
  s = -real(solve(chi',lambda')');
  k = solve( r, nn'+b'*s );

  return << k=k; s=s >>;

};

// Now run a little test problem.

k = 1; m = 1; c = .1;
a = [0     ,1    ,0    , 0;
    -k/m, -c/m,  k/m,  c/m;
     0,     0,    0,    1;
     k/m,  c/m, -k/m, -c/m ];
b = [ 0; 1/m; 0; 0 ];
qxx = diag( [0, 0, 100, 0] );
ruu = [1];
K = lqr( a, b, qxx, ruu ).k;

dot = function( t, x )
{
  global (a, b, K)
  return (a-b*K)*x + b*K*([1,0,1,0]');
};

x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );

m = maxi( x[;2] );

if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
     any (abs( x[x.nr;2,4] - 1 ) > 0.001) )
{
  printf( "\tfailed***\n" );
  else
  printf( "\tpassed the lqr/ode test...\n" );
}

printf("Elapsed time = %10.3f seconds\n", toc() );
"FINISHED TESTS"