home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / test < prev    next >
Encoding:
Text File  |  1995-01-26  |  69.5 KB  |  3,119 lines  |  [TEXT/ttxt]

  1. //
  2. // Test Suite for RLaB.
  3. // This test suite is supposed to test as much as possible, and still
  4. // be runable on most platforms. This means we cannot do graphics or
  5. // other platform specific stuff like piping. However, that is OK,
  6. // since we are mostly interested in assuring the builder that RLaB
  7. // will produce correct numerical answers.
  8.  
  9. // We use randsvd, which in turn uses RLaB's random number
  10. // generator. We try to be carefull and seed the random number
  11. // generator so that each user will get similar results (or at least
  12. // similar inputs to the tests).
  13.  
  14. // Test the other style comments 
  15. 1 + 2; // A simple statement with trailing comment.
  16. #  Optional RLaB comment style.
  17. 1 + 2; #  A simple statement with trailing comment.
  18. %  Optional Matlab comment style.
  19. 1 + 2; %  A simple statement with trailing comment.
  20.  
  21. srand (SEED = 10);   // Seems to produce reasonable results
  22. rand("default");
  23.  
  24. tic();    // Start timing the tests...
  25.  
  26. //
  27. // Test Parameters and some functions we will need later
  28. //
  29.  
  30. pi = 4.0*atan(1.0);
  31. X = 3;    // should be 3 (heuristic).
  32.  
  33. //
  34. // Compute machine epsilon
  35. //
  36.  
  37. epsilon = function() 
  38. {
  39.   eps = 1.0;
  40.   while((eps + 1.0) != 1.0) 
  41.   {
  42.     eps = eps/2.0;
  43.   }
  44.   return 2*eps;
  45. };
  46.  
  47. eps = epsilon();
  48.  
  49. eye = function( m , n ) 
  50. {
  51.   if (!exist (n))
  52.   {
  53.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  54.     new = zeros (m[1], m[2]);
  55.     N = min ([m[1], m[2]]);
  56.   else
  57.     if (class (m) == "string" || class (n) == "string") {
  58.       error ("eye(), string arguments not allowed");
  59.     }
  60.     if (max (size (m)) == 1 && max (size (n)) == 1)
  61.     {
  62.       new = zeros (m[1], n[1]);
  63.       N = min ([m[1], n[1]]);
  64.     else
  65.       error ("matrix arguments to eye() must be 1x1");
  66.     }
  67.   }
  68.   for(i in 1:N)
  69.   {
  70.     new[i;i] = 1.0;
  71.   }
  72.   return new;
  73. };
  74.  
  75. symm = function( A )
  76. {
  77.   return (A + A')./2;
  78. };
  79.  
  80. //-------------------------------------------------------------------//
  81.  
  82. // Synopsis:    Pascal matrix.
  83.  
  84. // Syntax:    P = pascal ( N )
  85.  
  86. // Description:
  87.  
  88. //    The Pascal matrix of order N: a symmetric positive definite
  89. //    matrix with integer entries taken from Pascal's triangle. The
  90. //    Pascal matrix is totally positive and its inverse has integer
  91. //    entries.  Its eigenvalues occur in reciprocal pairs. COND(P)
  92. //    is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
  93. //    lower triangular Cholesky factor (up to signs of columns) of
  94. //    the Pascal matrix.   It is involutary (is its own
  95. //    inverse). PASCAL(N,2) is a transposed and permuted version of
  96. //    PASCAL(N,1) which is a cube root of the identity.
  97.  
  98. //      References:
  99. //      R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
  100. //           Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
  101. //           gives a factorization of L = PASCAL(N,1) and a formula for the
  102. //           elements of L^k).
  103. //      S. Karlin, Total Positivity, Volume 1, Stanford University Press,
  104. //           1968.  (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
  105. //                   PASCAL(N) is a submatrix of this matrix.)
  106. //      M. Newman and J. Todd, The evaluation of matrix inversion programs,
  107. //           J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
  108. //      H. Rutishauser, On test matrices, Programmation en Mathematiques
  109. //           Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  110. //           1966, pp. 349-365.  (Gives an integral formula for the
  111. //           elements of PASCAL(N).)
  112. //      J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  113. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
  114. //      H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
  115. //           Blackie, London and Glasgow, 1929.  (PASCAL(N,2) on page 332.)
  116.  
  117. //    This file is a translation of pascal.m from version 2.0 of
  118. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  119. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  120.  
  121. //-------------------------------------------------------------------//
  122.  
  123. pascal = function ( n , k )
  124. {
  125.   local(n, k)
  126.  
  127.   if (!exist (k)) { k = 0; }
  128.  
  129.   P = diag( (-1).^[0:n-1] );
  130.   P[; 1] = ones(n,1);
  131.  
  132.   //  Generate the Pascal Cholesky factor (up to signs).
  133.  
  134.   for (j in 2:n-1)
  135.   {
  136.     for (i in j+1:n)
  137.     {
  138.       P[i;j] = P[i-1;j] - P[i-1;j-1];
  139.     }
  140.   }
  141.  
  142.   if (k == 0)
  143.   {
  144.     P = P*P';
  145.   else if (k == 2) {
  146.     P = rot90(P,3);
  147.     if (n/2 == round(n/2)) { P = -P; }
  148.   }}
  149.  
  150.   return P;
  151. };
  152.  
  153. //-------------------------------------------------------------------//
  154.  
  155. // Synopsis:    Random, orthogonal upper Hessenberg matrix.
  156.  
  157. // Syntax:      H = ohess ( N )
  158.  
  159. // Description:
  160.  
  161. //      H is an N-by-N real, random, orthogonal upper Hessenberg
  162. //      matrix. Alternatively, H = OHESS(X), where X is an arbitrary
  163. //      real N-vector (N > 1) constructs H non-randomly using the
  164. //      elements of X as parameters. In both cases H is constructed
  165. //      via a product of N-1 Givens rotations. 
  166.  
  167. //      Note: See Gragg (1986) for how to represent an N-by-N
  168. //      (complex) unitary Hessenberg matrix with positive subdiagonal
  169. //      elements in terms of 2N-1 real parameters (the Schur
  170. //      parametrization). This M-file handles the real case only and
  171. //      is intended simply as a convenient way to generate random or
  172. //      non-random orthogonal Hessenberg matrices.
  173.  
  174. //      Reference:
  175. //      W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
  176. //      J. Comp. Appl. Math., 16 (1986), pp. 1-8.
  177.  
  178. //    This file is a translation of ohess.m from version 2.0 of
  179. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  180. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  181.  
  182. //-------------------------------------------------------------------//
  183.  
  184. ohess = function ( x )
  185. {
  186.   local (x)
  187.   global (pi)
  188.  
  189.   if (any (imag (x))) { error("Parameter must be real."); }
  190.  
  191.   n = max(size(x));
  192.  
  193.   if (n == 1)
  194.   {
  195.     //  Handle scalar x.
  196.     n = x;
  197.     x = rand(n-1, 1)*2*pi;
  198.     H = eye(n,n);
  199.     H[n;n] = sign(rand());
  200.   else
  201.     H = eye(n,n);
  202.     H[n;n] = sign(x[n]) + (x[n]==0);   // Second term ensures H[n;n] nonzero.
  203.   }
  204.  
  205.   for (i in n:2:-1)
  206.   {
  207.     // Apply Givens rotation through angle x(i-1).
  208.     theta = x[i-1];
  209.     c = cos(theta);
  210.     s = sin(theta);
  211.     H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
  212.                        -s*H[i-1;]+c*H[i;] ];
  213.   }
  214.  
  215.   return H;
  216. };
  217.  
  218. //-------------------------------------------------------------------//
  219.  
  220. // Synopsis:    Random matrix with pre-assigned singular values.
  221.  
  222. // Syntax:    R = randsvd (N, KAPPA, MODE, KL, KU)
  223.  
  224. // Description:
  225.  
  226. //    R is a (banded) random matrix of order N with COND(A) = KAPPA
  227. //    and singular values from the distribution MODE.
  228.  
  229. //      N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
  230. //      Available types:
  231. //             MODE = 1:   one large singular value,
  232. //             MODE = 2:   one small singular value,
  233. //             MODE = 3:   geometrically distributed singular values,
  234. //             MODE = 4:   arithmetically distributed singular values,
  235. //             MODE = 5:   random singular values with unif. dist. logarithm.
  236.  
  237. //      If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
  238. //      If MODE < 0 then the effect is as for ABS(MODE) except that in the
  239. //      original matrix of singular values the order of the diagonal entries
  240. //      is reversed: small to large instead of large to small.
  241. //      KL and KU are the lower and upper bandwidths respectively; if they
  242. //      are omitted a full matrix is produced.
  243. //      If only KL is present, KU defaults to KL.
  244. //      Special case: if KAPPA < 0 then a random full symmetric positive
  245. //                    definite matrix is produced with COND(A) = -KAPPA and
  246. //                    eigenvalues distributed according to MODE.
  247. //                    KL and KU, if present, are ignored.
  248.  
  249. //      This routine is similar to the more comprehensive Fortran routine xLATMS
  250. //      in the following reference:
  251. //      J.W. Demmel and A. McKenney, A test matrix generation suite,
  252. //      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
  253. //      New York, 1989.
  254.  
  255. //    This file is a translation of randsvd.m from version 2.0 of
  256. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  257. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  258.  
  259. // Dependencies
  260. #   rfile bandred
  261. #   rfile qmult
  262.  
  263. //-------------------------------------------------------------------//
  264.  
  265. randsvd = function (n, kappa, mode, kl, ku)
  266. {
  267.   local (n, kappa, mode, kl, ku)
  268.  
  269.   if (!exist (kappa)) { kappa = sqrt(1/eps); }
  270.   if (!exist (mode))  { mode = 3; }
  271.   if (!exist (kl))    { kl = n-1; }    // Full matrix.
  272.   if (!exist (ku))    { ku = kl; }     // Same upper and lower bandwidths.
  273.  
  274.   if (abs(kappa) < 1) {
  275.     error("Condition number must be at least 1!");
  276.   }
  277.  
  278.   posdef = 0;
  279.   if (kappa < 0)            // Special case.
  280.   {
  281.     posdef = 1;
  282.     kappa = -kappa;
  283.   }
  284.  
  285.   p = min(n);
  286.   m = n[1];        // Parameter n specifies dimension: m-by-n.
  287.   n = n[max(size(n))];
  288.  
  289.   if (p == 1)        // Handle case where A is a vector.
  290.   {
  291.     rand("normal", -10, 10);
  292.     A = rand(m, n);
  293.     A = A/norm(A, "2");
  294.     return A;
  295.   }
  296.  
  297.   j = abs(mode);
  298.  
  299.   // Set up vector sigma of singular values.
  300.   if (j == 3)
  301.   {
  302.     factor = kappa^(-1/(p-1));
  303.     sigma = factor.^[0:p-1];
  304.  
  305.   else if (j == 4) {
  306.     sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
  307.  
  308.   else if (j == 5) {    // In this case cond(A) <= kappa.
  309.     rand("uniform", 0, 1)
  310.     sigma = exp( -rand(p,1)*log(kappa) );
  311.  
  312.   else if (j == 2) {
  313.     sigma = ones(p,1);
  314.     sigma[p] = 1/kappa;
  315.  
  316.   else if (j == 1) {
  317.     sigma = ones(p,1)./kappa;
  318.     sigma[1] = 1;
  319.  
  320.   }}}}}
  321.  
  322.  
  323.   // Convert to diagonal matrix of singular values.
  324.   if (mode < 0) {
  325.     sigma = sigma[p:1:-1];
  326.   }
  327.  
  328.   sigma = diag(sigma);
  329.  
  330.   if (posdef)        // Handle special case.
  331.   {
  332.     Q = qmult(p);
  333.     A = Q'*sigma*Q;
  334.     A = (A + A')/2;    // Ensure matrix is symmetric.
  335.     return A;
  336.   }
  337.  
  338.   if (m != n) 
  339.   {
  340.     sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  341.   }
  342.  
  343.   if (kl == 0 && ku == 0) // Diagonal matrix requested - nothing more to do.
  344.   {
  345.     A = sigma;
  346.     return A;
  347.   }
  348.  
  349.   // A = U*sigma*V, where U, V are random orthogonal matrices from the
  350.   // Haar distribution.
  351.   A = qmult(sigma');
  352.   A = qmult(A');
  353.  
  354.   if (kl < n-1 || ku < n-1)    // Bandwidth reduction.
  355.   {
  356.    A = bandred(A, kl, ku);
  357.   }
  358.  
  359.   rand("default");
  360.   return A;
  361. };
  362.  
  363. //-------------------------------------------------------------------//
  364.  
  365. // Synopsis:    Band reduction by two-sided unitary transformations.
  366.  
  367. // Syntax:    bandred ( A , KL, KU )
  368.  
  369. // Description:
  370.  
  371. //    bandred(A, KL, KU) is a matrix unitarily equivalent to A with
  372. //    lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
  373. //    i > j+KL or j > i+KU).  The reduction is performed using
  374. //    Householder transformations. If KU is omitted it defaults to
  375. //    KL. 
  376.  
  377. //    Called by randsvd.
  378. //    This is a `standard' reduction.  Cf. reduction to bidiagonal
  379. //    form prior to computing the SVD.  This code is a little
  380. //    wasteful in that it computes certain elements which are
  381. //    immediately set to zero! 
  382. //
  383. //      Reference:
  384. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  385. //      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  386. //      Section 5.4.3.
  387.  
  388. //    This file is a translation of bandred.m from version 2.0 of
  389. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  390. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  391.  
  392. // Dependencies
  393. #   rfile house
  394.  
  395. //-------------------------------------------------------------------//
  396.  
  397. bandred = function ( A , kl , ku )
  398. {
  399.   local (A, kl, ku)
  400.  
  401.   if (!exist (ku)) { ku = kl; else ku = ku; }
  402.  
  403.   if (kl == 0 && ku == 0) {
  404.     error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  405.   }
  406.  
  407.   // Check for special case where order of left/right transformations matters.
  408.   // Easiest approach is to work on the transpose, flipping back at the end.
  409.  
  410.   flip = 0;
  411.   if (ku == 0)
  412.   {
  413.     A = A';
  414.     temp = kl; kl = ku; ku = temp; flip = 1;
  415.   }
  416.  
  417.   m = A.nr; n = A.nc; 
  418.  
  419.   for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  420.   {
  421.     if (j+kl+1 <= m)
  422.     {
  423.        ltmp = house(A[j+kl:m;j]);
  424.        beta = ltmp.beta; v = ltmp.v;
  425.        temp = A[j+kl:m;j:n];
  426.        A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
  427.        A[j+kl+1:m;j] = zeros(m-j-kl,1);
  428.     }
  429.  
  430.     if (j+ku+1 <= n)
  431.     {
  432.        ltmp = house(A[j;j+ku:n]');
  433.        beta = ltmp.beta; v = ltmp.v;
  434.        temp = A[j:m;j+ku:n];
  435.        A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
  436.        A[j;j+ku+1:n] = zeros(1,n-j-ku);
  437.     }
  438.   }
  439.  
  440.   if (flip) {
  441.     A = A';
  442.   }
  443.  
  444.   return A;
  445. };
  446.  
  447. //-------------------------------------------------------------------//
  448.  
  449. // Synopsis:    Lehmer matrix - symmetric positive definite.
  450.  
  451. // Syntax:      A = lehmer ( N )
  452.  
  453. // Description:
  454.  
  455. //      A is the symmetric positive definite N-by-N matrix with
  456. //                     A(i,j) = i/j for j >= i.
  457. //      A is totally nonnegative.  INV(A) is tridiagonal, and explicit
  458. //      formulas are known for its entries. 
  459.  
  460. //      N <= COND(A) <= 4*N*N.
  461.  
  462. //      References:
  463. //        M. Newman and J. Todd, The evaluation of matrix inversion
  464. //           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
  465. //        Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
  466. //           of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
  467. //        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  468. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.
  469.  
  470. //    This file is a translation of lehmer.m from version 2.0 of
  471. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  472. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  473.  
  474. //-------------------------------------------------------------------//
  475.  
  476. lehmer = function ( n )
  477. {
  478.   local (n)
  479.   global (tril)
  480.  
  481.   A = ones(n,1)*(1:n);
  482.   A = A./A';
  483.   A = tril(A) + tril(A,-1)';
  484.  
  485.   return A;
  486. };
  487.  
  488. //-------------------------------------------------------------------//
  489.  
  490. // Synopsis:    Pre-multiply by random orthogonal matrix.
  491.  
  492. // Syntax:    B = qmult ( A )
  493.  
  494. // Description:
  495.  
  496. //    B is Q*A where Q is a random real orthogonal matrix from the
  497. //    Haar distribution, of dimension the number of rows in
  498. //    A. Special case: if A is a scalar then QMULT(A) is the same as
  499.  
  500. //        qmult(eye(a))
  501.  
  502. //       Called by RANDSVD.
  503.  
  504. //       Reference:
  505. //       G.W. Stewart, The efficient generation of random
  506. //       orthogonal matrices with an application to condition estimators,
  507. //       SIAM J. Numer. Anal., 17 (1980), 403-409.
  508.  
  509. //    This file is a translation of qmult.m from version 2.0 of
  510. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  511. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  512.  
  513. //-------------------------------------------------------------------//
  514.  
  515. qmult = function ( A )
  516. {
  517.   local (A)
  518.  
  519.   n = A.nr; m = A.nc;
  520.  
  521.   //  Handle scalar A.
  522.   if (max(n,m) == 1)
  523.   {
  524.     n = A;
  525.     A = eye(n,n);
  526.   }
  527.  
  528.   d = zeros(n,n);
  529.  
  530.   for (k in n-1:1:-1)
  531.   {
  532.     // Generate random Householder transformation.
  533.     rand("normal", 0, 1);
  534.     x = rand(n-k+1,1);
  535.     s = norm(x, "2");
  536.     sgn = sign(x[1]) + (x[1]==0);    // Modification for sign(1)=1.
  537.     s = sgn*s;
  538.     d[k] = -sgn;
  539.     x[1] = x[1] + s;
  540.     beta = s*x[1];
  541.  
  542.     // Apply the transformation to A.
  543.     y = x'*A[k:n;];
  544.     A[k:n;] = A[k:n;] - x*(y/beta);
  545.   }
  546.  
  547.   // Tidy up signs.
  548.   for (i in 1:n-1)
  549.   {
  550.     A[i;] = d[i]*A[i;];
  551.   }
  552.  
  553.   A[n;] = A[n;]*sign(rand());
  554.   B = A;
  555.  
  556.   rand("default");
  557.   return B;
  558. };
  559.  
  560. //-------------------------------------------------------------------//
  561.  
  562. // Synopsis:    Create a Householder matrix.
  563.  
  564. // Syntax:    house ( X )
  565.  
  566. //    If HOUSE(x), which returns a list containing elements `v' and
  567. //    `beta',  then H = EYE - beta*v*v' is a Householder matrix such
  568. //    that Hx = -sign(x(1))*norm(x)*e_1. 
  569. //    NB: If x = 0 then v = 0, beta = 1 is returned.
  570. //            x can be real or complex.
  571. //            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
  572.  
  573. //    Theory: (textbook references Golub & Van Loan 1989, 38-43;
  574. //             Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
  575. //      Hx = y: (I - beta*v*v')x = -s*e_1.
  576. //      Must have |s| = norm(x), v = x+s*e_1, and
  577. //      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
  578. //      So take s = sign(x(1))*norm(x) (which avoids cancellation).
  579. //      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
  580. //          = 2*norm(x)*(norm(x) + |x(1)|).
  581.  
  582. //    References:
  583. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  584. //           Johns Hopkins University Press, Baltimore, Maryland, 1989.
  585. //        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  586. //           New York, 1973,
  587. //        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  588. //           Press, 1965.
  589.  
  590. //    This file is a translation of house.m from version 2.0 of
  591. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  592. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  593.  
  594. //-------------------------------------------------------------------//
  595.  
  596. house = function ( x )
  597. {
  598.   local (x)
  599.  
  600.   m = x.nr; n = x.nc;
  601.   if (n > 1) { error ("Argument must be a column vector."); }
  602.  
  603.   s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  604.   v = x;
  605.   if (s == 0)    // Quit if x is the zero vector.
  606.   {
  607.     beta = 1; 
  608.     return << beta = beta ; v = v >>;
  609.   }
  610.  
  611.   v[1] = v[1] + s;
  612.   beta = 1/(s'*v[1]);        // NB the conjugated s.
  613.  
  614.   // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  615.   // But beta as above can be non-real (due to rounding) only 
  616.   // when x is complex.
  617.  
  618.   return << beta = beta ; v = v >>;
  619. };
  620.  
  621. //-------------------------------------------------------------------//
  622.  
  623. //  Syntax:    tril ( A )
  624. //        tril ( A , K )
  625.  
  626. //  Description:
  627.  
  628. //  tril(x) returns the lower triangular part of A.
  629.  
  630. //  tril(A,K) returns the elements on and below the K-th diagonal of
  631. //  A.
  632.  
  633. //  K = 0: main diagonal
  634. //  K > 0: above the main diag.
  635. //  K < 0: below the main diag.
  636.  
  637. //  See Also: triu
  638. //-------------------------------------------------------------------//
  639.  
  640. tril = function(x, k) 
  641. {
  642.   local(x, k)
  643.  
  644.   if (!exist (k)) { k = 0; }
  645.   nr = x.nr; nc = x.nc;
  646.   if(k > 0) 
  647.   { 
  648.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  649.   else
  650.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  651.   }
  652.  
  653.   y = zeros(nr, nc);
  654.  
  655.   for(i in max( [1,1-k] ):nr) {
  656.     j = 1:min( [nc, i+k] );
  657.     y[i;j] = x[i;j];
  658.   }
  659.  
  660.   return y;
  661. };
  662.  
  663. //-------------------------------------------------------------------//
  664.  
  665. //  Syntax:    triu ( A )
  666. //        triu ( A , K )
  667.  
  668. //  Description:
  669.  
  670. //  triu(x) returns the upper triangular part of A.
  671.  
  672. //  tril(x; k) returns the elements on and above the k-th diagonal of
  673. //  A. 
  674.  
  675. //  K = 0: main diagonal
  676. //  K > 0: above the main diag.
  677. //  K < 0: below the main diag.
  678.  
  679. //  See Also: tril
  680. //-------------------------------------------------------------------//
  681.  
  682. triu = function(x, k) 
  683. {
  684.   local(x, k)
  685.  
  686.   if (!exist (k)) { k = 0; }
  687.   nr = x.nr; nc = x.nc;
  688.  
  689.   if(k > 0) 
  690.   { 
  691.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  692.   else
  693.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  694.   }
  695.  
  696.   y = zeros(nr, nc);
  697.  
  698.   for(j in max( [1,1+k] ):nc) {
  699.     i = 1:min( [nr, j-k] );
  700.     y[i;j] = x[i;j];
  701.   }
  702.  
  703.   return y;
  704. };
  705.  
  706. //
  707. //-------------------- Test Relational Expressions -------------------
  708. //
  709. printf("\tstart scalar tests...\n");
  710. printf("\tstart relational tests...\n");
  711.  
  712. //    SCALAR CONSTANTS (REAL)
  713. if( !(1<2) ) { error(); }
  714. if( !(1<=2) ) { error(); }
  715. if( 1>2 ) { error(); }
  716. if( 1>=2 ) { error(); }
  717. if( 1==2 ) { error(); }
  718. if( !(1!=2) ) { error(); }
  719. if( !1 ) { error(); }
  720. if( !!!1) { error(); }
  721.  
  722. if( !([1]<[2]) ) { error(); }
  723. if( !([1]<=[2]) ) { error(); }
  724. if( [1]>[2] ) { error(); }
  725. if( [1]>=[2] ) { error(); }
  726. if( [1]==[2] ) { error(); }
  727. if( !([1]!=[2]) ) { error(); }
  728. if( ![1] ) { error(); }
  729. if( !!![1]) { error(); }
  730.  
  731. //    SCALAR CONSTANTS (COMPLEX)
  732. if( !(1+2i<2+3i) ) { error(); }
  733. if( !(1+2i<=2+3i) ) { error(); }
  734. if( 1+2i>2+3i ) { error(); }
  735. if( 1+2i>=2+3i ) { error(); }
  736. if( 1+2i==2+3i ) { error(); }
  737. if( !(1+2i!=2+3i) ) { error(); }
  738. if( !1+2i ) { error(); }
  739. if( !!!1+2i) { error(); }
  740.  
  741. if( !([1+2i]<[2+3i]) ) { error(); }
  742. if( !([1+2i]<=[2+3i]) ) { error(); }
  743. if( [1+2i]>[2+3i] ) { error(); }
  744. if( [1+2i]>=[2+3i] ) { error(); }
  745. if( [1+2i]==[2+3i] ) { error(); }
  746. if( !([1+2i]!=[2+3i]) ) { error(); }
  747. if( ![1+2i] ) { error(); }
  748. if( !!![1+2i]) { error(); }
  749.  
  750. //    SCALAR ENTITIES (REAL)
  751. a=1;b=2;
  752. if( !(a<b) ) { error(); }
  753. if( !(a<=b) ) { error(); }
  754. if( a>b ) { error(); }
  755. if( a>=b ) { error(); }
  756. if( a==b ) { error(); }
  757. if( !(a!=b) ) { error(); }
  758. if( !a ) { error(); }
  759. if( !!!a) { error(); }
  760.  
  761. if( !([a]<[b]) ) { error(); }
  762. if( !([a]<=[b]) ) { error(); }
  763. if( [a]>[b] ) { error(); }
  764. if( [a]>=[b] ) { error(); }
  765. if( [a]==[b] ) { error(); }
  766. if( !([a]!=[b]) ) { error(); }
  767. if( ![a] ) { error(); }
  768. if( !!![a]) { error(); }
  769.  
  770. //    SCALAR ENTITIES (COMPLEX)
  771. a=1+2i;b=2+3i;
  772. if( !(a<b) ) { error(); }
  773. if( !(a<=b) ) { error(); }
  774. if( a>b ) { error(); }
  775. if( a>=b ) { error(); }
  776. if( a==b ) { error(); }
  777. if( !(a!=b) ) { error(); }
  778. if( !a ) { error(); }
  779. if( !!!a) { error(); }
  780.  
  781. if( !([a]<[b]) ) { error(); }
  782. if( !([a]<=[b]) ) { error(); }
  783. if( [a]>[b] ) { error(); }
  784. if( [a]>=[b] ) { error(); }
  785. if( [a]==[b] ) { error(); }
  786. if( !([a]!=[b]) ) { error(); }
  787. if( ![a] ) { error(); }
  788. if( !!![a]) { error(); }
  789.  
  790. if (! all (all (-2 <  rand(4,4)))) { error (); }
  791. if (! all (all (-2 <= rand(4,4)))) { error (); }
  792. if (! all (all ( 2 >  rand(4,4)))) { error (); }
  793. if (! all (all ( 2 >= rand(4,4)))) { error (); }
  794.  
  795. if (! all (all (rand(4,4) >  -2))) { error (); }
  796. if (! all (all (rand(4,4) >= -2))) { error (); }
  797. if (! all (all (rand(4,4) <   2))) { error (); }
  798. if (! all (all (rand(4,4) <=  2))) { error (); }
  799.  
  800. if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
  801. if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
  802. if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
  803. if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
  804.  
  805. //------------------------- Test REAL SCALARS ------------------------
  806. //
  807. //    CONSTANTS
  808. //      Addition
  809. if(1+2 != 3) { error(); }
  810. //      Subtraction
  811. if(1-2 != -1) { error(); }
  812. //      Multiply
  813. if(1*2 != 2) { error(); }
  814. //      Divide
  815. if(1/2 != 0.5) { error(); }
  816. //      Power
  817. if(2^2 != 4) { error(); }
  818. if(4^0 != 1) { error(); }
  819. //      Unary Minus
  820. if(2-3 != -1) { error(); }
  821. //
  822. //    ENTITIES
  823. //
  824. a = 1; b = 2; c = 3; d = 0.5;
  825. //      Addition
  826. if(a+b != c) { error(); }
  827. //      Subtraction
  828. if(a-b != -a) { error(); }
  829. //      Multiply
  830. if(a*b != b) { error(); }
  831. //      Divide
  832. if(a/b != d) { error(); }
  833. //      Power
  834. if(b^b != b*b) { error(); }
  835. if(b^0 != 1) { error (); }
  836. //      Unary Minus
  837. if(b-c != -a) { error(); }
  838. //
  839. //  ENTITIES & CONSTANTS
  840. //
  841. if(a+2 != c) { error(); }
  842. //      Subtraction
  843. if(1-b != -a) { error(); }
  844. //      Multiply
  845. if(a*2 != 2) { error(); }
  846. //      Divide
  847. if(1/b != d) { error(); }
  848. //      Power
  849. if(2^b != b*b) { error(); }
  850. //      Unary Minus
  851. if(b-3 != -a) { error(); }
  852. //
  853. //------------------------Test COMPLEX SCALARS -------------------------
  854. //
  855. //    CONSTANTS
  856. if(sqrt(-1) != 1i) { error(); }
  857. //      Addition
  858. if((1+2i)+(2+3i) != (3+5i)) { error(); }
  859. //      Subtraction
  860. if((1+2i)-(3+4i) != (-2-2i)) { error(); }
  861. //      Multiply
  862. if((1+2i)*(3+4i) != (-5+10i)) { error(); }
  863. //      Divide
  864. if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
  865. //      Power
  866. //      Precision problems prevent us from testing these. Have to
  867. //      be checked by hand.
  868. //  (1+2i)^2 = -3 + 4i
  869. //  (1+2i)^.5 = 1.272 + 7.862e-1i
  870. //  if((1+2i)^2 != (-3+4i)) { error(); }
  871. //  if((1+2i)^10 != (237+3116i)) { error(); }
  872. //      Unary Minus
  873. if(-(1+2i) != -1-2i) { error(); }
  874. //
  875. //    ENTITIES
  876. //
  877. a = 1+2i; b = 3+4i; c = 4+6i;
  878. if(a+b != c) { error(); }
  879. //      Subtraction
  880. if(a-b != -2-2i) { error(); }
  881. //      Multiply
  882. if(a*b != -5+10i) { error(); }
  883. //      Divide
  884. //if(a/(3-4i) != -.2+.4i) { error(); }
  885. //      Power
  886. //  if(b^b != b*b) { error(); }
  887. //      Unary Minus
  888. if(-a != -1-2i) { error(); }
  889. //
  890. //    ENTITIES & CONSTANTS
  891. //
  892. if(a+(3+4i) != c) { error(); }
  893. //      Subtraction
  894. if((1+2i)-b != -2-2i) { error(); }
  895. //      Multiply
  896. if(a*(3+4i) != -5+10i) { error(); }
  897. //      Divide
  898. //if(a/(3-4i) != -.2+.4i) { error(); }
  899. //      Power
  900. //if(b^b != b*b) { error(); }
  901. //      Unary Minus
  902. if(-(1+2i) != -1-2i) { error(); }
  903. //
  904. // String - Numerical Equalities
  905. //
  906. if ((1 == "1")) { error(); }
  907. if (([1] == "1")) { error(); }
  908. if (("1" == 1)) { error(); }
  909. if (("1" == [1])) { error(); }
  910.  
  911. if (rand(3,3) == "str") { error(); }
  912. if ("str" == rand(3,3)) { error(); }
  913. if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
  914.  
  915. //
  916. //----------------------- Test REAL MATRICES ---------------------------
  917. //
  918.  
  919. printf("\tstart matrix tests...\n");
  920. printf("\t\treal-matrices\n");
  921.  
  922. //  Read in test matrices
  923. //
  924.  
  925. read("test.input");
  926.  
  927. //
  928. //  Matrix construction
  929. //
  930.  
  931. if(any([1;2;3] != [1,2,3]')) {
  932.   error();
  933. }
  934.  
  935. if(any (any (m0 != zeros(2,2)))) { error(); }
  936. if(any (any (m1 != 1+zeros(2,2)))) { error(); }
  937. if(any (any (m2 != [1,2;3,4]))) { error(); }
  938. if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
  939. if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
  940. if(any (any ([m2] != m2))) { error(); }
  941.  
  942. //
  943. //  Matrix indexing
  944. //
  945.  
  946. p = pascal(6);
  947. if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
  948. if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
  949. if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  950. if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  951. if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }
  952.  
  953. //
  954. //  Sub-Matrix promotion
  955. //
  956.  
  957. if(m2[2;2] != 4) { error(); }
  958. if(any(m2[2;] != [3,4])) { error(); }
  959. if(any(m2[;2] != [2,4]')) { error(); }
  960. i=2;j=1;
  961. if(m2[i;j] != 3) { error(); }
  962. i=1;j=2;
  963. if(m2[i;j] != 2) { error(); }
  964. m = [1,2,3;4,5,6;7,8,9];
  965.  
  966. if(any(m[1;1,2] != [1,2])) 
  967. {
  968.   error();
  969. }
  970.  
  971. if(any(m[1,2;1] != [1;4])) 
  972. {
  973.   error();
  974. }
  975.  
  976. if(any (any (m[1,2;1,2] != [1,2;4,5]))) 
  977. {
  978.   error();
  979. }
  980.  
  981. //
  982.  
  983. if(m3[2;2] != (5+6i)) { error(); }
  984. if(any(m3[2;] != [3+4i,5+6i])) { error(); }
  985. if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
  986.  
  987. //
  988. //  Automatic creation, extension
  989. //
  990.  
  991. if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
  992. a=[1,2,3;4,5,6];
  993. if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
  994. a=[1,2;3,4];
  995. 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]))) 
  996. {
  997.   error();
  998. }
  999.  
  1000. mmm[2;] = a[1;];
  1001.  
  1002. //
  1003. //  Matrix binary operations
  1004. //
  1005.  
  1006. a = m2; b = [5,6;7,8];
  1007. if(any (any (a+a != [2,4;6,8]))) { error(); }
  1008. if(any (any (a-a != zeros(2,2)))) { error(); }
  1009. if(any (any (2+a != [3,4;5,6]))) { error(); }
  1010. if(any (any (2-a != [1,0;-1,-2]))) { error(); }
  1011. if(any (any (a-2 != [-1,0;1,2]))) { error(); }
  1012. if(any (any (2*a != [2,4;6,8]))) { error(); }
  1013. if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
  1014. if(any (any (a*a != [7,10;15,22]))) { error(); }
  1015. if(any (any (a*a*a != [37,54;81,118]))) { error(); }
  1016. if(any (any (a .* a != [1,4;9,16]))) { error(); }
  1017. if(any (any (a./a != [1,1;1,1]))) { error(); }
  1018. if(any (any (a' != [1,3;2,4]))) { error(); }
  1019.  
  1020. if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
  1021. if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
  1022. if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
  1023. if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
  1024. if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
  1025. if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
  1026.  
  1027. if(any ([1;2;3]+[4;5;6] != [5;7;9])) 
  1028. {
  1029.   error();
  1030. }
  1031.  
  1032. if(any ([1;2;3]-[4;5;6] != [-3;-3;-3])) 
  1033. {
  1034.   error();
  1035. }
  1036.  
  1037. if(any ([2;2;2] ./ [1;1;1] != [2;2;2])) 
  1038. {
  1039.   error();
  1040. }
  1041.  
  1042. if(any ([1;2;3] .* [4;5;6] != [4;10;18])) 
  1043. {
  1044.   error();
  1045. }
  1046.  
  1047. if (type (1^.33333) != "real") { error (); }
  1048. if (type (1^(1/3)) != "real") { error (); }
  1049. if (type ([1]^.33333) != "real") { error (); }
  1050. if (type (1^[.33333]) != "real") { error (); }
  1051. if (type ([1]^[.33333]) != "real") { error (); }
  1052.  
  1053. //
  1054. // Test row-wise matrix addition
  1055. //
  1056.  
  1057. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1058. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1059. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1060.  
  1061. a = [1,2,3] + [1,2,3]*1i; 
  1062. 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;
  1063. 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;
  1064. if (!all (all (a .+ b == c))) { error (); }
  1065. if (!all (all (b .+ a == c))) { error (); }
  1066.  
  1067. printf("\t\tpassed matrix row-wise add test...\n");
  1068.  
  1069. //
  1070. // Test row-wise matrix subtraction
  1071. //
  1072.  
  1073. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1074. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1075. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1076.  
  1077. a = [1,1,1] + [1,1,1]*1i;
  1078. 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;
  1079. 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;
  1080. if (!all (all (a .- b == -c))) { error (); }
  1081. if (!all (all (b .- a ==  c))) { error (); }
  1082.  
  1083. printf("\t\tpassed matrix row-wise subtraction test...\n");
  1084.  
  1085. //
  1086. // Test col-wise matrix addition
  1087. //
  1088.  
  1089. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1090. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1091. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1092.  
  1093. a = [1;1;1;1] + [1;1;1;1]*1i;
  1094. 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;
  1095. 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;
  1096. if (!all (all (a .+ b == c))) { error (); }
  1097. if (!all (all (b .+ a == c))) { error (); }
  1098.  
  1099. printf("\t\tpassed matrix col-wise add test...\n");
  1100.  
  1101. //
  1102. // Test col-wise matrix subtraction
  1103. //
  1104.  
  1105. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1106. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1107. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1108.  
  1109. a = [1;1;1;1] + [1;1;1;1]*1i;
  1110. 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;
  1111. 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;
  1112. if (!all (all (a .- b == -c))) { error (); }
  1113. if (!all (all (b .- a ==  c))) { error (); }
  1114.  
  1115. printf("\t\tpassed matrix col-wise subtraction test...\n");
  1116.  
  1117. a = [1,2,3];
  1118. b = [1,2,3;4,5,6;7,8,9];
  1119. c = [1,4,9;4,10,18;7,16,27];
  1120.  
  1121. if (!all (all (a .* b == c))) { error (); }
  1122. if (!all (all (b .* a == c))) { error (); }
  1123.  
  1124. za = a + rand (size (a))*1j;
  1125. zb = b + rand (size (b))*1j;
  1126.  
  1127. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  1128. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  1129. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  1130. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  1131. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  1132. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  1133.  
  1134. printf("\t\tpassed matrix row-wise multiplication test...\n");
  1135.  
  1136. a = [1,2,3];
  1137. b = [1,2,3;4,6,6;7,8,9];
  1138. c = [1,1,1;4,3,2;7,4,3];
  1139.  
  1140. if (!all (all (b ./ a == c))) { error (); }
  1141. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  1142. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  1143.  
  1144. za = a + rand (size (a))*1j;
  1145. zb = b + rand (size (b))*1j;
  1146.  
  1147. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  1148. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  1149. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  1150. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  1151. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  1152. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  1153.  
  1154. printf("\t\tpassed matrix row-wise division test...\n");
  1155.  
  1156. a = [1;2;3];
  1157. b = [1,2,3;4,5,6;7,8,9];
  1158.  
  1159. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  1160. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  1161.  
  1162. za = a + rand (size (a))*1j;
  1163. zb = b + rand (size (b))*1j;
  1164.  
  1165. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  1166. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  1167. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  1168. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  1169. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  1170. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  1171.  
  1172. printf("\t\tpassed matrix column-wise multiplication test...\n");
  1173.  
  1174. a = [1;2;3];
  1175. b = [1,2,3;4,6,6;7,8,9];
  1176.  
  1177. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  1178. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  1179.  
  1180. za = a + rand (size(a))*1j;
  1181. zb = b + rand (size(b))*1j;
  1182.  
  1183. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  1184. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  1185. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  1186. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  1187. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  1188. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  1189.  
  1190. printf("\t\tpassed matrix column-wise division test...\n");
  1191.  
  1192.  
  1193. //
  1194. //--------------------- Test COMPLEX MATRICES --------------------------
  1195. //
  1196. //  Automatic creation, extension
  1197. //
  1198. printf("\t\tcomplex-matrices\n");
  1199. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  1200. a=[1,2,3;4,5,6];
  1201. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  1202. //
  1203. a = m3;
  1204. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1205. if(any (any (a-a != zeros(2,2)))) { error(); }
  1206. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
  1207. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  1208. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  1209. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1210. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  1211. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  1212. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  1213. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  1214. //
  1215. // The following test may not work on some machines
  1216. //
  1217. if(any (any (a./a != [1,1;1,1]))) { 
  1218.   printf("\t\t***complex division inaccuracy, check manually***\n");
  1219. }
  1220.  
  1221. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  1222. //  
  1223. //--------------------- Test NULL MATRICES -------------------------
  1224. //
  1225. printf("\t\tnull-matrices\n");
  1226. // Create a NULL matrix
  1227. mnull = [];
  1228. // Test it with SCALARS
  1229. if( any([1,mnull] != 1)) {
  1230.   error();
  1231. }
  1232. if( any([mnull,1] != 1)) {
  1233.   error();
  1234. }
  1235. // Test with MATRIX construction
  1236. m = [1,2;3,4;5,6];
  1237. if( any([mnull;1] != [1])) {
  1238.   error();
  1239. }
  1240. if( any([1;mnull] != [1])) {
  1241.   error();
  1242. }
  1243. if( any([mnull;1,2,3] != [1,2,3])) {
  1244.   error();
  1245. }
  1246. if( any([1,2,3;mnull] != [1,2,3])) {
  1247.   error();
  1248. }
  1249. if(any (any ([mnull,m] != m))) {
  1250.   error();
  1251. }
  1252. if(any (any ([m,mnull] != m))) {
  1253.   error();
  1254. }
  1255. if(any (any ([mnull;m] != m))) {
  1256.   error();
  1257. }
  1258. if(any (any ([m;mnull] != m))) {
  1259.   error();
  1260.  
  1261. mnull = matrix();
  1262. mnull[1] = [1];
  1263. }
  1264.  
  1265. //--------------------- Test Matrix Multiply --------------------------
  1266.  
  1267. i = sqrt(-1);
  1268. a = [1,2,3;4,5,6;7,8,9];
  1269. b = [4,5,6;7,8,9;10,11,12];
  1270. c = [ 48,  54,  60 ;
  1271.      111, 126, 141 ;
  1272.      174, 198, 222 ];
  1273.  
  1274. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1275.  
  1276. az = a + b*i;
  1277. bz = b + a*i;
  1278.  
  1279. cz = [-18+141*i , -27+162*i , -36+183*i ;
  1280.         9+240*i ,   0+279*i ,  -9+318*i ;
  1281.        36+339*i ,  27+396*i ,  18+453*i ];
  1282.  
  1283. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  1284.        111+66*i , 126+81*i  , 141+96*i ;
  1285.        174+102*i, 198+126*i , 222+150*i ];
  1286.  
  1287. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  1288.         111+174*i , 126+198*i , 141+222*i ;
  1289.         174+237*i , 198+270*i , 222+303*i ];
  1290.  
  1291. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1292. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1293. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1294.  
  1295. a = [a,a];
  1296. b = [b;b];
  1297. c = [  96 , 108 , 120 ;
  1298.       222 , 252 , 282 ;
  1299.       348 , 396 , 444 ];
  1300.  
  1301. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1302.  
  1303. az = [az,az];
  1304. bz = [bz;bz];
  1305.  
  1306. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  1307.          18+480*i ,    0+558*i ,  -18+636*i ;
  1308.          72+678*i ,   54+792*i ,   36+906*i ];
  1309.  
  1310. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  1311.         222+132*i , 252+162*i , 282+192*i ;
  1312.         348+204*i , 396+252*i , 444+300*i ];
  1313.  
  1314. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  1315.          222+348*i , 252+396*i , 282+444*i ;
  1316.          348+474*i , 396+540*i , 444+606*i ];
  1317.  
  1318. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1319. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1320. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1321.  
  1322. printf("\t\tpassed matrix multiply test...\n");
  1323.  
  1324. //--------------------- Test STRING MATRICES --------------------------
  1325. //
  1326. printf("\t\tstring-matrices\n");
  1327. sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
  1328. if(sm[1] != "s1") { error(); }
  1329. if( sm[1;3] != "sm3" ) { error(); }
  1330. if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
  1331. if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
  1332. {
  1333.   error();
  1334. }
  1335. if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
  1336.  
  1337. //  Test string-matrix add functionality
  1338.  
  1339. sm = sm[1,2;1,2];
  1340. if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
  1341.  
  1342. if ("c"+["1"] != "c1") { error (); }
  1343. if (["c"]+"1" != "c1") { error (); }
  1344.  
  1345. printf("\tpassed matrix test...\n");
  1346.  
  1347. //
  1348. //---------------------------- List Tests --------------------------
  1349. //
  1350. //  List creation
  1351. listest = << << 11; 12 >>; << 21; 22>> >>;
  1352. if( listest.[1].[2] != 12 ) { error(); }
  1353. if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
  1354. mlist.[0] = m;
  1355. if(any(any(mlist.[0] != m))) { error(); }
  1356.  
  1357. // Test list functions...
  1358.  
  1359. listtest.fun = function ( a ) { return sin(a); };
  1360. if (!(listtest.fun(0.5) == sin(0.5))) { error() ; }
  1361.  
  1362. // Test open-list assignment...
  1363. </v;d/> = eig(rand(3,3));
  1364.  
  1365. printf("\tpassed list test...\n");
  1366.  
  1367. //
  1368. // Reset random number generator seed...
  1369. //
  1370.  
  1371. rand("default");
  1372.  
  1373. //
  1374. //-------------------------- Test printf () --------------------------
  1375. //
  1376.  
  1377. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
  1378. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1379. printf("\tpassed sprintf test...\n");
  1380.  
  1381. //
  1382. //------------------------- Test strtod()  ----------------------------
  1383. //
  1384.  
  1385. if (123.456 != strtod ("123.456")) { error (); }
  1386. if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  1387.   { error (); }
  1388. printf("\tpassed strtod test...\n");
  1389.  
  1390. //
  1391. //------------------------- Test getline()  ---------------------------
  1392. //
  1393. //
  1394.  
  1395. close( "test.getline" );
  1396.  
  1397. x = getline( "test.getline" );
  1398. if ( x.[1] !=  123.456 ) { error(); }
  1399. if ( x.[2] != -123.456 ) { error(); }
  1400. if ( x.[3] !=  123.456 ) { error(); }
  1401.  
  1402. x = getline( "test.getline" );
  1403. if ( x.[1] !=  .123 ) { error(); }
  1404. if ( x.[2] != -.123 ) { error(); }
  1405. if ( x.[3] !=  .123 ) { error(); }
  1406.  
  1407. x = getline( "test.getline" );
  1408. if ( x.[1] !=  123 ) { error(); }
  1409. if ( x.[2] != -123 ) { error(); }
  1410. if ( x.[3] !=  123 ) { error(); }
  1411.  
  1412. x = getline( "test.getline" );
  1413. if ( x.[1] !=  1e6 ) { error(); }
  1414. if ( x.[2] != -1e6 ) { error(); }
  1415. if ( x.[3] !=  1e6 ) { error(); }
  1416.  
  1417. x = getline( "test.getline" );
  1418. if ( x.[1] !=  1e5 ) { error(); }
  1419. if ( x.[2] != -1e5 ) { error(); }
  1420. if ( x.[3] !=  1e5 ) { error(); }
  1421.  
  1422. x = getline( "test.getline" );
  1423. if ( x.[1] !=  123.456e3 ) { error(); }
  1424. if ( x.[2] != -123.456e3 ) { error(); }
  1425. if ( x.[3] !=  123.456e3 ) { error(); }
  1426.  
  1427. x = getline( "test.getline" );
  1428. if ( x.[1] !=  123.456e3 ) { error(); }
  1429. if ( x.[2] != -123.456e3 ) { error(); }
  1430. if ( x.[3] !=  123.456e3 ) { error(); }
  1431.  
  1432. x = getline( "test.getline" );
  1433. if ( x.[1] !=  123.456e-3 ) { error(); }
  1434. if ( x.[2] != -123.456e-3 ) { error(); }
  1435. if ( x.[3] !=  123.456e-3 ) { error(); }
  1436.  
  1437. x = getline( "test.getline" );
  1438. if ( x.[1] !=  .123e3 ) { error(); }
  1439. if ( x.[2] != -.123e3 ) { error(); }
  1440. if ( x.[3] !=  .123e3 ) { error(); }
  1441.  
  1442. x = getline( "test.getline" );
  1443. if ( x.[1] !=  123e3 ) { error(); }
  1444. if ( x.[2] != -123e3 ) { error(); }
  1445. if ( x.[3] !=  123e3 ) { error(); }
  1446.  
  1447. x = getline( "test.getline" );
  1448. if ( x.[1] != "123abc" ) { error(); }
  1449. if ( x.[2] != "abc123" ) { error(); }
  1450. if ( x.[3] != "123.abc" ) { error(); }
  1451.  
  1452. x = getline( "test.getline" );
  1453. if ( x.[1] != "quoted string" ) { error(); }
  1454. if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
  1455.  
  1456. x = getline( "test.getline" );
  1457. if ( x.[1] != "quoted string" ) { error(); }
  1458. if ( x.[2] !=  1.23e3 ) { error(); }
  1459. if ( x.[3] !=  100 ) { error(); }
  1460. if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
  1461. if ( x.[5] !=  200 ) { error(); }
  1462.  
  1463. printf("\tpassed getline() test...\n");
  1464.  
  1465. //
  1466. //---------------------- Test readb()/writeb()  --------------------
  1467. //
  1468.  
  1469. a = rand (5,5);
  1470. z = rand(3,3) + rand(3,3)*1j;
  1471. strm = what()[1:5;1:5];
  1472. pi = 4*atan(1.0);
  1473. sc = 2*pi;
  1474. str = "this is a sample string\ttab";
  1475. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1476.  
  1477. writeb ("jnk.rb", a, z, strm, sc, str, l);
  1478.  
  1479. # Set aside matrices for later tests
  1480.  
  1481. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1482.  
  1483. clear (a, z, strm, sc, str, l);
  1484.  
  1485. close ("jnk.rb");
  1486.  
  1487. readb ("jnk.rb");
  1488.  
  1489. #
  1490. # Now do checks
  1491. #
  1492.  
  1493. if (!all (all (a == check.a))) { error (); }
  1494. if (!all (all (z == check.z))) { error (); }
  1495. if (!all (all (strm == check.strm))) { error (); }
  1496. if (sc != check.sc) { error (); }
  1497. if (str != check.str) { error (); }
  1498.  
  1499. if (length (l) != 5) { error (); }
  1500. if (!all (all (l.a == check.a))) { error (); }
  1501. if (!all (all (l.z == check.z))) { error (); }
  1502. if (!all (all (l.strm == check.strm))) { error (); }
  1503. if (l.sc != check.sc) { error (); }
  1504. if (l.str != check.str) { error (); }
  1505.  
  1506.  
  1507. printf("\tpassed binary I/O test...\n");
  1508.  
  1509. //
  1510. //---------------------- Test read()/write()  --------------------
  1511. //
  1512.  
  1513. a = rand (5,5);
  1514. z = rand(3,3) + rand(3,3)*1j;
  1515. strm = what()[1:5;1:5];
  1516. pi = 4*atan(1.0);
  1517. sc = 2*pi;
  1518. str = "this is a sample string\ttab";
  1519. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1520.  
  1521. write ("jnk.ra", a, z, strm, sc, str, l);
  1522.  
  1523. # Set aside matrices for later tests
  1524.  
  1525. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1526.  
  1527. clear (a, z, strm, sc, str, l);
  1528.  
  1529. close ("jnk.ra");
  1530.  
  1531. read ("jnk.ra");
  1532.  
  1533. #
  1534. # Now do checks
  1535. #
  1536.  
  1537. if (a.nr != 5 || a.nc != 5) { error (); }
  1538. if (z.nr != 3 || z.nc != 3) { error (); }
  1539. if (strm.nr != 5 || strm.nc != 5) { error (); }
  1540. if (str != check.str) { error (); }
  1541.  
  1542. if (length (l) != 5) { error (); }
  1543. if (l.a.nr != 5 || l.a.nc != 5) { error (); }
  1544. if (l.z.nr != 3 || l.z.nc != 3) { error (); }
  1545. if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
  1546. if (l.str != check.str) { error (); }
  1547.  
  1548. printf("\tpassed ascii I/O test...\n");
  1549.  
  1550. //
  1551. //-------------------------- Test eval () ------------------------------
  1552. //
  1553.  
  1554. if (1 + 2 != eval("1 + 2")) { error ("eval() error"); }
  1555. x = function (s, a) { return eval(s); };
  1556. str = "yy = 2 + x(\"2*a\",3)";
  1557. z = eval(str);
  1558. if (z != yy) { error ("eval() error"); }
  1559. printf("\tpassed eval test...\n");
  1560.  
  1561.  
  1562. //-------------------------- ------------------------------
  1563.  
  1564. printf ("\ttest builtin function correctness...\n");
  1565.  
  1566. //-------------------------- ------------------------------
  1567.  
  1568. //
  1569. //-------------------------- Test abs () ------------------------------
  1570. //
  1571.  
  1572. A = rand(5,5);
  1573. T = ( A == abs (A));
  1574. if (!all (all (A))) { error ("abs() incorrect"); }
  1575. printf("\tabs()");
  1576.  
  1577.  
  1578. //
  1579. //-------------------------- Test max () ------------------------------
  1580. //
  1581.  
  1582. A = [1,10,100;2,20,200;1,2,3];
  1583. B = A./2;
  1584. ZA = A + rand (3,3)*A*1i;
  1585. ZB = B + rand (3,3)*B*1i;
  1586. if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
  1587. if (max (max(A)) != 200) { error ("max() incorrect"); }
  1588. if (any (any (max (A, B) != A))) { error (); }
  1589. if (any (any (max (B, A) != max (A, B)))) { error (); }
  1590. if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
  1591. if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
  1592. if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
  1593. printf("\tmax()");
  1594.  
  1595. //
  1596. //-------------------------- Test min () ------------------------------
  1597. //
  1598.  
  1599. if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
  1600. if (min (min(A)) != 1) { error ("min() incorrect"); }
  1601. if (any (any (min (A, B) != B))) { error (); }
  1602. if (any (any (min (B, A) != min (A, B)))) { error (); }
  1603. if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
  1604. if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
  1605. if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
  1606. printf("\tmin()");
  1607.  
  1608. //
  1609. //-------------------------- Test maxi () -----------------------------
  1610. //
  1611.  
  1612. if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
  1613. if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
  1614. printf("\tmaxi()");
  1615.  
  1616. //
  1617. //-------------------------- Test mini () -----------------------------
  1618. //
  1619.  
  1620. if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
  1621. if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
  1622. printf("\tmini()");
  1623.  
  1624. //
  1625. //-------------------------- Test floor () ----------------------------
  1626. //
  1627.  
  1628. if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
  1629. if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1630.   error ("floor output incorrect");
  1631. }
  1632. printf("\tfloor()");
  1633.  
  1634. //
  1635. //-------------------------- Test ceil () 0----------------------------
  1636. //
  1637.  
  1638. if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
  1639. if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  1640.   error ("ceil output incorrect");
  1641. }
  1642. printf("\tceil()\n");
  1643.  
  1644. //
  1645. //-------------------------- Test round () ----------------------------
  1646. //
  1647.  
  1648. if (round (1.8) != 2) { error ("round() output incorrect"); }
  1649. if (round (1.4) != 1) { error ("round() output incorrect"); }
  1650. if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  1651.   error ("round output incorrect");
  1652. }
  1653. printf("\tround()");
  1654.  
  1655. //
  1656. //-------------------------- Test sum () ------------------------------
  1657. //
  1658.  
  1659. S = [1:4; 4:7; 8:11];
  1660. if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
  1661. if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
  1662. if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
  1663. printf("\tsum()");
  1664.  
  1665. //
  1666. //-------------------------- Test int () ------------------------------
  1667. //
  1668.  
  1669. if (int (1.9999) != 1) { error ("int() output incorrect"); }
  1670. if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1671.   error ("int() output incorrect");
  1672. }
  1673. printf("\tint()");
  1674.  
  1675. //
  1676. //-------------------------- Test mod () ------------------------------
  1677. //
  1678.  
  1679. if (mod (1,1) != 0) { error ("mod() output incorrect"); }
  1680. if (mod (4,2) != 0) { error ("mod() output incorrect"); }
  1681. if (mod (3,2) != 1) { error ("mod() output incorrect"); }
  1682. if (mod (5,3) != 2) { error ("mod() output incorrect"); }
  1683. printf("\tmod()");
  1684.  
  1685. //
  1686. //-------------------------- Test find () ------------------------------
  1687. //
  1688.  
  1689. if (find ([0,1]) != 2) { error ("find() output incorrect"); }
  1690. if (find ([1,0]) != 1) { error ("find() output incorrect"); }
  1691. if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
  1692. if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
  1693. if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }
  1694.  
  1695. printf("\tfind()\n");
  1696.  
  1697.  
  1698. //-------------------------- ------------------------------
  1699.  
  1700. printf ("\ttest builtin function operation...\n");
  1701.  
  1702. //-------------------------- ------------------------------
  1703.  
  1704. //
  1705. // At least use most of the builtins....
  1706. //
  1707.  
  1708. A = rand(5,5);
  1709. Z = rand(5,5) + rand(5,5)*1j;
  1710. S = what()[3:6;2:4];
  1711.  
  1712. abs (A);
  1713. abs ([A]);
  1714. abs (3.14);
  1715.  
  1716. abs (Z);
  1717. abs ([Z]);
  1718. abs (3.14j);
  1719.  
  1720. printf ("\tabs");
  1721.  
  1722. acos (A);
  1723. acos ([A]);
  1724. acos (3.14);
  1725.  
  1726. acos (Z);
  1727. acos ([Z]);
  1728. acos (3.14j);
  1729.  
  1730. printf ("\tacos");
  1731.  
  1732. all (A);
  1733. all ([A]);
  1734. all (3.14);
  1735.  
  1736. all (Z);
  1737. all ([Z]);
  1738. all (3.14j);
  1739.  
  1740. printf ("\tall");
  1741.  
  1742. any (A);
  1743. any ([A]);
  1744. any (3.14);
  1745.  
  1746. any (Z);
  1747. any ([Z]);
  1748. any (3.14j);
  1749.  
  1750. printf ("\tany");
  1751.  
  1752. asin (A);
  1753. asin ([A]);
  1754. asin (3.14);
  1755.  
  1756. asin (Z);
  1757. asin ([Z]);
  1758. asin (3.14j);
  1759.  
  1760. printf ("\tasin");
  1761.  
  1762. atan (A);
  1763. atan ([A]);
  1764. atan (3.14);
  1765.  
  1766. atan (Z);
  1767. atan ([Z]);
  1768. atan (3.14j);
  1769.  
  1770. printf ("\tatan");
  1771.  
  1772. atan2 (A,A);
  1773. atan2 ([A],[A]);
  1774. atan2 (3.14,3.14);
  1775.  
  1776. #atan2 (Z,Z);
  1777. #atan2 ([Z],[Z]);
  1778. #atan2 (3.14j,3.14j);
  1779.  
  1780. printf ("\tatan2");
  1781. printf ("\n");
  1782.  
  1783. cd (":");
  1784. printf ("\tcd");
  1785.  
  1786. ceil (A);
  1787. ceil ([A]);
  1788. ceil (3.14);
  1789.  
  1790. ceil (Z);
  1791. ceil ([Z]);
  1792. ceil (3.14j);
  1793.  
  1794. printf ("\tceil");
  1795.  
  1796. class (A);
  1797. class ([A]);
  1798. class (3.14);
  1799. class ("string");
  1800. class (S);
  1801. class ([S]);
  1802. class (<< rand(3,3); rand(4,4) >>);
  1803.  
  1804. class (Z);
  1805. class ([Z]);
  1806. class (3.14j);
  1807.  
  1808. printf ("\tclass");
  1809.  
  1810. clear (rand(3,3));
  1811. clear ([rand(3,3)]);
  1812.  
  1813. printf ("\tclear");
  1814.  
  1815. conj (A);
  1816. conj ([A]);
  1817. conj (3.14);
  1818.  
  1819. conj (Z);
  1820. conj ([Z]);
  1821. conj (3.14j);
  1822.  
  1823. printf ("\tconj");
  1824.  
  1825. cos (A);
  1826. cos ([A]);
  1827. cos (3.14);
  1828.  
  1829. cos (Z);
  1830. cos ([Z]);
  1831. cos (3.14j);
  1832.  
  1833. printf ("\tcos");
  1834.  
  1835. diag (A);
  1836. diag ([A]);
  1837. diag (3.14);
  1838.  
  1839. diag (Z);
  1840. diag ([Z]);
  1841. diag (3.14j);
  1842.  
  1843. printf ("\tdiag");
  1844. printf ("\n");
  1845.  
  1846. exist (A);
  1847. exist (Z);
  1848.  
  1849. printf ("\texist");
  1850.  
  1851. exp (A);
  1852. exp ([A]);
  1853. exp (3.14);
  1854.  
  1855. exp (Z);
  1856. exp ([Z]);
  1857. exp (3.14j);
  1858.  
  1859. printf ("\texp");
  1860.  
  1861. find (A);
  1862. find ([A]);
  1863. find (3.14);
  1864.  
  1865. find (Z);
  1866. find ([Z]);
  1867. find (3.14j);
  1868.  
  1869. printf ("\tfind");
  1870.  
  1871. floor (A);
  1872. floor ([A]);
  1873. floor (3.14);
  1874.  
  1875. floor (Z);
  1876. floor ([Z]);
  1877. floor (3.14j);
  1878.  
  1879. printf ("\tfloor");
  1880.  
  1881. format (1,1);
  1882. format ([2],[3]);
  1883. format ();
  1884. format ([3], 3);
  1885. format ();
  1886.  
  1887. printf ("\tformat");
  1888.  
  1889. imag (A);
  1890. imag ([A]);
  1891. imag (3.14);
  1892.  
  1893. imag (Z);
  1894. imag ([Z]);
  1895. imag (3.14j);
  1896.  
  1897. printf ("\timag");
  1898.  
  1899. inf ();
  1900. printf ("\tinf");
  1901. printf ("\n");
  1902.  
  1903. int (A);
  1904. int ([A]);
  1905. int (3.14);
  1906.  
  1907. int (Z);
  1908. int ([Z]);
  1909. int (3.14j);
  1910.  
  1911. printf ("\tint");
  1912.  
  1913. issymm (A);
  1914. issymm ([A]);
  1915. issymm (3.14);
  1916.  
  1917. issymm (Z);
  1918. issymm ([Z]);
  1919. issymm (3.14j);
  1920.  
  1921. printf ("\tissymm");
  1922.  
  1923. length (A);
  1924. length ([A]);
  1925. length (3.14);
  1926.  
  1927. length (Z);
  1928. length ([Z]);
  1929. length (3.14j);
  1930.  
  1931. length (3.14);
  1932. length (S);
  1933. length (<< rand(2,2); rand(10,10); 3.14 >>);
  1934.  
  1935. printf ("\tlength");
  1936.  
  1937. log (A);
  1938. log ([A]);
  1939. log (3.14);
  1940.  
  1941. log (Z);
  1942. log ([Z]);
  1943. log (3.14j);
  1944.  
  1945. printf ("\tlog");
  1946.  
  1947. log10 (A);
  1948. log10 ([A]);
  1949. log10 (3.14);
  1950.  
  1951. log10 (Z);
  1952. log10 ([Z]);
  1953. log10 (3.14j);
  1954.  
  1955. printf ("\tlog10");
  1956.  
  1957. max (A);
  1958. max ([A]);
  1959. max (3.14);
  1960.  
  1961. max (Z);
  1962. max ([Z]);
  1963. max (3.14j);
  1964.  
  1965. printf ("\tmax");
  1966.  
  1967. maxi (A);
  1968. maxi ([A]);
  1969. maxi (3.14);
  1970.  
  1971. maxi (Z);
  1972. maxi ([Z]);
  1973. maxi (3.14j);
  1974.  
  1975. printf ("\tmaxi");
  1976. printf ("\n");
  1977.  
  1978. members (<< rand(2,2); rand(); rand(3,3) >>);
  1979.  
  1980. printf ("\tmembers");
  1981.  
  1982. min (A);
  1983. min ([A]);
  1984. min (3.14);
  1985.  
  1986. min (Z);
  1987. min ([Z]);
  1988. min (3.14j);
  1989.  
  1990. printf ("\tmin");
  1991.  
  1992. mini (A);
  1993. mini ([A]);
  1994. mini (3.14);
  1995.  
  1996. mini (Z);
  1997. mini ([Z]);
  1998. mini (3.14j);
  1999.  
  2000. printf ("\tmini");
  2001.  
  2002. mod (A,A);
  2003. mod ([A],A);
  2004. mod (3.14,2);
  2005.  
  2006. mod (Z,Z);
  2007. mod ([Z],Z);
  2008. mod (3.14j,2);
  2009.  
  2010. printf ("\tmod");
  2011.  
  2012. nan ();
  2013.  
  2014. ones(3,3);
  2015. ones([4,4]);
  2016.  
  2017. printf ("\tones");
  2018.  
  2019. prod (A);
  2020. prod ([A]);
  2021. prod (3.14);
  2022.  
  2023. prod (Z);
  2024. prod ([Z]);
  2025. prod (3.14j);
  2026.  
  2027. printf ("\tprod");
  2028.  
  2029. real (A);
  2030. real ([A]);
  2031. real (3.14);
  2032.  
  2033. real (Z);
  2034. real ([Z]);
  2035. real (3.14j);
  2036.  
  2037. printf ("\treal");
  2038. printf ("\n");
  2039.  
  2040. reshape (A, A.nr*A.nc, 1);
  2041. reshape (A, [A.nr*A.nc], [1]);
  2042. reshape (Z, Z.nr*Z.nc, 1);
  2043. reshape (Z, [Z.nr*Z.nc], [1]);
  2044.  
  2045. printf ("\treshape");
  2046.  
  2047. round (A);
  2048. round ([A]);
  2049. round (3.14);
  2050.  
  2051. round (Z);
  2052. round ([Z]);
  2053. round (3.14j);
  2054.  
  2055. printf ("\tround");
  2056.  
  2057. show(A);
  2058. show(Z);
  2059. show([A]);
  2060. show([Z]);
  2061.  
  2062. printf ("\tshow");
  2063.  
  2064. sign (A);
  2065. sign ([A]);
  2066. sign (3.14);
  2067.  
  2068. sign (Z);
  2069. sign ([Z]);
  2070. sign (3.14j);
  2071.  
  2072. printf ("\tsign");
  2073.  
  2074. sin (A);
  2075. sin ([A]);
  2076. sin (3.14);
  2077.  
  2078. sin (Z);
  2079. sin ([Z]);
  2080. sin (3.14j);
  2081.  
  2082. printf ("\tsin");
  2083.  
  2084. size (A);
  2085. size ([A]);
  2086. size (3.14);
  2087.  
  2088. size (Z);
  2089. size ([Z]);
  2090. size (3.14j);
  2091.  
  2092. printf ("\tsize");
  2093.  
  2094. sizeof (A);
  2095. sizeof ([A]);
  2096. sizeof (3.14);
  2097.  
  2098. sizeof (Z);
  2099. sizeof ([Z]);
  2100. sizeof (3.14j);
  2101.  
  2102. printf ("\tsizeof");
  2103. printf ("\n");
  2104.  
  2105. sort (A);
  2106. sort ([A]);
  2107. sort (3.14);
  2108.  
  2109. sort (Z);
  2110. sort ([Z]);
  2111. sort (3.14j);
  2112.  
  2113. printf ("\tsort");
  2114.  
  2115. sqrt (A);
  2116. sqrt ([A]);
  2117. sqrt (3.14);
  2118.  
  2119. sqrt (Z);
  2120. sqrt ([Z]);
  2121. sqrt (3.14j);
  2122.  
  2123. printf ("\tsqrt");
  2124.  
  2125. srand ();
  2126. srand ("clock");
  2127. srand (SEED);
  2128.  
  2129. printf ("\tsrand");
  2130.  
  2131. strsplt (S);
  2132. printf ("\tstrsplt");
  2133.  
  2134. strtod ("string");
  2135. strtod (S);
  2136. strtod ("1.2");
  2137. strtod (["1.2", "1e3"]);
  2138.  
  2139. printf ("\tstrtod");
  2140.  
  2141. sum (A);
  2142. sum ([A]);
  2143. sum (3.14);
  2144.  
  2145. sum (Z);
  2146. sum ([Z]);
  2147. sum (3.14j);
  2148.  
  2149. printf ("\tsum");
  2150.  
  2151. tan (A);
  2152. tan ([A]);
  2153. tan (3.14);
  2154.  
  2155. tan (Z);
  2156. tan ([Z]);
  2157. tan (3.14j);
  2158.  
  2159. printf ("\ttan");
  2160. printf ("\n");
  2161.  
  2162. tmpnam ();
  2163. printf ("\ttmpnam");
  2164.  
  2165. type (A);
  2166. type ([A]);
  2167. type (Z);
  2168. type ([Z]);
  2169. type (3.14);
  2170. type (S);
  2171. type ("string");
  2172.  
  2173. printf ("\ttype");
  2174.  
  2175. zeros (3,3);
  2176. zeros ([3], [3]);
  2177. zeros ([3,3]);
  2178.  
  2179. printf ("\tzeros");
  2180.  
  2181. printf ("\n");
  2182.  
  2183. srand(SEED);
  2184. rand("default");
  2185.  
  2186. //
  2187. //-------------------------- Test rand () -----------------------------
  2188. //
  2189. rand ("normal", 5, 1);
  2190. xrand = rand(4000,1);
  2191.  
  2192. mean = function(x)
  2193. {
  2194.   local(m);
  2195.  
  2196.   m = size (x)[1];
  2197.   if( m == 1 ) 
  2198.   { 
  2199.     m = size (x)[2];
  2200.   }
  2201.  
  2202.   return sum( x ) / m;
  2203. };
  2204.  
  2205. std = function(x)
  2206. {
  2207.   local(i, m, s);
  2208.  
  2209.   if(class(x) != "num") { error("std() requires NUMERICAL input"); }
  2210.  
  2211.   m = x.nr;
  2212.   if( m == 1 ) 
  2213.   { 
  2214.     return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  2215.   else
  2216.     for( i in 1:x.nc) {
  2217.       s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
  2218.     }
  2219.     return s;
  2220.   }
  2221. };
  2222.  
  2223. if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1)) 
  2224.   { error ("error in random"); }
  2225. if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  2226.   { error ("error in random"); }
  2227. printf("\tpassed rand test...\n");
  2228.  
  2229. rand("default");
  2230.  
  2231. //
  2232. //-------------------------- Test norm () -----------------------------
  2233. //
  2234.  
  2235. tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
  2236. if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
  2237. if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
  2238. if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
  2239. printf("\tpassed norm test...\n");
  2240.  
  2241. //
  2242. //-------------------------- Test qr () ------------------------------
  2243. //
  2244.  
  2245. a = ohess(4);
  2246. qa = qr (a);
  2247. if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  2248.   { error ("possible qr() problems"); }
  2249.  
  2250. z = ohess (4) + ohess(4)*1i;
  2251. qz = qr (z);
  2252. if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  2253.   { error ("possible qr() problems"); }
  2254.  
  2255. printf("\tpassed qr test...\n");
  2256.  
  2257. //
  2258. //-------------------------- Test schur () ----------------------------
  2259. //
  2260.  
  2261. a = randsvd (10, 10);
  2262. sa = schur (a);
  2263. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  2264.   { error ("possible schur() problems"); }
  2265.  
  2266. z = rand (4,4) + rand(4,4)*1i;
  2267. sz = schur (z);
  2268. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  2269.   { error ("possible schur() problems"); }
  2270.  
  2271. a = randsvd (10, -10);
  2272. sa = schur (a);
  2273. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  2274.   { error ("possible schur() problems"); }
  2275.  
  2276. z = rand (4,4) + rand(4,4)*1i;
  2277. sz = schur (z);
  2278. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  2279.   { error ("possible schur() problems"); }
  2280.  
  2281. printf("\tpassed schur test...\n");
  2282.  
  2283. //
  2284. //-------------------------- Test schord () ----------------------------
  2285. //
  2286.  
  2287. s2a = schord (sa, 2, 4);
  2288. if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
  2289.   error ("possible schord() problems"); 
  2290. }
  2291.  
  2292. s2z = schord (sz, 3, 1);
  2293. if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
  2294.   error ("possible schord() problems"); 
  2295. }
  2296.  
  2297. printf("\tpassed schord test...\n");
  2298.  
  2299. //
  2300. //-------------------------- Test chol () -----------------------------
  2301. //
  2302.  
  2303. c = lehmer(10);
  2304. u = chol (c);
  2305. if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
  2306.   error ("possible chol() problems"); 
  2307. }
  2308.  
  2309. cz = lehmer(10) + lehmer(10)*1j;
  2310. cz = symm(cz);
  2311. uz = chol (cz);
  2312. if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
  2313.   error ("possible chol() problems"); 
  2314. }
  2315.  
  2316. printf("\tpassed chol test...\n");
  2317.  
  2318. //
  2319. //--------------------------- Test inv () --------------------------------
  2320. //
  2321.  
  2322. a = randsvd(10,10);
  2323. b = ones(10,1);
  2324. x = inv(a) * b;
  2325. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2326.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2327.   printf ("\tA*X - B:\n");
  2328.   abs (a*x - b)
  2329.   error ("possible inv() problems\n");
  2330. }
  2331.  
  2332. az = randsvd(10,10) + randsvd(10,10)*1j;
  2333. bz = rand(10,1) + rand(10,1)*1j;
  2334. xz = inv (az)*bz;
  2335. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2336.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2337.   printf ("\tA*X - B:\n");
  2338.   abs (az*xz - bz)
  2339.   error ("possible inv() problems\n");
  2340. }
  2341.  
  2342. printf("\tpassed inv test...\n");
  2343.  
  2344. //
  2345. //-------------------------- Test solve () -----------------------------
  2346. //
  2347.  
  2348. //
  2349. // Real - General case
  2350. //
  2351.  
  2352. a = randsvd(10,10);
  2353. b = ones(10,1);
  2354. x = solve (a,b);
  2355. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2356.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2357.   printf ("\tA*X - B:\n");
  2358.   abs (a*x - b)
  2359.   error ("possible solve() problems\n");
  2360. }
  2361.  
  2362. //
  2363. // Real - Symmetric case
  2364. //
  2365.  
  2366. s = symm (randsvd(10,10));
  2367. b = ones(10,1);
  2368. x = solve (s,b);
  2369. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2370.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2371.   printf ("\tA*X - B:\n");
  2372.   abs (s*x - b)
  2373.   error ("possible solve() problems\n");
  2374. }
  2375.  
  2376. //
  2377. // Complex - General  case
  2378. //
  2379.  
  2380. az = randsvd(10,10) + randsvd(10,10)*1j;
  2381. bz = rand(10,1) + rand(10,1)*1j;
  2382. xz = solve (az,bz);
  2383. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2384.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2385.   printf ("\tA*X - B:\n");
  2386.   abs (az*xz - bz)
  2387.   error ("possible solve() problems\n");
  2388. }
  2389.  
  2390. //
  2391. // Complex - Symmetric case
  2392. //
  2393.  
  2394. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2395. bz = rand(10,1) + rand(10,1)*1j;
  2396. xz = solve (sz,bz);
  2397. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2398.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2399.   printf ("\tA*X - B:\n");
  2400.   abs (sz*xz - bz)
  2401.   error ("possible solve() problems\n");
  2402. }
  2403.  
  2404. printf("\tpassed solve test...\n");
  2405.  
  2406. //
  2407. //-------------------------- Test factor() / backsub() -----------------------------
  2408. //
  2409.  
  2410. //
  2411. // Real - General case
  2412. //
  2413.  
  2414. a = randsvd(10,10);
  2415. b = ones(10,1);
  2416. f = factor (a);
  2417. x = backsub (f,b);
  2418. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2419.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2420.   printf ("\tA*X - B:\n");
  2421.   abs (a*x - b)
  2422.   error ("possible factor/backsub problems\n");
  2423. }
  2424.  
  2425. //
  2426. // Real - Symmetric case
  2427. //
  2428.  
  2429. s = symm (randsvd(10,10));
  2430. f = factor (s);
  2431. x = backsub (f,b);
  2432. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2433.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2434.   printf ("\tA*X - B:\n");
  2435.   abs (s*x - b)
  2436.   error ("possible factor/backsub problems\n");
  2437. }
  2438.  
  2439. s = symm (randsvd(10,10));
  2440. f = factor (s, "s");
  2441. x = backsub (f,b);
  2442. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2443.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2444.   printf ("\tA*X - B:\n");
  2445.   abs (s*x - b)
  2446.   error ("possible factor/backsub problems\n");
  2447. }
  2448.  
  2449. //
  2450. // Complex - General  case
  2451. //
  2452.  
  2453. az = randsvd(10,10) + randsvd(10,10)*1j;
  2454. bz = rand(10,1) + rand(10,1)*1j;
  2455. fz = factor (az);
  2456. xz = backsub (fz,bz);
  2457. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2458.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2459.   printf ("\tA*X - B:\n");
  2460.   abs (az*xz - bz)
  2461.   error ("possible factor/backsub problems\n");
  2462. }
  2463.  
  2464. az = randsvd(10,10) + randsvd(10,10)*1j;
  2465. bz = rand(10,1) + rand(10,1)*1j;
  2466. fz = factor (az, "g");
  2467. xz = backsub (fz,bz);
  2468. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2469.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2470.   printf ("\tA*X - B:\n");
  2471.   abs (az*xz - bz)
  2472.   error ("possible factor/backsub problems\n");
  2473. }
  2474.  
  2475. //
  2476. // Complex - Symmetric case
  2477. //
  2478.  
  2479. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2480. fz = factor(sz);
  2481. xz = backsub (fz,bz);
  2482. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2483.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2484.   printf ("\tA*X - B:\n");
  2485.   abs (sz*xz - bz)
  2486.   error ("possible factor/backsub problems\n");
  2487. }
  2488.  
  2489. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2490. fz = factor(sz, "s");
  2491. xz = backsub (fz,bz);
  2492. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2493.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2494.   printf ("\tA*X - B:\n");
  2495.   abs (sz*xz - bz)
  2496.   error ("possible factor/backsub problems\n");
  2497. }
  2498.  
  2499. printf("\tpassed factor/backsub test...\n");
  2500.  
  2501. //
  2502. //------------------------------ Test lu() ---------------------------------
  2503. //
  2504.  
  2505. static (swap);
  2506.  
  2507. lu = function ( A )
  2508. {
  2509.   local (i, l, u, pvt, x)
  2510.  
  2511.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  2512.  
  2513.   x = factor (A, "g");    // Do the factorization
  2514.  
  2515.   //
  2516.   // Now create l, u, and pvt from lu and pvt.
  2517.   //
  2518.  
  2519.   l = tril (x.lu, -1) + eye (size (x.lu));
  2520.   u = triu (x.lu);
  2521.   pvt = eye (size (x.lu));
  2522.  
  2523.   //
  2524.   // Now re-arange the columns of pvt
  2525.   //
  2526.  
  2527.   for (i in 1:max (size (x.lu)))
  2528.   {
  2529.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  2530.   }
  2531.   return << l = l; u = u; pvt = pvt >>;
  2532. };
  2533.  
  2534. //
  2535. //  In vector V, swap elements I, J
  2536. //
  2537.  
  2538. swap = function ( V, I, J )
  2539. {
  2540.   local (v, tmp);
  2541.   v = V;
  2542.   tmp = v[I];
  2543.   v[I] = v[J];
  2544.   v[J] = tmp;
  2545.   return v;
  2546. };
  2547.  
  2548. a = randsvd(10,10);
  2549. lua = lu (a);
  2550. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  2551.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2552.   printf ("\tA - p*l*u:\n");
  2553.   abs (a - lua.pvt*lua.l*lua.u)
  2554.   error ("possible lu()/factor() problems\n");
  2555. }
  2556.  
  2557. //
  2558. // Real
  2559. az = randsvd(10,10) + randsvd(10,10)*1j;
  2560. luz = lu (az);
  2561. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  2562.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2563.   printf ("\tA - p*l*u:\n");
  2564.   abs (az - luz.ovt*luz.l*luz.u)
  2565.   error ("possible lu()/factor()() problems\n");
  2566. }
  2567.  
  2568. printf("\tpassed lu/factor test...\n");
  2569.  
  2570. //
  2571. //-------------------------- Test svd ()   -----------------------------
  2572. //
  2573.  
  2574. a = randsvd(10,10);
  2575. s = svd (a);
  2576. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
  2577. {
  2578.   error ("possible svd() problems"); 
  2579. }
  2580.  
  2581. z = randsvd(10,10) + rand(10,10)*1j;
  2582. sz = svd (z);
  2583. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
  2584.   error ("possible svd() problems"); 
  2585. }
  2586.  
  2587. printf("\tpassed svd test...\n");
  2588.  
  2589. //
  2590. //-------------------------- Test hess ()  -----------------------------
  2591. //
  2592.  
  2593. a = randsvd(10,10);
  2594. h = hess (a);
  2595. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  2596.   error ("possible hess() problems");
  2597. }
  2598.  
  2599. z = randsvd(10,10) + randsvd(10,10)*1j;
  2600. hz = hess (z);
  2601. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  2602.   error ("possible hess() problems"); 
  2603. }
  2604.  
  2605. printf("\tpassed hess test...\n");
  2606.  
  2607. //
  2608. //-------------------------- Test lyap () ------------------------------
  2609. //
  2610.  
  2611. lyap = function ( A, B, C )
  2612. {
  2613.   local (A, B, C)
  2614.  
  2615.   if (!exist (B)) 
  2616.   { 
  2617.     B = A';    // Solve the special form: A*X + X*A' = -C
  2618.   }
  2619.  
  2620.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  2621.     error ("Dimensions do not agree.");
  2622.   }
  2623.  
  2624.   //
  2625.   // Schur decomposition on A and B
  2626.   //
  2627.  
  2628.   sa = schur (A);
  2629.   sb = schur (B);
  2630.  
  2631.   //
  2632.   // transform C
  2633.   //
  2634.  
  2635.   tc = sa.z' * C * sb.z;
  2636.  
  2637.   X = sylv (sa.t, sb.t, tc);
  2638.  
  2639.   //
  2640.   // Undo the transformation
  2641.   //
  2642.  
  2643.   X = sa.z * X * sb.z';
  2644.  
  2645.   return X;
  2646. };
  2647.  
  2648. a = randsvd (10,10);
  2649. b = rand (10,10);
  2650. c = rand (10,10);
  2651.  
  2652. x = lyap (a, b, c);
  2653. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  2654.   error ("possible problems with lyap() or sylv()"); 
  2655. }
  2656.  
  2657. printf("\tpassed lyap test...\n");
  2658.  
  2659. //
  2660. //-------------------------- Test eig () ------------------------------
  2661. //
  2662.  
  2663. trace = function(m) 
  2664. {
  2665.   local(i, tr);
  2666.  
  2667.   if(m.class != "num") { 
  2668.     error("must provide NUMERICAL input to trace()");
  2669.   }
  2670.  
  2671.   tr = 0;
  2672.   for(i in 1:min( [m.nr, m.nc] )) {
  2673.     tr = tr + m[i;i];
  2674.   }
  2675.  
  2676.   return tr;
  2677. };
  2678.  
  2679. eye = function( m , n ) 
  2680. {
  2681.   local(i, N, new);
  2682.  
  2683.   if (!exist (n))
  2684.   {
  2685.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  2686.     new = zeros (m[1], m[2]);
  2687.     N = min ([m[1], m[2]]);
  2688.   else
  2689.     if (class (m) == "string" || class (n) == "string") {
  2690.       error ("eye(), string arguments not allowed");
  2691.     }
  2692.     if (max (size (m)) == 1 && max (size (n)) == 1)
  2693.     {
  2694.       new = zeros (m[1], n[1]);
  2695.       N = min ([m[1], n[1]]);
  2696.     else
  2697.       error ("matrix arguments to eye() must be 1x1");
  2698.     }
  2699.   }
  2700.   for(i in 1:N)
  2701.   {
  2702.     new[i;i] = 1.0;
  2703.   }
  2704.   return new;
  2705. };
  2706.  
  2707. //
  2708. // Standard eigenvalue problem
  2709. //
  2710.  
  2711. a = randsvd(10,10);
  2712. ta = trace (a);
  2713. sa = symm (a);
  2714. tsa = trace (sa);
  2715.  
  2716. z = randsvd(10,10) + randsvd(10,10)*1i;
  2717. tz = trace (z);
  2718. sz = symm (z);
  2719. tsz = trace (sz);
  2720.  
  2721. tol = 1.e-6;
  2722.  
  2723. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  2724. {
  2725.   error ("error in eig");
  2726. }
  2727. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  2728. {
  2729.   error ("error in eig");
  2730. }
  2731. if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
  2732. {
  2733.   error ("error in eig");
  2734. }
  2735. if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
  2736. {
  2737.   error ("error in eig");
  2738. }
  2739.  
  2740. //
  2741. // Generalized eigenvalue problem
  2742. //
  2743.  
  2744. b = randsvd(10,10);
  2745. sb = symm (b) + eye(size(b))*3;
  2746. tb = trace (b);
  2747. tsb = trace (sb);
  2748.  
  2749. zb = randsvd(10,10) + randsvd(10,10)*1i;
  2750. szb = symm (zb) + eye(size(zb))*3;
  2751. tzb = trace (zb);
  2752. tszb = trace (szb);
  2753.  
  2754. eig(a,b);    // not sure of a good way to check these yet
  2755. eigs(sa,sb);
  2756. eigs(sa,sb);
  2757.  
  2758. eig(z, zb);
  2759. eigs(sz, szb);
  2760. eigs(sz, szb);
  2761.  
  2762. printf("\tpassed eig test...\n");
  2763.  
  2764. //
  2765. //-------------------------- Test fft () -----------------------------
  2766. //
  2767.  
  2768. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  2769. printf("\tpassed fft test...\n");
  2770.  
  2771. //
  2772. //------------------------- Fibonacci Test -------------------------
  2773. //
  2774. //  Calculate Fibonacci numbers
  2775. //
  2776.  
  2777. i=1; 
  2778. while ( i < 2 ) { 
  2779.   i=i+1;
  2780.   a=0; b=1;
  2781.   while ( b < 10000 ) {
  2782.     c = b;
  2783.     b = a+b;
  2784.     a = c;
  2785.   }
  2786. }
  2787. if ( b != 10946 ) {
  2788.   error("failed fibonacci test");
  2789. else
  2790.   printf("\tpassed fibonacci test...\n");
  2791. }
  2792.  
  2793. //
  2794. //------------------------- Factorial Test -------------------------
  2795. //
  2796.  
  2797. fac = function(a) 
  2798. {
  2799.   if(a <= 1) 
  2800.   {
  2801.     return 1
  2802.   else
  2803.     return a*$self(a-1)
  2804.   }
  2805. };
  2806.  
  2807. if(fac(10) != 3628800)
  2808.   error(); else printf("\tpassed factorial test...\n"); 
  2809. }
  2810.  
  2811. //
  2812. //--------------------------- ACK Test ----------------------------
  2813. //
  2814.  
  2815. ack = function(a, b) 
  2816. {
  2817.   if(a == 0) { return b + 1; }
  2818.   if(b == 0) { return $self(a-1, 1); }
  2819.   return $self(a-1, $self(a, b-1));
  2820. };
  2821.  
  2822. if(ack(2,2) != 7) 
  2823. {
  2824.   error("error in ack() test");
  2825.   else
  2826.   printf("\tpassed ACK test...\n");
  2827. }
  2828.  
  2829. //
  2830. //------------------------- Prime Test -----------------------------
  2831. //
  2832. // An example that finds all primes less than limit
  2833. //
  2834.  
  2835. primes = function (limit) 
  2836. {
  2837.   local(prime, cnt, i, j, k);
  2838.   
  2839.   i = 1; cnt = 0;
  2840.   for(k in 2:limit) 
  2841.   {
  2842.     j = 2;
  2843.     while(mod(k,j) != 0) 
  2844.     {
  2845.       j++;
  2846.     }
  2847.     if(j == k)             // Found prime
  2848.     {
  2849.       cnt++;
  2850.       prime[i;1] = k;
  2851.       i++;
  2852.     }
  2853.   }
  2854.   return prime;
  2855. };
  2856.  
  2857. if(max(size(primes(100))) != 25) 
  2858. {  
  2859.   error("error in prime test");
  2860.   else
  2861.   printf("\tpassed prime test...\n");
  2862. }
  2863.  
  2864. //
  2865. //--------------------------- Fib Min Test -----------------------------
  2866. //
  2867. //    fibmin() will minimize an arbitrary function 
  2868. //    in 1D using Fibonacci search
  2869.  
  2870. f065 = function(x)
  2871. {
  2872.   return (x - 0.65) * (x - 0.65);
  2873. };
  2874.  
  2875. fib = function(x)
  2876. {
  2877.   local (n, a, b);
  2878.   
  2879.   a = 1;
  2880.   b = 1;
  2881.   if (x >= 2)
  2882.   {
  2883.     n = x - 1;
  2884.     for (n in n:1:-1)
  2885.     {
  2886.       c = b;
  2887.       b = a + b;
  2888.       a = c;
  2889.       n = n - 1;
  2890.     }
  2891.   }
  2892.   return b;
  2893. };
  2894.  
  2895. //  Minimize a 1D function using Fibonacci search
  2896. //  f = function to minimize
  2897. //  xlo = lower bound
  2898. //  xhi = upper bound
  2899. //  n = number of iterations (the bigger the more accurate)
  2900.  
  2901. fibmin = function(f, xlo, xhi, n) 
  2902. {
  2903.   local(a, b, x, y, ex, ey, k, lo, hi);
  2904.   
  2905.   lo = xlo;
  2906.   hi = xhi;
  2907.   k = n;
  2908.   for (k in k:2:-1)
  2909.   {
  2910.     a = fib(k - 2) / fib(k);
  2911.     b = fib(k - 1) / fib(k);
  2912.     x = lo + (hi - lo) * a;
  2913.     y = lo + (hi - lo) * b;
  2914.     ex = f(x);
  2915.     ey = f(y);
  2916.     if (ex >= ey)
  2917.     {
  2918.       lo = x;
  2919.       else
  2920.       hi = y;
  2921.     }
  2922.     //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  2923.   }
  2924.   return (lo + hi) / 2;
  2925. };
  2926.  
  2927. //
  2928. // Simple example using above function to mimize f065. Answer is 0.65
  2929. //
  2930.  
  2931. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  2932. if (abs(x - 0.65) > 1e-6)
  2933. {
  2934.   printf("x = %f\n", x);
  2935.   error("failed fibmin test");
  2936. }
  2937.  
  2938. printf("\tpassed fibmin test...\n");
  2939.  
  2940. //
  2941. //--------------------- Nasty Function Test ------------------------
  2942. //
  2943.  
  2944. printf("\tStarting Nasty Function Test...");
  2945. printf("\tthis will take awhile\n");
  2946. check = function( a, b, c, d, e, f, g, h ) {
  2947.   if ( a+b+c+d == e+f+g+h && ...
  2948.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  2949.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
  2950.     return 1;
  2951.   else
  2952.     return 0;
  2953.   }
  2954. };
  2955.  
  2956. for(a in 8:10) {
  2957.   for(b in 7:(a-1)) {
  2958.     for(c in 6:(b-1)) {
  2959.       for(d in 5:(c-1)) {
  2960.         for(e in 4:(d-1)) {
  2961.           for(f in 3:(e-1)) {
  2962.             for(g in 2:(f-1)) {
  2963.               for(h in 1:(g-1)) {                  
  2964.               if(check( a, b, c, d,  e, f, g, h ) || ...
  2965.                      check( a, e, c, d,  b, f, g, h ) || ...
  2966.                      check( a, f, c, d,  e, b, g, h ) || ...
  2967.                      check( a, g, c, d,  e, f, b, h ) || ...
  2968.                      check( a, h, c, d,  e, f, g, b ) || ...
  2969.                      check( a, b, e, d,  c, f, g, h ) || ...
  2970.                      check( a, b, f, d,  e, c, g, h ) || ...
  2971.                      check( a, b, g, d,  e, f, c, h ) || ...
  2972.                      check( a, b, h, d,  e, f, g, c ) || ...
  2973.                      check( a, b, c, e,  d, f, g, h ) || ...
  2974.                      check( a, b, c, f,  e, d, g, h ) || ...
  2975.                      check( a, b, c, g,  e, f, d, h ) || ...
  2976.                      check( a, b, c, h,  e, f, g, d ) || ...
  2977.                      check( a, e, f, d,  b, c, g, h ) || ...
  2978.                      check( a, e, g, d,  b, f, c, h ) || ...
  2979.                      check( a, e, h, d,  b, f, g, c ) || ...
  2980.                      check( a, f, g, d,  e, b, c, h ) || ...
  2981.                      check( a, f, h, d,  e, b, g, c ) || ...
  2982.                      check( a, g, h, d,  e, f, b, c ) || ...
  2983.                      check( a, b, e, f,  c, d, g, h ) || ...
  2984.                      check( a, b, e, g,  c, f, d, h ) || ...
  2985.                      check( a, b, e, h,  c, f, g, d ) || ...
  2986.                      check( a, b, f, g,  e, c, d, h ) || ...
  2987.                      check( a, b, f, h,  e, c, g, d ) || ...
  2988.                      check( a, b, g, h,  e, f, c, d ) || ...
  2989.                      check( a, e, f, g,  e, f, g, h ) || ...
  2990.                      check( a, e, f, h,  e, f, g, h ) || ...
  2991.                      check( a, e, g, h,  e, f, g, h ) || ...
  2992.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  2993.               }
  2994.             }
  2995.           }
  2996.         }
  2997.       }
  2998.     }
  2999.   }
  3000. }
  3001.  
  3002. if(1) {  // figure out the value of cnt, and check!
  3003.   printf("\tpassed nasty function test...\n");
  3004. else
  3005.   error();
  3006. }
  3007.  
  3008. //
  3009. //------------------ Test More Advanced Functions --------------------
  3010. //
  3011.  
  3012. printf( "\tStarting the lqr/ode test..." );
  3013. printf( "\tthis will take awhile\n" );
  3014.  
  3015. lqr = function( a, b, q, r ) 
  3016. {
  3017.   local( k, s,...
  3018.          m, n, mb, nb, mq, nq,...
  3019.          e, v, d );
  3020.  
  3021.   m = size(a)[1]; n = size(a)[2];
  3022.   mb = size(b)[1]; nb = size(b)[2];
  3023.   mq = size(q)[1]; nq = size(q)[2];
  3024.   
  3025.   if ( m != mq || n != nq ) 
  3026.   {
  3027.     fprintf( "stderr", "A and Q must be the same size.\n" );
  3028.     quit
  3029.   }
  3030.     
  3031.   mr = size(r)[1]; nr = size(r)[2];
  3032.   if ( mr != nr || nb != mr ) 
  3033.   {
  3034.     fprintf( "stderr", "B and R must be consistent.\n" );
  3035.     quit
  3036.   }
  3037.     
  3038.   nn = zeros( m, nb );
  3039.     
  3040.   // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  3041.     
  3042.   e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  3043.   v = e.vec; d = e.val;
  3044.     
  3045.   index = sort( real( d ) ).ind;
  3046.   d = real( d[ index ] );
  3047.   
  3048.   if ( !( d[n] < 0 && d[n+1] > 0 ) ) 
  3049.   {
  3050.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  3051.     quit
  3052.   }
  3053.   
  3054.   chi = v[ 1:n; index[1:n] ];
  3055.   lambda = v[ (n+1):(2*n); index[1:n] ];
  3056.   s = -real(solve(chi',lambda')');
  3057.   k = solve( r, nn'+b'*s );
  3058.   
  3059.   return << k=k; s=s >>;
  3060.   
  3061. };
  3062.  
  3063. // Now run a little test problem.
  3064.  
  3065. k = 1; m = 1; c = .1;
  3066. a = [0     ,1    ,0    , 0;
  3067.     -k/m, -c/m,  k/m,  c/m;
  3068.      0,     0,    0,    1;
  3069.      k/m,  c/m, -k/m, -c/m ];
  3070. b = [ 0; 1/m; 0; 0 ];
  3071. qxx = diag( [0, 0, 100, 0] );
  3072. ruu = [1];
  3073. K = lqr( a, b, qxx, ruu ).k;
  3074.  
  3075. dot = function( t, x ) 
  3076. {
  3077.   global (a, b, K)
  3078.   return (a-b*K)*x + b*K*([1,0,1,0]');
  3079. };
  3080.  
  3081. x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
  3082.  
  3083. m = maxi( x[;2] );
  3084.  
  3085. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  3086.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) 
  3087. {
  3088.   printf( "\tfailed***\n" );
  3089.   else
  3090.   printf( "\tpassed the lqr/ode test...\n" );
  3091. }
  3092.  
  3093. printf("Elapsed time = %10.3f seconds\n", toc() );
  3094. "FINISHED TESTS"
  3095.