home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / test < prev    next >
Encoding:
Text File  |  1995-02-25  |  70.9 KB  |  1,445 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  that in the
  107. //      original matrix of singular values the order of the diagonal entries
  108. //      is reversed: small to large instead of large to small.
  109. //      KL and KU are the lower and upper bandwidths respectively; if they
  110. //      are omitted a full matrix is produced.
  111. //      If only KL is present, KU defaults to KL.
  112. //      Special case: if KAPPA < 0 then a random full symmetric positive
  113. //                    definite matrix is produced with COND(A) = -KAPPA and
  114. //                    eigenvalues distributed according to MODE.
  115. //                    KL and KU, if present, are ignored.
  116.  
  117. //      This routine is similar to the more comprehensive Fortran routine xLATMS
  118. //      in the following reference:
  119. //      J.W. Demmel and A. McKenney, A test matrix generation suite,
  120. //      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
  121. //      New York, 1989.
  122.  
  123. //    This file is a translation of randsvd.m from version 2.0 of
  124. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  125. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  126.  
  127. // Dependencies
  128. #   rfile bandred
  129. #   rfile qmult
  130.  
  131. //-------------------------------------------------------------------//
  132.  
  133. randsvd = function (n, kappa, mode, kl, ku)
  134. {
  135.   local (n, kappa, mode, kl, ku)
  136.  
  137.   if (!exist (kappa)) { kappa = sqrt(1/eps); }
  138.   if (!exist (mode))  { mode = 3; }
  139.   if (!exist (kl))    { kl = n-1; }    // Full matrix.
  140.   if (!exist (ku))    { ku = kl; }     // Same upper and lower bandwidths.
  141.  
  142.   if (abs(kappa) < 1) {
  143.     error("Condition number must be at least 1!");
  144.   }
  145.  
  146.   posdef = 0;
  147.   if (kappa < 0)            // Special case.
  148.   {
  149.     posdef = 1;
  150.     kappa = -kappa;
  151.   }
  152.  
  153.   p = min(n);
  154.   m = n[1];        // Parameter n specifies dimension: m-by-n.
  155.   n = n[max(size(n))];
  156.  
  157.   if (p == 1)        // Handle case where A is a vector.
  158.   {
  159.     rand("normal", -10, 10);
  160.     A = rand(m, n);
  161.     A = A/norm(A, "2");
  162.     return A;
  163.   }
  164.  
  165.   j = abs(mode);
  166.  
  167.   // Set up vector sigma of singular values.
  168.   if (j == 3)
  169.   {
  170.     factor = kappa^(-1/(p-1));
  171.     sigma = factor.^[0:p-1];
  172.  
  173.   else if (j == 4) {
  174.     sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
  175.  
  176.   else if (j == 5) {    // In this case cond(A) <= kappa.
  177.     rand("uniform", 0, 1)
  178.     sigma = exp( -rand(p,1)*log(kappa) );
  179.  
  180.   else if (j == 2) {
  181.     sigma = ones(p,1);
  182.     sigma[p] = 1/kappa;
  183.  
  184.   else if (j == 1) {
  185.     sigma = ones(p,1)./kappa;
  186.     sigma[1] = 1;
  187.  
  188.   }}}}}
  189.  
  190.  
  191.   // Convert to diagonal matrix of singular values.
  192.   if (mode < 0) {
  193.     sigma = sigma[p:1:-1];
  194.   }
  195.  
  196.   sigma = diag(sigma);
  197.  
  198.   if (posdef)        // Handle special case.
  199.   {
  200.     Q = qmult(p);
  201.     A = Q'*sigma*Q;
  202.     A = (A + A')/2;    // Ensure matrix is symmetric.
  203.     return A;
  204.   }
  205.  
  206.   if (m != n) 
  207.   {
  208.     sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  209.   }
  210.  
  211.   if (kl == 0 && ku == 0) // Diagonal matrix requested - nothing more to do.
  212.   {
  213.     A = sigma;
  214.     return A;
  215.   }
  216.  
  217.   // A = U*sigma*V, where U, V are random orthogonal matrices from the
  218.   // Haar distribution.
  219.   A = qmult(sigma');
  220.   A = qmult(A');
  221.  
  222.   if (kl < n-1 || ku < n-1)    // Bandwidth reduction.
  223.   {
  224.    A = bandred(A, kl, ku);
  225.   }
  226.  
  227.   rand("default");
  228.   return A;
  229. };
  230.  
  231. //-------------------------------------------------------------------//
  232.  
  233. // Synopsis:    Band reduction by two-sided unitary transformations.
  234.  
  235. // Syntax:    bandred ( A , KL, KU )
  236.  
  237. // Description:
  238.  
  239. //    bandred(A, KL, KU) is a matrix unitarily equivalent to A with
  240. //    lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
  241. //    i > j+KL or j > i+KU).  The reduction is performed using
  242. //    Householder transformations. If KU is omitted it defaults to
  243. //    KL. 
  244.  
  245. //    Called by randsvd.
  246. //    This is a `standard' reduction.  Cf. reduction to bidiagonal
  247. //    form prior to computing the SVD.  This code is a little
  248. //    wasteful in that it computes certain elements which are
  249. //    immediately set to zero! 
  250. //
  251. //      Reference:
  252. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  253. //      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  254. //      Section 5.4.3.
  255.  
  256. //    This file is a translation of bandred.m from version 2.0 of
  257. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  258. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  259.  
  260. // Dependencies
  261. #   rfile house
  262.  
  263. //-------------------------------------------------------------------//
  264.  
  265. bandred = function ( A , kl , ku )
  266. {
  267.   local (A, kl, ku)
  268.  
  269.   if (!exist (ku)) { ku = kl; else ku = ku; }
  270.  
  271.   if (kl == 0 && ku == 0) {
  272.     error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  273.   }
  274.  
  275.   // Check for special case where order of left/right transformations matters.
  276.   // Easiest approach is to work on the transpose, flipping back at the end.
  277.  
  278.   flip = 0;
  279.   if (ku == 0)
  280.   {
  281.     A = A';
  282.     temp = kl; kl = ku; ku = temp; flip = 1;
  283.   }
  284.  
  285.   m = A.nr; n = A.nc; 
  286.  
  287.   for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  288.   {
  289.     if (j+kl+1 <= m)
  290.     {
  291.        ltmp = house(A[j+kl:m;j]);
  292.        beta = ltmp.beta; v = ltmp.v;
  293.        temp = A[j+kl:m;j:n];
  294.        A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
  295.        A[j+kl+1:m;j] = zeros(m-j-kl,1);
  296.     }
  297.  
  298.     if (j+ku+1 <= n)
  299.     {
  300.        ltmp = house(A[j;j+ku:n]');
  301.        beta = ltmp.beta; v = ltmp.v;
  302.        temp = A[j:m;j+ku:n];
  303.        A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
  304.        A[j;j+ku+1:n] = zeros(1,n-j-ku);
  305.     }
  306.   }
  307.  
  308.   if (flip) {
  309.     A = A';
  310.   }
  311.  
  312.   return A;
  313. };
  314.  
  315. //-------------------------------------------------------------------//
  316.  
  317. // Synopsis:    Lehmer matrix - symmetric positive definite.
  318.  
  319. // Syntax:      A = lehmer ( N )
  320.  
  321. // Description:
  322.  
  323. //      A is the symmetric positive definite N-by-N matrix with
  324. //                     A(i,j) = i/j for j >= i.
  325. //      A is totally nonnegative.  INV(A) is tridiagonal, and explicit
  326. //      formulas are known for its entries. 
  327.  
  328. //      N <= COND(A) <= 4*N*N.
  329.  
  330. //      References:
  331. //        M. Newman and J. Todd, The evaluation of matrix inversion
  332. //           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
  333. //        Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
  334. //           of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
  335. //        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  336. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.
  337.  
  338. //    This file is a translation of lehmer.m from version 2.0 of
  339. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  340. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  341.  
  342. //-------------------------------------------------------------------//
  343.  
  344. lehmer = function ( n )
  345. {
  346.   local (n)
  347.   global (tril)
  348.  
  349.   A = ones(n,1)*(1:n);
  350.   A = A./A';
  351.   A = tril(A) + tril(A,-1)';
  352.  
  353.   return A;
  354. };
  355.  
  356. //-------------------------------------------------------------------//
  357.  
  358. // Synopsis:    Pre-multiply by random orthogonal matrix.
  359.  
  360. // Syntax:    B = qmult ( A )
  361.  
  362. // Description:
  363.  
  364. //    B is Q*A where Q is a random real orthogonal matrix from the
  365. //    Haar distribution, of dimension the number of rows in
  366. //    A. Special case: if A is a scalar then QMULT(A) is the same as
  367.  
  368. //        qmult(eye(a))
  369.  
  370. //       Called by RANDSVD.
  371.  
  372. //       Reference:
  373. //       G.W. Stewart, The efficient generation of random
  374. //       orthogonal matrices with an application to condition estimators,
  375. //       SIAM J. Numer. Anal., 17 (1980), 403-409.
  376.  
  377. //    This file is a translation of qmult.m from version 2.0 of
  378. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  379. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  380.  
  381. //-------------------------------------------------------------------//
  382.  
  383. qmult = function ( A )
  384. {
  385.   local (A)
  386.  
  387.   n = A.nr; m = A.nc;
  388.  
  389.   //  Handle scalar A.
  390.   if (max(n,m) == 1)
  391.   {
  392.     n = A;
  393.     A = eye(n,n);
  394.   }
  395.  
  396.   d = zeros(n,n);
  397.  
  398.   for (k in n-1:1:-1)
  399.   {
  400.     // Generate random Householder transformation.
  401.     rand("normal", 0, 1);
  402.     x = rand(n-k+1,1);
  403.     s = norm(x, "2");
  404.     sgn = sign(x[1]) + (x[1]==0);    // Modification for sign(1)=1.
  405.     s = sgn*s;
  406.     d[k] = -sgn;
  407.     x[1] = x[1] + s;
  408.     beta = s*x[1];
  409.  
  410.     // Apply the transformation to A.
  411.     y = x'*A[k:n;];
  412.     A[k:n;] = A[k:n;] - x*(y/beta);
  413.   }
  414.  
  415.   // Tidy up signs.
  416.   for (i in 1:n-1)
  417.   {
  418.     A[i;] = d[i]*A[i;];
  419.   }
  420.  
  421.   A[n;] = A[n;]*sign(rand());
  422.   B = A;
  423.  
  424.   rand("default");
  425.   return B;
  426. };
  427.  
  428. //-------------------------------------------------------------------//
  429.  
  430. // Synopsis:    Create a Householder matrix.
  431.  
  432. // Syntax:    house ( X )
  433.  
  434. //    If HOUSE(x), which returns a list containing elements `v' and
  435. //    `beta',  then H = EYE - beta*v*v' is a Householder matrix such
  436. //    that Hx = -sign(x(1))*norm(x)*e_1. 
  437. //    NB: If x = 0 then v = 0, beta = 1 is returned.
  438. //            x can be real or complex.
  439. //            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
  440.  
  441. //    Theory: (textbook references Golub & Van Loan 1989, 38-43;
  442. //             Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
  443. //      Hx = y: (I - beta*v*v')x = -s*e_1.
  444. //      Must have |s| = norm(x), v = x+s*e_1, and
  445. //      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
  446. //      So take s = sign(x(1))*norm(x) (which avoids cancellation).
  447. //      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
  448. //          = 2*norm(x)*(norm(x) + |x(1)|).
  449.  
  450. //    References:
  451. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  452. //           Johns Hopkins University Press, Baltimore, Maryland, 1989.
  453. //        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  454. //           New York, 1973,
  455. //        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  456. //           Press, 1965.
  457.  
  458. //    This file is a translation of house.m from version 2.0 of
  459. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  460. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  461.  
  462. //-------------------------------------------------------------------//
  463.  
  464. house = function ( x )
  465. {
  466.   local (x)
  467.  
  468.   m = x.nr; n = x.nc;
  469.   if (n > 1) { error ("Argument must be a column vector."); }
  470.  
  471.   s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  472.   v = x;
  473.   if (s == 0)    // Quit if x is the zero vector.
  474.   {
  475.     beta = 1; 
  476.     return << beta = beta ; v = v >>;
  477.   }
  478.  
  479.   v[1] = v[1] + s;
  480.   beta = 1/(s'*v[1]);        // NB the conjugated s.
  481.  
  482.   // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  483.   // But beta as above can be non-real (due to rounding) only 
  484.   // when x is complex.
  485.  
  486.   return << beta = beta ; v = v >>;
  487. };
  488.  
  489. //-------------------------------------------------------------------//
  490.  
  491. //  Syntax:    tril ( A )
  492. //        tril ( A , K )
  493.  
  494. //  Description:
  495.  
  496. //  tril(x) returns the lower triangular part of A.
  497.  
  498. //  tril(A,K) returns the elements on and below the K-th diagonal of
  499. //  A.
  500.  
  501. //  K = 0: main diagonal
  502. //  K > 0: above the main diag.
  503. //  K < 0: below the main diag.
  504.  
  505. //  See Also: triu
  506. //-------------------------------------------------------------------//
  507.  
  508. tril = function(x, k) 
  509. {
  510.   local(x, k)
  511.  
  512.   if (!exist (k)) { k = 0; }
  513.   nr = x.nr; nc = x.nc;
  514.   if(k > 0) 
  515.   { 
  516.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  517.   else
  518.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  519.   }
  520.  
  521.   y = zeros(nr, nc);
  522.  
  523.   for(i in max( [1,1-k] ):nr) {
  524.     j = 1:min( [nc, i+k] );
  525.     y[i;j] = x[i;j];
  526.   }
  527.  
  528.   return y;
  529. };
  530.  
  531. //-------------------------------------------------------------------//
  532.  
  533. //  Syntax:    triu ( A )
  534. //        triu ( A , K )
  535.  
  536. //  Description:
  537.  
  538. //  triu(x) returns the upper triangular part of A.
  539.  
  540. //  tril(x; k) returns the elements on and above the k-th diagonal of
  541. //  A. 
  542.  
  543. //  K = 0: main diagonal
  544. //  K > 0: above the main diag.
  545. //  K < 0: below the main diag.
  546.  
  547. //  See Also: tril
  548. //-------------------------------------------------------------------//
  549.  
  550. triu = function(x, k) 
  551. {
  552.   local(x, k)
  553.  
  554.   if (!exist (k)) { k = 0; }
  555.   nr = x.nr; nc = x.nc;
  556.  
  557.   if(k > 0) 
  558.   { 
  559.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  560.   else
  561.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  562.   }
  563.  
  564.   y = zeros(nr, nc);
  565.  
  566.   for(j in max( [1,1+k] ):nc) {
  567.     i = 1:min( [nr, j-k] );
  568.     y[i;j] = x[i;j];
  569.   }
  570.  
  571.   return y;
  572. };
  573.  
  574. //
  575. //-------------------- Test Relational Expressions -------------------
  576. //
  577. printf("\tstart scalar tests...\n");
  578. printf("\tstart relational tests...\n");
  579.  
  580. //    SCALAR CONSTANTS (REAL)
  581. if( !(1<2) ) { error(); }
  582. if( !(1<=2) ) { error(); }
  583. if( 1>2 ) { error(); }
  584. if( 1>=2 ) { error(); }
  585. if( 1==2 ) { error() != "real") { error (); }
  586.  
  587. //
  588. // Test row-wise matrix addition
  589. //
  590.  
  591. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  592. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  593. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  594.  
  595. a = [1,2,3] + [1,2,3]*1i; 
  596. 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;
  597. 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;
  598. if (!all (all (a .+ b == c))) { error (); }
  599. if (!all (all (b .+ a == c))) { error (); }
  600.  
  601. printf("\t\tpassed matrix row-wise add test...\n");
  602.  
  603. //
  604. // Test row-wise matrix subtraction
  605. //
  606.  
  607. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  608. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  609. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  610.  
  611. a = [1,1,1] + [1,1,1]*1i;
  612. 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;
  613. 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;
  614. if (!all (all (a .- b == -c))) { error (); }
  615. if (!all (all (b .- a ==  c))) { error (); }
  616.  
  617. printf("\t\tpassed matrix row-wise subtraction test...\n");
  618.  
  619. //
  620. // Test col-wise matrix addition
  621. //
  622.  
  623. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  624. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  625. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  626.  
  627. a = [1;1;1;1] + [1;1;1;1]*1i;
  628. 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;
  629. 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;
  630. if (!all (all (a .+ b == c))) { error (); }
  631. if (!all (all (b .+ a == c))) { error (); }
  632.  
  633. printf("\t\tpassed matrix col-wise add test...\n");
  634.  
  635. //
  636. // Test col-wise matrix subtraction
  637. //
  638.  
  639. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  640. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  641. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  642.  
  643. a = [1;1;1;1] + [1;1;1;1]*1i;
  644. 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;
  645. 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;
  646. if (!all (all (a .- b == -c))) { error (); }
  647. if (!all (all (b .- a ==  c))) { error (); }
  648.  
  649. printf("\t\tpassed matrix col-wise subtraction test...\n");
  650.  
  651. a = [1,2,3];
  652. b = [1,2,3;4,5,6;7,8,9];
  653. c = [1,4,9;4,10,18;7,16,27];
  654.  
  655. if (!all (all (a .* b == c))) { error (); }
  656. if (!all (all (b .* a == c))) { error (); }
  657.  
  658. za = a + rand (size (a))*1j;
  659. zb = b + rand (size (b))*1j;
  660.  
  661. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  662. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  663. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  664. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  665. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  666. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  667.  
  668. printf("\t\tpassed matrix row-wise multiplication test...\n");
  669.  
  670. a = [1,2,3];
  671. b = [1,2,3;4,6,6;7,8,9];
  672. c = [1,1,1;4,3,2;7,4,3];
  673.  
  674. if (!all (all (b ./ a == c))) { error (); }
  675. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  676. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  677.  
  678. za = a + rand (size (a))*1j;
  679. zb = b + rand (size (b))*1j;
  680.  
  681. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  682. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  683. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  684. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  685. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  686. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  687.  
  688. printf("\t\tpassed matrix row-wise division test...\n");
  689.  
  690. a = [1;2;3];
  691. b = [1,2,3;4,5,6;7,8,9];
  692.  
  693. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  694. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  695.  
  696. za = a + rand (size (a))*1j;
  697. zb = b + rand (size (b))*1j;
  698.  
  699. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  700. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  701. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  702. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  703. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  704. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  705.  
  706. printf("\t\tpassed matrix column-wise multiplication test...\n");
  707.  
  708. a = [1;2;3];
  709. b = [1,2,3;4,6,6;7,8,9];
  710.  
  711. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  712. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  713.  
  714. za = a + rand (size(a))*1j;
  715. zb = b + rand (size(b))*1j;
  716.  
  717. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  718. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  719. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  720. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  721. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  722. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  723.  
  724. printf("\t\tpassed matrix column-wise division test...\n");
  725.  
  726.  
  727. //
  728. //--------------------- Test COMPLEX MATRICES --------------------------
  729. //
  730. //  Automatic creation, extension
  731. //
  732. printf("\t\tcomplex-matrices\n");
  733. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  734. a=[1,2,3;4,5,6];
  735. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  736. //
  737. a = m3;
  738. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  739. if(any (any (a-a != zeros(2,2)))) { error(); }
  740. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
  741. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  742. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  743. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  744. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  745. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  746. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  747. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  748. //
  749. // The following test may not work on some machines
  750. //
  751. if(any (any (a./a != [1,1;1,1]))) { 
  752.   printf("\t\t***complex division inaccuracy, check manually***\n");
  753. }
  754.  
  755. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  756. //  
  757. //--------------------- Test NULL MATRICES -------------------------
  758. //
  759. printf("\t\tnull-matrices\n");
  760. // Create a NULL matrix
  761. mnull = [];
  762. // Test it with SCALARS
  763. if( any([1,mnull] != 1)) {
  764.   error();
  765. }
  766. if( any([mnull,1] != 1)) {
  767.   error();
  768. }
  769. // Test with MATRIX construction
  770. m = [1,2;3,4;5,6];
  771. if( any([mnull;1] != [1])) {
  772.   error();
  773. }
  774. if( any([1;mnull] != [1])) {
  775.   error();
  776. }
  777. if( any([mnull;1,2,3] != [1,2,3])) {
  778.   error();
  779. }
  780. if( any([1,2,3;mnull] != [1,2,3])) {
  781.   error();
  782. }
  783. if(any (any ([mnull,m] != m))) {
  784.   error();
  785. }
  786. if(any (any ([m,mnull] != m))) {
  787.   error();
  788. }
  789. if(any (any ([mnull;m] != m))) {
  790.   error();
  791. }
  792. if(any (any ([m;mnull] != m))) {
  793.   error();
  794.  
  795. mnull = matrix();
  796. mnull[1] = [1];
  797. }
  798.  
  799. //--------------------- Test Matrix Multiply --------------------------
  800.  
  801. i = sqrt(-1);
  802. a = [1,2,3;4,5,6;7,8,9];
  803. b = [4,5,6;7,8,9;10,11,12];
  804. c = [ 48,  54,  60 ;
  805.      111, 126, 141 ;
  806.      174, 198, 222 ];
  807.  
  808. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  809.  
  810. az = a + b*i;
  811. bz = b + a*i;
  812.  
  813. cz = [-18+141*i , -27+162*i , -36+183*i ;
  814.         9+240*i ,   0+279*i ,  -9+318*i ;
  815.        36+339*i ,  27+396*i ,  18+453*i ];
  816.  
  817. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  818.        111+66*i , 126+81*i  , 141+96*i ;
  819.        174+102*i, 198+126*i , 222+150*i ];
  820.  
  821. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  822.         111+174*i , 126+198*i , 141+222*i ;
  823.         174+237*i , 198+270*i , 222+303*i ];
  824.  
  825. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  826. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  827. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  828.  
  829. a = [a,a];
  830. b = [b;b];
  831. c = [  96 , 108 , 120 ;
  832.       222 , 252 , 282 ;
  833.       348 , 396 , 444 ];
  834.  
  835. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  836.  
  837. az = [az,az];
  838. bz = [bz;bz];
  839.  
  840. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  841.          18+480*i ,    0+558*i ,  -18+636*i ;
  842.          72+678*i ,   54+792*i ,   36+906*i ];
  843.  
  844. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  845.         222+132*i , 252+162*i , 282+192*i ;
  846.         348+204*i , 396+252*i , 444+300*i ];
  847.  
  848. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  849.          222+348*i , 252+396*i , 282+444*i ;
  850.          348+474*i , 396+540*i , 444+606*ion ( A )
  851. {
  852.   local (i, l, u, pvt, x)
  853.  
  854.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  855.  
  856.   x = factor (A, "g");    // Do the factorization
  857.  
  858.   //
  859.   // Now create l, u, and pvt from lu and pvt.
  860.   //
  861.  
  862.   l = tril (x.lu, -1) + eye (size (x.lu));
  863.   u = triu (x.lu);
  864.   pvt = eye (size (x.lu));
  865.  
  866.   //
  867.   // Now re-arange the columns of pvt
  868.   //
  869.  
  870.   for (i in 1:max (size (x.lu)))
  871.   {
  872.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  873.   }
  874.   return << l = l; u = u; pvt = pvt >>;
  875. };
  876.  
  877. //
  878. //  In vector V, swap elements I, J
  879. //
  880.  
  881. swap = function ( V, I, J )
  882. {
  883.   local (v, tmp);
  884.   v = V;
  885.   tmp = v[I];
  886.   v[I] = v[J];
  887.   v[J] = tmp;
  888.   return v;
  889. };
  890.  
  891. a = randsvd(10,10);
  892. lua = lu (a);
  893. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  894. { 
  895.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  896.   printf ("\tA - p*l*u:\n");
  897.   abs (a - lua.pvt*lua.l*lua.u)
  898.   error ("possible lu()/factor() problems\n");
  899. }
  900.  
  901. //
  902. // Real
  903. az = randsvd(10,10) + randsvd(10,10)*1j;
  904. luz = lu (az);
  905. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  906. { 
  907.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  908.   printf ("\tA - p*l*u:\n");
  909.   abs (az - luz.ovt*luz.l*luz.u)
  910.   error ("possible lu()/factor()() problems\n");
  911. }
  912.  
  913. printf("\tpassed lu/factor test...\n");
  914.  
  915. //
  916. //-------------------------- Test svd ()   -----------------------------
  917. //
  918.  
  919. a = randsvd(10,10);
  920. s = svd (a);
  921. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
  922. {
  923.   error ("possible svd() problems"); 
  924. }
  925.  
  926. z = randsvd(10,10) + rand(10,10)*1j;
  927. sz = svd (z);
  928. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
  929. { 
  930.   error ("possible svd() problems"); 
  931. }
  932.  
  933. printf("\tpassed svd test...\n");
  934.  
  935. //
  936. //-------------------------- Test hess ()  -----------------------------
  937. //
  938.  
  939. a = randsvd(10,10);
  940. h = hess (a);
  941. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  942. { 
  943.   error ("possible hess() problems");
  944. }
  945.  
  946. z = randsvd(10,10) + randsvd(10,10)*1j;
  947. hz = hess (z);
  948. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  949. { 
  950.   error ("possible hess() problems"); 
  951. }
  952.  
  953. printf("\tpassed hess test...\n");
  954.  
  955. //
  956. //-------------------------- Test lyap () ------------------------------
  957. //
  958.  
  959. lyap = function ( A, B, C )
  960. {
  961.   local (A, B, C)
  962.  
  963.   if (!exist (B)) 
  964.   { 
  965.     B = A';    // Solve the special form: A*X + X*A' = -C
  966.   }
  967.  
  968.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  969.     error ("Dimensions do not agree.");
  970.   }
  971.  
  972.   //
  973.   // Schur decomposition on A and B
  974.   //
  975.  
  976.   sa = schur (A);
  977.   sb = schur (B);
  978.  
  979.   //
  980.   // transform C
  981.   //
  982.  
  983.   tc = sa.z' * C * sb.z;
  984.  
  985.   X = sylv (sa.t, sb.t, tc);
  986.  
  987.   //
  988.   // Undo the transformation
  989.   //
  990.  
  991.   X = sa.z * X * sb.z';
  992.  
  993.   return X;
  994. };
  995.  
  996. a = randsvd (10,10);
  997. b = rand (10,10);
  998. c = rand (10,10);
  999.  
  1000. x = lyap (a, b, c);
  1001. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  1002. { 
  1003.   error ("possible problems with lyap() or sylv()"); 
  1004. }
  1005.  
  1006. printf("\tpassed lyap test...\n");
  1007.  
  1008. //
  1009. //-------------------------- Test eig () ------------------------------
  1010. //
  1011.  
  1012. trace = function(m) 
  1013. {
  1014.   local(i, tr);
  1015.  
  1016.   if(m.class != "num") { 
  1017.     error("must provide NUMERICAL input to trace()");
  1018.   }
  1019.  
  1020.   tr = 0;
  1021.   for(i in 1:min( [m.nr, m.nc] )) {
  1022.     tr = tr + m[i;i];
  1023.   }
  1024.  
  1025.   return tr;
  1026. };
  1027.  
  1028. eye = function( m , n ) 
  1029. {
  1030.   local(i, N, new);
  1031.  
  1032.   if (!exist (n))
  1033.   {
  1034.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  1035.     new = zeros (m[1], m[2]);
  1036.     N = min ([m[1], m[2]]);
  1037.   else
  1038.     if (class (m) == "string" || class (n) == "string") {
  1039.       error ("eye(), string arguments not allowed");
  1040.     }
  1041.     if (max (size (m)) == 1 && max (size (n)) == 1)
  1042.     {
  1043.       new = zeros (m[1], n[1]);
  1044.       N = min ([m[1], n[1]]);
  1045.     else
  1046.       error ("matrix arguments to eye() must be 1x1");
  1047.     }
  1048.   }
  1049.   for(i in 1:N)
  1050.   {
  1051.     new[i;i] = 1.0;
  1052.   }
  1053.   return new;
  1054. };
  1055.  
  1056. //
  1057. // Standard eigenvalue problem
  1058. //
  1059.  
  1060. a = randsvd(10,10);
  1061. ta = trace (a);
  1062. sa = symm (a);
  1063. tsa = trace (sa);
  1064.  
  1065. z = randsvd(10,10) + randsvd(10,10)*1i;
  1066. tz = trace (z);
  1067. sz = symm (z);
  1068. tsz = trace (sz);
  1069.  
  1070. tol = 1.e-6;
  1071.  
  1072. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  1073. {
  1074.   error ("error in eig");
  1075. }
  1076. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  1077. {
  1078.   error ("error in eig");
  1079. }
  1080. if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
  1081. {
  1082.   error ("error in eig");
  1083. }
  1084. if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
  1085. {
  1086.   error ("error in eig");
  1087. }
  1088.  
  1089. //
  1090. // Generalized eigenvalue problem
  1091. //
  1092.  
  1093. b = randsvd(10,10);
  1094. sb = symm (b) + eye(size(b))*3;
  1095. tb = trace (b);
  1096. tsb = trace (sb);
  1097.  
  1098. zb = randsvd(10,10) + randsvd(10,10)*1i;
  1099. szb = symm (zb) + eye(size(zb))*3;
  1100. tzb = trace (zb);
  1101. tszb = trace (szb);
  1102.  
  1103. eig(a,b);    // not sure of a good way to check these yet
  1104. eigs(sa,sb);
  1105. eigs(sa,sb);
  1106.  
  1107. eig(z, zb);
  1108. eigs(sz, szb);
  1109. eigs(sz, szb);
  1110.  
  1111. printf("\tpassed eig test...\n");
  1112.  
  1113. //
  1114. //-------------------------- Test fft () -----------------------------
  1115. //
  1116.  
  1117. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  1118. printf("\tpassed fft test...\n");
  1119.  
  1120. //
  1121. //------------------------- Fibonacci Test -------------------------
  1122. //
  1123. //  Calculate Fibonacci numbers
  1124. //
  1125.  
  1126. i=1; 
  1127. while ( i < 2 ) { 
  1128.   i=i+1;
  1129.   a=0; b=1;
  1130.   while ( b < 10000 ) {
  1131.     c = b;
  1132.     b = a+b;
  1133.     a = c;
  1134.   }
  1135. }
  1136. if ( b != 10946 ) {
  1137.   error("failed fibonacci test");
  1138. else
  1139.   printf("\tpassed fibonacci test...\n");
  1140. }
  1141.  
  1142. //
  1143. //------------------------- Factorial Test -------------------------
  1144. //
  1145.  
  1146. fac = function(a) 
  1147. {
  1148.   if(a <= 1) 
  1149.   {
  1150.     return 1
  1151.   else
  1152.     return a*$self(a-1)
  1153.   }
  1154. };
  1155.  
  1156. if(fac(10) != 3628800)
  1157. { 
  1158.   error(); else printf("\tpassed factorial test...\n"); 
  1159. }
  1160.  
  1161. //
  1162. //--------------------------- ACK Test ----------------------------
  1163. //
  1164.  
  1165. ack = function(a, b) 
  1166. {
  1167.   if(a == 0) { return b + 1; }
  1168.   if(b == 0) { return $self(a-1, 1); }
  1169.   return $self(a-1, $self(a, b-1));
  1170. };
  1171.  
  1172. if(ack(2,2) != 7) 
  1173. {
  1174.   error("error in ack() test");
  1175.   else
  1176.   printf("\tpassed ACK test...\n");
  1177. }
  1178.  
  1179. //
  1180. //------------------------- Prime Test -----------------------------
  1181. //
  1182. // An example that finds all primes less than limit
  1183. //
  1184.  
  1185. primes = function (limit) 
  1186. {
  1187.   local(prime, cnt, i, j, k);
  1188.   
  1189.   i = 1; cnt = 0;
  1190.   for(k in 2:limit) 
  1191.   {
  1192.     j = 2;
  1193.     while(mod(k,j) != 0) 
  1194.     {
  1195.       j++;
  1196.     }
  1197.     if(j == k)             // Found prime
  1198.     {
  1199.       cnt++;
  1200.       prime[i;1] = k;
  1201.       i++;
  1202.     }
  1203.   }
  1204.   return prime;
  1205. };
  1206.  
  1207. if(max(size(primes(100))) != 25) 
  1208. {  
  1209.   error("error in prime test");
  1210.   else
  1211.   printf("\tpassed prime test...\n");
  1212. }
  1213.  
  1214. //
  1215. //--------------------------- Fib Min Test -----------------------------
  1216. //
  1217. //    fibmin() will minimize an arbitrary function 
  1218. //    in 1D using Fibonacci search
  1219.  
  1220. f065 = function(x)
  1221. {
  1222.   return (x - 0.65) * (x - 0.65);
  1223. };
  1224.  
  1225. fib = function(x)
  1226. {
  1227.   local (n, a, b);
  1228.   
  1229.   a = 1;
  1230.   b = 1;
  1231.   if (x >= 2)
  1232.   {
  1233.     n = x - 1;
  1234.     for (n in n:1:-1)
  1235.     {
  1236.       c = b;
  1237.       b = a + b;
  1238.       a = c;
  1239.       n = n - 1;
  1240.     }
  1241.   }
  1242.   return b;
  1243. };
  1244.  
  1245. //  Minimize a 1D function using Fibonacci search
  1246. //  f = function to minimize
  1247. //  xlo = lower bound
  1248. //  xhi = upper bound
  1249. //  n = number of iterations (the bigger the more accurate)
  1250.  
  1251. fibmin = function(f, xlo, xhi, n) 
  1252. {
  1253.   local(a, b, x, y, ex, ey, k, lo, hi);
  1254.   
  1255.   lo = xlo;
  1256.   hi = xhi;
  1257.   k = n;
  1258.   for (k in k:2:-1)
  1259.   {
  1260.     a = fib(k - 2) / fib(k);
  1261.     b = fib(k - 1) / fib(k);
  1262.     x = lo + (hi - lo) * a;
  1263.     y = lo + (hi - lo) * b;
  1264.     ex = f(x);
  1265.     ey = f(y);
  1266.     if (ex >= ey)
  1267.     {
  1268.       lo = x;
  1269.       else
  1270.       hi = y;
  1271.     }
  1272.     //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  1273.   }
  1274.   return (lo + hi) / 2;
  1275. };
  1276.  
  1277. //
  1278. // Simple example using above function to mimize f065. Answer is 0.65
  1279. //
  1280.  
  1281. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  1282. if (abs(x - 0.65) > 1e-6)
  1283. {
  1284.   printf("x = %f\n", x);
  1285.   error("failed fibmin test");
  1286. }
  1287.  
  1288. printf("\tpassed fibmin test...\n");
  1289.  
  1290. //
  1291. //--------------------- Nasty Function Test ------------------------
  1292. //
  1293.  
  1294. printf("\tStarting Nasty Function Test...");
  1295. printf("\tthis will take awhile\n");
  1296. check = function( a, b, c, d, e, f, g, h ) {
  1297.   if ( a+b+c+d == e+f+g+h && ...
  1298.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  1299.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
  1300.     return 1;
  1301.   else
  1302.     return 0;
  1303.   }
  1304. };
  1305.  
  1306. for(a in 8:10) {
  1307.   for(b in 7:(a-1)) {
  1308.     for(c in 6:(b-1)) {
  1309.       for(d in 5:(c-1)) {
  1310.         for(e in 4:(d-1)) {
  1311.           for(f in 3:(e-1)) {
  1312.             for(g in 2:(f-1)) {
  1313.               for(h in 1:(g-1)) {                  
  1314.               if(check( a, b, c, d,  e, f, g, h ) || ...
  1315.                      check( a, e, c, d,  b, f, g, h ) || ...
  1316.                      check( a, f, c, d,  e, b, g, h ) || ...
  1317.                      check( a, g, c, d,  e, f, b, h ) || ...
  1318.                      check( a, h, c, d,  e, f, g, b ) || ...
  1319.                      check( a, b, e, d,  c, f, g, h ) || ...
  1320.                      check( a, b, f, d,  e, c, g, h ) || ...
  1321.                      check( a, b, g, d,  e, f, c, h ) || ...
  1322.                      check( a, b, h, d,  e, f, g, c ) || ...
  1323.                      check( a, b, c, e,  d, f, g, h ) || ...
  1324.                      check( a, b, c, f,  e, d, g, h ) || ...
  1325.                      check( a, b, c, g,  e, f, d, h ) || ...
  1326.                      check( a, b, c, h,  e, f, g, d ) || ...
  1327.                      check( a, e, f, d,  b, c, g, h ) || ...
  1328.                      check( a, e, g, d,  b, f, c, h ) || ...
  1329.                      check( a, e, h, d,  b, f, g, c ) || ...
  1330.                      check( a, f, g, d,  e, b, c, h ) || ...
  1331.                      check( a, f, h, d,  e, b, g, c ) || ...
  1332.                      check( a, g, h, d,  e, f, b, c ) || ...
  1333.                      check( a, b, e, f,  c, d, g, h ) || ...
  1334.                      check( a, b, e, g,  c, f, d, h ) || ...
  1335.                      check( a, b, e, h,  c, f, g, d ) || ...
  1336.                      check( a, b, f, g,  e, c, d, h ) || ...
  1337.                      check( a, b, f, h,  e, c, g, d ) || ...
  1338.                      check( a, b, g, h,  e, f, c, d ) || ...
  1339.                      check( a, e, f, g,  e, f, g, h ) || ...
  1340.                      check( a, e, f, h,  e, f, g, h ) || ...
  1341.                      check( a, e, g, h,  e, f, g, h ) || ...
  1342.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  1343.               }
  1344.             }
  1345.           }
  1346.         }
  1347.       }
  1348.     }
  1349.   }
  1350. }
  1351.  
  1352. if(1) {  // figure out the value of cnt, and check!
  1353.   printf("\tpassed nasty function test...\n");
  1354. else
  1355.   error();
  1356. }
  1357.  
  1358. //
  1359. //------------------ Test More Advanced Functions --------------------
  1360. //
  1361.  
  1362. printf( "\tStarting the lqr/ode test..." );
  1363. printf( "\tthis will take awhile\n" );
  1364.  
  1365. lqr = function( a, b, q, r ) 
  1366. {
  1367.   local( k, s,...
  1368.          m, n, mb, nb, mq, nq,...
  1369.          e, v, d );
  1370.  
  1371.   m = size(a)[1]; n = size(a)[2];
  1372.   mb = size(b)[1]; nb = size(b)[2];
  1373.   mq = size(q)[1]; nq = size(q)[2];
  1374.   
  1375.   if ( m != mq || n != nq ) 
  1376.   {
  1377.     fprintf( "stderr", "A and Q must be the same size.\n" );
  1378.     quit
  1379.   }
  1380.     
  1381.   mr = size(r)[1]; nr = size(r)[2];
  1382.   if ( mr != nr || nb != mr ) 
  1383.   {
  1384.     fprintf( "stderr", "B and R must be consistent.\n" );
  1385.     quit
  1386.   }
  1387.     
  1388.   nn = zeros( m, nb );
  1389.     
  1390.   // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  1391.     
  1392.   e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  1393.   v = e.vec; d = e.val;
  1394.     
  1395.   index = sort( real( d ) ).ind;
  1396.   d = real( d[ index ] );
  1397.   
  1398.   if ( !( d[n] < 0 && d[n+1] > 0 ) ) 
  1399.   {
  1400.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  1401.     quit
  1402.   }
  1403.   
  1404.   chi = v[ 1:n; index[1:n] ];
  1405.   lambda = v[ (n+1):(2*n); index[1:n] ];
  1406.   s = -real(solve(chi',lambda')');
  1407.   k = solve( r, nn'+b'*s );
  1408.   
  1409.   return << k=k; s=s >>;
  1410.   
  1411. };
  1412.  
  1413. // Now run a little test problem.
  1414.  
  1415. k = 1; m = 1; c = .1;
  1416. a = [0     ,1    ,0    , 0;
  1417.     -k/m, -c/m,  k/m,  c/m;
  1418.      0,     0,    0,    1;
  1419.      k/m,  c/m, -k/m, -c/m ];
  1420. b = [ 0; 1/m; 0; 0 ];
  1421. qxx = diag( [0, 0, 100, 0] );
  1422. ruu = [1];
  1423. K = lqr( a, b, qxx, ruu ).k;
  1424.  
  1425. dot = function( t, x ) 
  1426. {
  1427.   global (a, b, K)
  1428.   return (a-b*K)*x + b*K*([1,0,1,0]');
  1429. };
  1430.  
  1431. x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
  1432.  
  1433. m = maxi( x[;2] );
  1434.  
  1435. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  1436.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) 
  1437. {
  1438.   printf( "\tfailed***\n" );
  1439.   else
  1440.   printf( "\tpassed the lqr/ode test...\n" );
  1441. }
  1442.  
  1443. printf("Elapsed time = %10.3f seconds\n", toc() );
  1444. "FINISHED TESTS"
  1445.