home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / liboctave / dMatrix.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  45KB  |  2,380 lines
  1. // Matrix manipulations.                              -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <sys/types.h>
  29. #include <iostream.h>
  30. #include <stdio.h>
  31. #include <float.h>
  32.  
  33. #include <Complex.h>
  34.  
  35. #include "mx-base.h"
  36. #include "dbleDET.h"
  37. #include "dbleSVD.h"
  38. #include "mx-inlines.cc"
  39. #include "lo-error.h"
  40. #include "f77-uscore.h"
  41.  
  42. // Fortran functions we call.
  43.  
  44. extern "C"
  45. {
  46.   int F77_FCN (dgemm) (const char*, const char*, const int*,
  47.                const int*, const int*, const double*,
  48.                const double*, const int*, const double*,
  49.                const int*, const double*, double*, const int*,
  50.                long, long);
  51.  
  52.   int F77_FCN (dgeco) (double*, const int*, const int*, int*, double*,
  53.                double*);
  54.  
  55.   int F77_FCN (dgesl) (const double*, const int*, const int*,
  56.                const int*, double*, const int*); 
  57.  
  58.   int F77_FCN (dgedi) (double*, const int*, const int*, const int*,
  59.                double*, double*, const int*);
  60.  
  61.   int F77_FCN (dgelss) (const int*, const int*, const int*, double*,
  62.             const int*, double*, const int*, double*,
  63.             const double*, int*, double*, const int*,
  64.             int*);
  65.  
  66. // Note that the original complex fft routines were not written for
  67. // double complex arguments.  They have been modified by adding an
  68. // implicit double precision (a-h,o-z) statement at the beginning of
  69. // each subroutine.
  70.  
  71.   int F77_FCN (cffti) (const int*, Complex*);
  72.  
  73.   int F77_FCN (cfftf) (const int*, Complex*, Complex*);
  74.  
  75.   int F77_FCN (cfftb) (const int*, Complex*, Complex*);
  76. }
  77.  
  78. #define KLUDGE_MATRICES
  79. #define TYPE double
  80. #define KL_MAT_TYPE Matrix
  81. #include "mx-kludge.cc"
  82. #undef KLUDGE_MATRICES
  83. #undef TYPE
  84. #undef KL_MAT_TYPE
  85.  
  86. /*
  87.  * Matrix class.
  88.  */
  89.  
  90. Matrix::Matrix (const DiagMatrix& a)
  91.   : Array2<double> (a.rows (), a.cols (), 0.0)
  92. {
  93.   for (int i = 0; i < a.length (); i++)
  94.     elem (i, i) = a.elem (i, i);
  95. }
  96.  
  97. int
  98. Matrix::operator == (const Matrix& a) const
  99. {
  100.   if (rows () != a.rows () || cols () != a.cols ())
  101.     return 0;
  102.  
  103.   return equal (data (), a.data (), length ());
  104. }
  105.  
  106. int
  107. Matrix::operator != (const Matrix& a) const
  108. {
  109.   return !(*this == a);
  110. }
  111.  
  112. Matrix&
  113. Matrix::insert (const Matrix& a, int r, int c)
  114. {
  115.   int a_rows = a.rows ();
  116.   int a_cols = a.cols ();
  117.   if (r < 0 || r + a_rows - 1 > rows ()
  118.       || c < 0 || c + a_cols - 1 > cols ())
  119.     {
  120.       (*current_liboctave_error_handler) ("range error for insert");
  121.       return *this;
  122.     }
  123.  
  124.   for (int j = 0; j < a_cols; j++)
  125.     for (int i = 0; i < a_rows; i++)
  126.       elem (r+i, c+j) = a.elem (i, j);
  127.  
  128.   return *this;
  129. }
  130.  
  131. Matrix&
  132. Matrix::insert (const RowVector& a, int r, int c)
  133. {
  134.   int a_len = a.length ();
  135.   if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ())
  136.     {
  137.       (*current_liboctave_error_handler) ("range error for insert");
  138.       return *this;
  139.     }
  140.  
  141.   for (int i = 0; i < a_len; i++)
  142.     elem (r, c+i) = a.elem (i);
  143.  
  144.   return *this;
  145. }
  146.  
  147. Matrix&
  148. Matrix::insert (const ColumnVector& a, int r, int c)
  149. {
  150.   int a_len = a.length ();
  151.   if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ())
  152.     {
  153.       (*current_liboctave_error_handler) ("range error for insert");
  154.       return *this;
  155.     }
  156.  
  157.   for (int i = 0; i < a_len; i++)
  158.     elem (r+i, c) = a.elem (i);
  159.  
  160.   return *this;
  161. }
  162.  
  163. Matrix&
  164. Matrix::insert (const DiagMatrix& a, int r, int c)
  165. {
  166.   if (r < 0 || r + a.rows () - 1 > rows ()
  167.       || c < 0 || c + a.cols () - 1 > cols ())
  168.     {
  169.       (*current_liboctave_error_handler) ("range error for insert");
  170.       return *this;
  171.     }
  172.  
  173.   for (int i = 0; i < a.length (); i++)
  174.     elem (r+i, c+i) = a.elem (i, i);
  175.  
  176.   return *this;
  177. }
  178.  
  179. Matrix&
  180. Matrix::fill (double val)
  181. {
  182.   int nr = rows ();
  183.   int nc = cols ();
  184.   if (nr > 0 && nc > 0)
  185.     for (int j = 0; j < nc; j++)
  186.       for (int i = 0; i < nr; i++)
  187.     elem (i, j) = val;
  188.  
  189.   return *this;
  190. }
  191.  
  192. Matrix&
  193. Matrix::fill (double val, int r1, int c1, int r2, int c2)
  194. {
  195.   int nr = rows ();
  196.   int nc = cols ();
  197.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  198.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  199.     {
  200.       (*current_liboctave_error_handler) ("range error for fill");
  201.       return *this;
  202.     }
  203.  
  204.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  205.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  206.  
  207.   for (int j = c1; j <= c2; j++)
  208.     for (int i = r1; i <= r2; i++)
  209.       elem (i, j) = val;
  210.  
  211.   return *this;
  212. }
  213.  
  214. Matrix
  215. Matrix::append (const Matrix& a) const
  216. {
  217.   int nr = rows ();
  218.   int nc = cols ();
  219.   if (nr != a.rows ())
  220.     {
  221.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  222.       return Matrix ();
  223.     }
  224.  
  225.   int nc_insert = nc;
  226.   Matrix retval (nr, nc + a.cols ());
  227.   retval.insert (*this, 0, 0);
  228.   retval.insert (a, 0, nc_insert);
  229.   return retval;
  230. }
  231.  
  232. Matrix
  233. Matrix::append (const RowVector& a) const
  234. {
  235.   int nr = rows ();
  236.   int nc = cols ();
  237.   if (nr != 1)
  238.     {
  239.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  240.       return Matrix ();
  241.     }
  242.  
  243.   int nc_insert = nc;
  244.   Matrix retval (nr, nc + a.length ());
  245.   retval.insert (*this, 0, 0);
  246.   retval.insert (a, 0, nc_insert);
  247.   return retval;
  248. }
  249.  
  250. Matrix
  251. Matrix::append (const ColumnVector& a) const
  252. {
  253.   int nr = rows ();
  254.   int nc = cols ();
  255.   if (nr != a.length ())
  256.     {
  257.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  258.       return Matrix ();
  259.     }
  260.  
  261.   int nc_insert = nc;
  262.   Matrix retval (nr, nc + 1);
  263.   retval.insert (*this, 0, 0);
  264.   retval.insert (a, 0, nc_insert);
  265.   return retval;
  266. }
  267.  
  268. Matrix
  269. Matrix::append (const DiagMatrix& a) const
  270. {
  271.   int nr = rows ();
  272.   int nc = cols ();
  273.   if (nr != a.rows ())
  274.     {
  275.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  276.       return *this;
  277.     }
  278.  
  279.   int nc_insert = nc;
  280.   Matrix retval (nr, nc + a.cols ());
  281.   retval.insert (*this, 0, 0);
  282.   retval.insert (a, 0, nc_insert);
  283.   return retval;
  284. }
  285.  
  286. Matrix
  287. Matrix::stack (const Matrix& a) const
  288. {
  289.   int nr = rows ();
  290.   int nc = cols ();
  291.   if (nc != a.cols ())
  292.     {
  293.       (*current_liboctave_error_handler)
  294.     ("column dimension mismatch for stack");
  295.       return Matrix ();
  296.     }
  297.  
  298.   int nr_insert = nr;
  299.   Matrix retval (nr + a.rows (), nc);
  300.   retval.insert (*this, 0, 0);
  301.   retval.insert (a, nr_insert, 0);
  302.   return retval;
  303. }
  304.  
  305. Matrix
  306. Matrix::stack (const RowVector& a) const
  307. {
  308.   int nr = rows ();
  309.   int nc = cols ();
  310.   if (nc != a.length ())
  311.     {
  312.       (*current_liboctave_error_handler)
  313.     ("column dimension mismatch for stack");
  314.       return Matrix ();
  315.     }
  316.  
  317.   int nr_insert = nr;
  318.   Matrix retval (nr + 1, nc);
  319.   retval.insert (*this, 0, 0);
  320.   retval.insert (a, nr_insert, 0);
  321.   return retval;
  322. }
  323.  
  324. Matrix
  325. Matrix::stack (const ColumnVector& a) const
  326. {
  327.   int nr = rows ();
  328.   int nc = cols ();
  329.   if (nc != 1)
  330.     {
  331.       (*current_liboctave_error_handler)
  332.     ("column dimension mismatch for stack");
  333.       return Matrix ();
  334.     }
  335.  
  336.   int nr_insert = nr;
  337.   Matrix retval (nr + a.length (), nc);
  338.   retval.insert (*this, 0, 0);
  339.   retval.insert (a, nr_insert, 0);
  340.   return retval;
  341. }
  342.  
  343. Matrix
  344. Matrix::stack (const DiagMatrix& a) const
  345. {
  346.   int nr = rows ();
  347.   int nc = cols ();
  348.   if (nc != a.cols ())
  349.     {
  350.       (*current_liboctave_error_handler)
  351.     ("column dimension mismatch for stack");
  352.       return Matrix ();
  353.     }
  354.  
  355.   int nr_insert = nr;
  356.   Matrix retval (nr + a.rows (), nc);
  357.   retval.insert (*this, 0, 0);
  358.   retval.insert (a, nr_insert, 0);
  359.   return retval;
  360. }
  361.  
  362. Matrix
  363. Matrix::transpose (void) const
  364. {
  365.   int nr = rows ();
  366.   int nc = cols ();
  367.   Matrix result (nc, nr);
  368.   if (length () > 0)
  369.     {
  370.       for (int j = 0; j < nc; j++)
  371.     for (int i = 0; i < nr; i++)
  372.       result.elem (j, i) = elem (i, j);
  373.     }
  374.   return result;
  375. }
  376.  
  377. Matrix
  378. real (const ComplexMatrix& a)
  379. {
  380.   int a_len = a.length ();
  381.   Matrix retval;
  382.   if (a_len > 0)
  383.     retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ());
  384.   return retval;
  385. }
  386.  
  387. Matrix
  388. imag (const ComplexMatrix& a)
  389. {
  390.   int a_len = a.length ();
  391.   Matrix retval;
  392.   if (a_len > 0)
  393.     retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ());
  394.   return retval;
  395. }
  396.  
  397. Matrix
  398. Matrix::extract (int r1, int c1, int r2, int c2) const
  399. {
  400.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  401.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  402.  
  403.   int new_r = r2 - r1 + 1;
  404.   int new_c = c2 - c1 + 1;
  405.  
  406.   Matrix result (new_r, new_c);
  407.  
  408.   for (int j = 0; j < new_c; j++)
  409.     for (int i = 0; i < new_r; i++)
  410.       result.elem (i, j) = elem (r1+i, c1+j);
  411.  
  412.   return result;
  413. }
  414.  
  415. // extract row or column i.
  416.  
  417. RowVector
  418. Matrix::row (int i) const
  419. {
  420.   int nc = cols ();
  421.   if (i < 0 || i >= rows ())
  422.     {
  423.       (*current_liboctave_error_handler) ("invalid row selection");
  424.       return RowVector ();
  425.     }
  426.  
  427.   RowVector retval (nc);
  428.   for (int j = 0; j < nc; j++)
  429.     retval.elem (j) = elem (i, j);
  430.  
  431.   return retval;
  432. }
  433.  
  434. RowVector
  435. Matrix::row (char *s) const
  436. {
  437.   if (! s)
  438.     {
  439.       (*current_liboctave_error_handler) ("invalid row selection");
  440.       return RowVector ();
  441.     }
  442.  
  443.   char c = *s;
  444.   if (c == 'f' || c == 'F')
  445.     return row (0);
  446.   else if (c == 'l' || c == 'L')
  447.     return row (rows () - 1);
  448.   else
  449.     {
  450.       (*current_liboctave_error_handler) ("invalid row selection");
  451.       return RowVector ();
  452.     }
  453. }
  454.  
  455. ColumnVector
  456. Matrix::column (int i) const
  457. {
  458.   int nr = rows ();
  459.   if (i < 0 || i >= cols ())
  460.     {
  461.       (*current_liboctave_error_handler) ("invalid column selection");
  462.       return ColumnVector ();
  463.     }
  464.  
  465.   ColumnVector retval (nr);
  466.   for (int j = 0; j < nr; j++)
  467.     retval.elem (j) = elem (j, i);
  468.  
  469.   return retval;
  470. }
  471.  
  472. ColumnVector
  473. Matrix::column (char *s) const
  474. {
  475.   if (! s)
  476.     {
  477.       (*current_liboctave_error_handler) ("invalid column selection");
  478.       return ColumnVector ();
  479.     }
  480.  
  481.   char c = *s;
  482.   if (c == 'f' || c == 'F')
  483.     return column (0);
  484.   else if (c == 'l' || c == 'L')
  485.     return column (cols () - 1);
  486.   else
  487.     {
  488.       (*current_liboctave_error_handler) ("invalid column selection");
  489.       return ColumnVector ();
  490.     }
  491. }
  492.  
  493. Matrix
  494. Matrix::inverse (void) const
  495. {
  496.   int info;
  497.   double rcond;
  498.   return inverse (info, rcond);
  499. }
  500.  
  501. Matrix
  502. Matrix::inverse (int& info) const
  503. {
  504.   double rcond;
  505.   return inverse (info, rcond);
  506. }
  507.  
  508. Matrix
  509. Matrix::inverse (int& info, double& rcond) const
  510. {
  511.   int nr = rows ();
  512.   int nc = cols ();
  513.   int len = length ();
  514.   if (nr != nc || nr == 0 || nc == 0)
  515.     {
  516.       (*current_liboctave_error_handler) ("inverse requires square matrix");
  517.       return Matrix ();
  518.     }
  519.  
  520.   info = 0;
  521.  
  522.   int *ipvt = new int [nr];
  523.   double *z = new double [nr];
  524.   double *tmp_data = dup (data (), len);
  525.  
  526.   F77_FCN (dgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z);
  527.  
  528.   volatile double tmp_rcond = rcond;
  529.   if (tmp_rcond + 1.0 == 1.0)
  530.     {
  531.       info = -1;
  532.       copy (tmp_data, data (), len);  // Restore matrix contents.
  533.     }
  534.   else
  535.     {
  536.       int job = 1;
  537.       double dummy;
  538.  
  539.       F77_FCN (dgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job);
  540.     }
  541.  
  542.   delete [] ipvt;
  543.   delete [] z;
  544.  
  545.   return Matrix (tmp_data, nr, nc);
  546. }
  547.  
  548. Matrix
  549. Matrix::pseudo_inverse (double tol)
  550. {
  551.   SVD result (*this);
  552.  
  553.   DiagMatrix S = result.singular_values ();
  554.   Matrix U = result.left_singular_matrix ();
  555.   Matrix V = result.right_singular_matrix ();
  556.  
  557.   ColumnVector sigma = S.diag ();
  558.  
  559.   int r = sigma.length () - 1;
  560.   int nr = rows ();
  561.   int nc = cols ();
  562.  
  563.   if (tol <= 0.0)
  564.     {
  565.       if (nr > nc)
  566.     tol = nr * sigma.elem (0) * DBL_EPSILON;
  567.       else
  568.     tol = nc * sigma.elem (0) * DBL_EPSILON;
  569.     }
  570.  
  571.   while (r >= 0 && sigma.elem (r) < tol)
  572.     r--;
  573.  
  574.   if (r < 0)
  575.     return Matrix (nc, nr, 0.0);
  576.   else
  577.     {
  578.       Matrix Ur = U.extract (0, 0, nr-1, r);
  579.       DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
  580.       Matrix Vr = V.extract (0, 0, nc-1, r);
  581.       return Vr * D * Ur.transpose ();
  582.     }
  583. }
  584.  
  585. ComplexMatrix
  586. Matrix::fourier (void) const
  587. {
  588.   int nr = rows ();
  589.   int nc = cols ();
  590.   int npts, nsamples;
  591.   if (nr == 1 || nc == 1)
  592.     {
  593.       npts = nr > nc ? nr : nc;
  594.       nsamples = 1;
  595.     }
  596.   else
  597.     {
  598.       npts = nr;
  599.       nsamples = nc;
  600.     }
  601.  
  602.   int nn = 4*npts+15;
  603.   Complex *wsave = new Complex [nn];
  604.   Complex *tmp_data = make_complex (data (), length ());
  605.  
  606.   F77_FCN (cffti) (&npts, wsave);
  607.  
  608.   for (int j = 0; j < nsamples; j++)
  609.     F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave);
  610.  
  611.   delete [] wsave;
  612.  
  613.   return ComplexMatrix (tmp_data, nr, nc);
  614. }
  615.  
  616. ComplexMatrix
  617. Matrix::ifourier (void) const
  618. {
  619.   int nr = rows ();
  620.   int nc = cols ();
  621.   int npts, nsamples;
  622.   if (nr == 1 || nc == 1)
  623.     {
  624.       npts = nr > nc ? nr : nc;
  625.       nsamples = 1;
  626.     }
  627.   else
  628.     {
  629.       npts = nr;
  630.       nsamples = nc;
  631.     }
  632.  
  633.   int nn = 4*npts+15;
  634.   Complex *wsave = new Complex [nn];
  635.   Complex *tmp_data = make_complex (data (), length ());
  636.  
  637.   F77_FCN (cffti) (&npts, wsave);
  638.  
  639.   for (int j = 0; j < nsamples; j++)
  640.     F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave);
  641.  
  642.   for (j = 0; j < npts*nsamples; j++)
  643.     tmp_data[j] = tmp_data[j] / (double) npts;
  644.  
  645.   delete [] wsave;
  646.  
  647.   return ComplexMatrix (tmp_data, nr, nc);
  648. }
  649.  
  650. ComplexMatrix
  651. Matrix::fourier2d (void) const
  652. {
  653.   int nr = rows ();
  654.   int nc = cols ();
  655.   int npts, nsamples;
  656.   if (nr == 1 || nc == 1)
  657.     {
  658.       npts = nr > nc ? nr : nc;
  659.       nsamples = 1;
  660.     }
  661.   else
  662.     {
  663.       npts = nr;
  664.       nsamples = nc;
  665.     }
  666.  
  667.   int nn = 4*npts+15;
  668.   Complex *wsave = new Complex [nn];
  669.   Complex *tmp_data = make_complex (data (), length ());
  670.  
  671.   F77_FCN (cffti) (&npts, wsave);
  672.  
  673.   for (int j = 0; j < nsamples; j++)
  674.     F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave);
  675.  
  676.   delete [] wsave;
  677.  
  678.   npts = nc;
  679.   nsamples = nr;
  680.   nn = 4*npts+15;
  681.   wsave = new Complex [nn];
  682.   Complex *row = new Complex[npts];
  683.  
  684.   F77_FCN (cffti) (&npts, wsave);
  685.  
  686.   for (j = 0; j < nsamples; j++)
  687.     {
  688.       for (int i = 0; i < npts; i++)
  689.     row[i] = tmp_data[i*nr + j];
  690.  
  691.       F77_FCN (cfftf) (&npts, row, wsave);
  692.  
  693.       for (i = 0; i < npts; i++)
  694.     tmp_data[i*nr + j] = row[i];
  695.     }
  696.  
  697.   delete [] wsave;
  698.   delete [] row;
  699.  
  700.   return ComplexMatrix (tmp_data, nr, nc);
  701. }
  702.  
  703. ComplexMatrix
  704. Matrix::ifourier2d (void) const
  705. {
  706.   int nr = rows ();
  707.   int nc = cols ();
  708.   int npts, nsamples;
  709.   if (nr == 1 || nc == 1)
  710.     {
  711.       npts = nr > nc ? nr : nc;
  712.       nsamples = 1;
  713.     }
  714.   else
  715.     {
  716.       npts = nr;
  717.       nsamples = nc;
  718.     }
  719.  
  720.   int nn = 4*npts+15;
  721.   Complex *wsave = new Complex [nn];
  722.   Complex *tmp_data = make_complex (data (), length ());
  723.  
  724.   F77_FCN (cffti) (&npts, wsave);
  725.  
  726.   for (int j = 0; j < nsamples; j++)
  727.     F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave);
  728.  
  729.   delete [] wsave;
  730.  
  731.   for (j = 0; j < npts*nsamples; j++)
  732.     tmp_data[j] = tmp_data[j] / (double) npts;
  733.  
  734.   npts = nc;
  735.   nsamples = nr;
  736.   nn = 4*npts+15;
  737.   wsave = new Complex [nn];
  738.   Complex *row = new Complex[npts];
  739.  
  740.   F77_FCN (cffti) (&npts, wsave);
  741.  
  742.   for (j = 0; j < nsamples; j++)
  743.     {
  744.       for (int i = 0; i < npts; i++)
  745.     row[i] = tmp_data[i*nr + j];
  746.  
  747.       F77_FCN (cfftb) (&npts, row, wsave);
  748.  
  749.       for (i = 0; i < npts; i++)
  750.     tmp_data[i*nr + j] = row[i] / (double) npts;
  751.     }
  752.  
  753.   delete [] wsave;
  754.   delete [] row;
  755.  
  756.   return ComplexMatrix (tmp_data, nr, nc);
  757. }
  758.  
  759. DET
  760. Matrix::determinant (void) const
  761. {
  762.   int info;
  763.   double rcond;
  764.   return determinant (info, rcond);
  765. }
  766.  
  767. DET
  768. Matrix::determinant (int& info) const
  769. {
  770.   double rcond;
  771.   return determinant (info, rcond);
  772. }
  773.  
  774. DET
  775. Matrix::determinant (int& info, double& rcond) const
  776. {
  777.   DET retval;
  778.  
  779.   int nr = rows ();
  780.   int nc = cols ();
  781.  
  782.   if (nr == 0 || nc == 0)
  783.     {
  784.       double d[2];
  785.       d[0] = 1.0;
  786.       d[1] = 0.0;
  787.       retval = DET (d);
  788.     }
  789.   else
  790.     {
  791.       info = 0;
  792.       int *ipvt = new int [nr];
  793.  
  794.       double *z = new double [nr];
  795.       double *tmp_data = dup (data (), length ());
  796.  
  797.       F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
  798.  
  799.       volatile double tmp_rcond = rcond;
  800.       if (tmp_rcond + 1.0 == 1.0)
  801.     {
  802.       info = -1;
  803.       retval = DET ();
  804.     }
  805.       else
  806.     {
  807.       int job = 10;
  808.       double d[2];
  809.       F77_FCN (dgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job);
  810.       retval = DET (d);
  811.     }
  812.  
  813.       delete [] tmp_data;
  814.       delete [] ipvt;
  815.       delete [] z;
  816.     }
  817.  
  818.   return retval;
  819. }
  820.  
  821. Matrix
  822. Matrix::solve (const Matrix& b) const
  823. {
  824.   int info;
  825.   double rcond;
  826.   return solve (b, info, rcond);
  827. }
  828.  
  829. Matrix
  830. Matrix::solve (const Matrix& b, int& info) const
  831. {
  832.   double rcond;
  833.   return solve (b, info, rcond);
  834. }
  835.  
  836. Matrix
  837. Matrix::solve (const Matrix& b, int& info, double& rcond) const
  838. {
  839.   Matrix retval;
  840.  
  841.   int nr = rows ();
  842.   int nc = cols ();
  843.   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
  844.     {
  845.       (*current_liboctave_error_handler)
  846.     ("matrix dimension mismatch solution of linear equations");
  847.       return Matrix ();
  848.     }
  849.  
  850.   info = 0;
  851.   int *ipvt = new int [nr];
  852.  
  853.   double *z = new double [nr];
  854.   double *tmp_data = dup (data (), length ());
  855.  
  856.   F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
  857.  
  858.   volatile double tmp_rcond = rcond;
  859.   if (tmp_rcond + 1.0 == 1.0)
  860.     {
  861.       info = -2;
  862.     }
  863.   else
  864.     {
  865.       int job = 0;
  866.  
  867.       double *result = dup (b.data (), b.length ());
  868.  
  869.       int b_nc = b.cols ();
  870.       for (int j = 0; j < b_nc; j++)
  871.     F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job);
  872.  
  873.       retval = Matrix (result, b.rows (), b_nc);
  874.     }
  875.  
  876.   delete [] tmp_data;
  877.   delete [] ipvt;
  878.   delete [] z;
  879.  
  880.   return retval;
  881. }
  882.  
  883. ComplexMatrix
  884. Matrix::solve (const ComplexMatrix& b) const
  885. {
  886.   ComplexMatrix tmp (*this);
  887.   return tmp.solve (b);
  888. }
  889.  
  890. ComplexMatrix
  891. Matrix::solve (const ComplexMatrix& b, int& info) const
  892. {
  893.   ComplexMatrix tmp (*this);
  894.   return tmp.solve (b, info);
  895. }
  896.  
  897. ComplexMatrix
  898. Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
  899. {
  900.   ComplexMatrix tmp (*this);
  901.   return tmp.solve (b, info, rcond);
  902. }
  903.  
  904. ColumnVector
  905. Matrix::solve (const ColumnVector& b) const
  906. {
  907.   int info; double rcond;
  908.   return solve (b, info, rcond);
  909. }
  910.  
  911. ColumnVector
  912. Matrix::solve (const ColumnVector& b, int& info) const
  913. {
  914.   double rcond;
  915.   return solve (b, info, rcond);
  916. }
  917.  
  918. ColumnVector
  919. Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
  920. {
  921.   ColumnVector retval;
  922.  
  923.   int nr = rows ();
  924.   int nc = cols ();
  925.   if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
  926.     {
  927.       (*current_liboctave_error_handler)
  928.     ("matrix dimension mismatch solution of linear equations");
  929.       return ColumnVector ();
  930.     }
  931.  
  932.   info = 0;
  933.   int *ipvt = new int [nr];
  934.  
  935.   double *z = new double [nr];
  936.   double *tmp_data = dup (data (), length ());
  937.  
  938.   F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
  939.  
  940.   volatile double tmp_rcond = rcond;
  941.   if (tmp_rcond + 1.0 == 1.0)
  942.     {
  943.       info = -2;
  944.     }
  945.   else
  946.     {
  947.       int job = 0;
  948.  
  949.       int b_len = b.length ();
  950.  
  951.       double *result = dup (b.data (), b_len);
  952.  
  953.       F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, result, &job);
  954.  
  955.       retval = ColumnVector (result, b_len);
  956.     }
  957.  
  958.   delete [] tmp_data;
  959.   delete [] ipvt;
  960.   delete [] z;
  961.  
  962.   return retval;
  963. }
  964.  
  965. ComplexColumnVector
  966. Matrix::solve (const ComplexColumnVector& b) const
  967. {
  968.   ComplexMatrix tmp (*this);
  969.   return tmp.solve (b);
  970. }
  971.  
  972. ComplexColumnVector
  973. Matrix::solve (const ComplexColumnVector& b, int& info) const
  974. {
  975.   ComplexMatrix tmp (*this);
  976.   return tmp.solve (b, info);
  977. }
  978.  
  979. ComplexColumnVector
  980. Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
  981. {
  982.   ComplexMatrix tmp (*this);
  983.   return tmp.solve (b, info, rcond);
  984. }
  985.  
  986. Matrix
  987. Matrix::lssolve (const Matrix& b) const
  988. {
  989.   int info;
  990.   int rank;
  991.   return lssolve (b, info, rank);
  992. }
  993.  
  994. Matrix
  995. Matrix::lssolve (const Matrix& b, int& info) const
  996. {
  997.   int rank;
  998.   return lssolve (b, info, rank);
  999. }
  1000.  
  1001. Matrix
  1002. Matrix::lssolve (const Matrix& b, int& info, int& rank) const
  1003. {
  1004.   int nrhs = b.cols ();
  1005.  
  1006.   int m = rows ();
  1007.   int n = cols ();
  1008.  
  1009.   if (m == 0 || n == 0 || m != b.rows ())
  1010.     {
  1011.       (*current_liboctave_error_handler)
  1012.     ("matrix dimension mismatch in solution of least squares problem");
  1013.       return Matrix ();
  1014.     }
  1015.  
  1016.   double *tmp_data = dup (data (), length ());
  1017.  
  1018.   int nrr = m > n ? m : n;
  1019.   Matrix result (nrr, nrhs);
  1020.  
  1021.   int i, j;
  1022.   for (j = 0; j < nrhs; j++)
  1023.     for (i = 0; i < m; i++)
  1024.       result.elem (i, j) = b.elem (i, j);
  1025.  
  1026.   double *presult = result.fortran_vec ();
  1027.  
  1028.   int len_s = m < n ? m : n;
  1029.   double *s = new double [len_s];
  1030.   double rcond = -1.0;
  1031.   int lwork;
  1032.   if (m < n)
  1033.     lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
  1034.   else
  1035.     lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
  1036.  
  1037.   double *work = new double [lwork];
  1038.  
  1039.   F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
  1040.             &rcond, &rank, work, &lwork, &info);
  1041.  
  1042.   Matrix retval (n, nrhs);
  1043.   for (j = 0; j < nrhs; j++)
  1044.     for (i = 0; i < n; i++)
  1045.       retval.elem (i, j) = result.elem (i, j);
  1046.  
  1047.   delete [] tmp_data;
  1048.   delete [] s;
  1049.   delete [] work;
  1050.  
  1051.   return retval;
  1052. }
  1053.  
  1054. ComplexMatrix
  1055. Matrix::lssolve (const ComplexMatrix& b) const
  1056. {
  1057.   ComplexMatrix tmp (*this);
  1058.   return tmp.lssolve (b);
  1059. }
  1060.  
  1061. ComplexMatrix
  1062. Matrix::lssolve (const ComplexMatrix& b, int& info) const
  1063. {
  1064.   ComplexMatrix tmp (*this);
  1065.   return tmp.lssolve (b);
  1066. }
  1067.  
  1068. ComplexMatrix
  1069. Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
  1070. {
  1071.   ComplexMatrix tmp (*this);
  1072.   return tmp.lssolve (b);
  1073. }
  1074.  
  1075. ColumnVector
  1076. Matrix::lssolve (const ColumnVector& b) const
  1077. {
  1078.   int info;
  1079.   int rank; return lssolve (b, info, rank);
  1080. }
  1081.  
  1082. ColumnVector
  1083. Matrix::lssolve (const ColumnVector& b, int& info) const
  1084. {
  1085.   int rank;
  1086.   return lssolve (b, info, rank);
  1087. }
  1088.  
  1089. ColumnVector
  1090. Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const
  1091. {
  1092.   int nrhs = 1;
  1093.  
  1094.   int m = rows ();
  1095.   int n = cols ();
  1096.  
  1097.   if (m == 0 || n == 0 || m != b.length ())
  1098.     {
  1099.       (*current_liboctave_error_handler)
  1100.     ("matrix dimension mismatch in solution of least squares problem");
  1101.       return ColumnVector ();
  1102.     }
  1103.  
  1104.   double *tmp_data = dup (data (), length ());
  1105.  
  1106.   int nrr = m > n ? m : n;
  1107.   ColumnVector result (nrr);
  1108.  
  1109.   int i;
  1110.   for (i = 0; i < m; i++)
  1111.     result.elem (i) = b.elem (i);
  1112.  
  1113.   double *presult = result.fortran_vec ();
  1114.  
  1115.   int len_s = m < n ? m : n;
  1116.   double *s = new double [len_s];
  1117.   double rcond = -1.0;
  1118.   int lwork;
  1119.   if (m < n)
  1120.     lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
  1121.   else
  1122.     lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
  1123.  
  1124.   double *work = new double [lwork];
  1125.  
  1126.   F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
  1127.             &rcond, &rank, work, &lwork, &info);
  1128.  
  1129.   ColumnVector retval (n);
  1130.   for (i = 0; i < n; i++)
  1131.     retval.elem (i) = result.elem (i);
  1132.  
  1133.   delete [] tmp_data;
  1134.   delete [] s;
  1135.   delete [] work;
  1136.  
  1137.   return retval;
  1138. }
  1139.  
  1140. ComplexColumnVector
  1141. Matrix::lssolve (const ComplexColumnVector& b) const
  1142. {
  1143.   ComplexMatrix tmp (*this);
  1144.   return tmp.lssolve (b);
  1145. }
  1146.  
  1147. ComplexColumnVector
  1148. Matrix::lssolve (const ComplexColumnVector& b, int& info) const
  1149. {
  1150.   ComplexMatrix tmp (*this);
  1151.   return tmp.lssolve (b, info);
  1152. }
  1153.  
  1154. ComplexColumnVector
  1155. Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
  1156. {
  1157.   ComplexMatrix tmp (*this);
  1158.   return tmp.lssolve (b, info, rank);
  1159. }
  1160.  
  1161. Matrix&
  1162. Matrix::operator += (const Matrix& a)
  1163. {
  1164.   int nr = rows ();
  1165.   int nc = cols ();
  1166.   if (nr != a.rows () || nc != a.cols ())
  1167.     {
  1168.       (*current_liboctave_error_handler)
  1169.     ("nonconformant matrix += operation attempted");
  1170.       return *this;
  1171.     }
  1172.  
  1173.   if (nr == 0 || nc == 0)
  1174.     return *this;
  1175.  
  1176.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1177.  
  1178.   add2 (d, a.data (), length ());
  1179.  
  1180.   return *this;
  1181. }
  1182.  
  1183. Matrix&
  1184. Matrix::operator -= (const Matrix& a)
  1185. {
  1186.   int nr = rows ();
  1187.   int nc = cols ();
  1188.   if (nr != a.rows () || nc != a.cols ())
  1189.     {
  1190.       (*current_liboctave_error_handler)
  1191.     ("nonconformant matrix -= operation attempted");
  1192.       return *this;
  1193.     }
  1194.  
  1195.   if (nr == 0 || nc == 0)
  1196.     return *this;
  1197.  
  1198.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1199.  
  1200.   subtract2 (d, a.data (), length ());
  1201.  
  1202.   return *this;
  1203. }
  1204.  
  1205. Matrix&
  1206. Matrix::operator += (const DiagMatrix& a)
  1207. {
  1208.   if (rows () != a.rows () || cols () != a.cols ())
  1209.     {
  1210.       (*current_liboctave_error_handler)
  1211.     ("nonconformant matrix += operation attempted");
  1212.       return *this;
  1213.     }
  1214.  
  1215.   for (int i = 0; i < a.length (); i++)
  1216.     elem (i, i) += a.elem (i, i);
  1217.  
  1218.   return *this;
  1219. }
  1220.  
  1221. Matrix&
  1222. Matrix::operator -= (const DiagMatrix& a)
  1223. {
  1224.   if (rows () != a.rows () || cols () != a.cols ())
  1225.     {
  1226.       (*current_liboctave_error_handler)
  1227.     ("nonconformant matrix += operation attempted");
  1228.       return *this;
  1229.     }
  1230.  
  1231.   for (int i = 0; i < a.length (); i++)
  1232.     elem (i, i) -= a.elem (i, i);
  1233.  
  1234.   return *this;
  1235. }
  1236.  
  1237. // unary operations
  1238.  
  1239. Matrix
  1240. Matrix::operator ! (void) const
  1241. {
  1242.   int nr = rows ();
  1243.   int nc = cols ();
  1244.  
  1245.   Matrix b (nr, nc);
  1246.  
  1247.   for (int j = 0; j < nc; j++)
  1248.     for (int i = 0; i < nr; i++)
  1249.       b.elem (i, j) = ! elem (i, j);
  1250.  
  1251.   return b;
  1252. }
  1253.  
  1254. // column vector by row vector -> matrix operations
  1255.  
  1256. Matrix
  1257. operator * (const ColumnVector& v, const RowVector& a)
  1258. {
  1259.   int len = v.length ();
  1260.   int a_len = a.length ();
  1261.   if (len != a_len)
  1262.     {
  1263.       (*current_liboctave_error_handler)
  1264.     ("nonconformant vector multiplication attempted");
  1265.       return Matrix ();
  1266.     }
  1267.  
  1268.   if (len == 0)
  1269.     return Matrix (len, len, 0.0);
  1270.  
  1271.   char transa = 'N';
  1272.   char transb = 'N';
  1273.   double alpha = 1.0;
  1274.   double beta  = 0.0;
  1275.   int anr = 1;
  1276.  
  1277.   double *c = new double [len * a_len];
  1278.  
  1279.   F77_FCN (dgemm) (&transa, &transb, &len, &a_len, &anr, &alpha,
  1280.            v.data (), &len, a.data (), &anr, &beta, c, &len,
  1281.            1L, 1L); 
  1282.  
  1283.   return Matrix (c, len, a_len);
  1284. }
  1285.  
  1286. // diagonal matrix by scalar -> matrix operations
  1287.  
  1288. Matrix
  1289. operator + (const DiagMatrix& a, double s)
  1290. {
  1291.   Matrix tmp (a.rows (), a.cols (), s);
  1292.   return a + tmp;
  1293. }
  1294.  
  1295. Matrix
  1296. operator - (const DiagMatrix& a, double s)
  1297. {
  1298.   Matrix tmp (a.rows (), a.cols (), -s);
  1299.   return a + tmp;
  1300. }
  1301.  
  1302. // scalar by diagonal matrix -> matrix operations
  1303.  
  1304. Matrix
  1305. operator + (double s, const DiagMatrix& a)
  1306. {
  1307.   Matrix tmp (a.rows (), a.cols (), s);
  1308.   return tmp + a;
  1309. }
  1310.  
  1311. Matrix
  1312. operator - (double s, const DiagMatrix& a)
  1313. {
  1314.   Matrix tmp (a.rows (), a.cols (), s);
  1315.   return tmp - a;
  1316. }
  1317.  
  1318. // matrix by diagonal matrix -> matrix operations
  1319.  
  1320. Matrix
  1321. operator + (const Matrix& m, const DiagMatrix& a)
  1322. {
  1323.   int nr = m.rows ();
  1324.   int nc = m.cols ();
  1325.   if (nr != a.rows () || nc != a.cols ())
  1326.     {
  1327.       (*current_liboctave_error_handler)
  1328.     ("nonconformant matrix addition attempted");
  1329.       return Matrix ();
  1330.     }
  1331.  
  1332.   if (nr == 0 || nc == 0)
  1333.     return Matrix (nr, nc);
  1334.  
  1335.   Matrix result (m);
  1336.   int a_len = a.length ();
  1337.   for (int i = 0; i < a_len; i++)
  1338.     result.elem (i, i) += a.elem (i, i);
  1339.  
  1340.   return result;
  1341. }
  1342.  
  1343. Matrix
  1344. operator - (const Matrix& m, const DiagMatrix& a)
  1345. {
  1346.   int nr = m.rows ();
  1347.   int nc = m.cols ();
  1348.   if (nr != a.rows () || nc != a.cols ())
  1349.     {
  1350.       (*current_liboctave_error_handler)
  1351.     ("nonconformant matrix subtraction attempted");
  1352.       return Matrix ();
  1353.     }
  1354.  
  1355.   if (nr == 0 || nc == 0)
  1356.     return Matrix (nr, nc);
  1357.  
  1358.   Matrix result (m);
  1359.   int a_len = a.length ();
  1360.   for (int i = 0; i < a_len; i++)
  1361.     result.elem (i, i) -= a.elem (i, i);
  1362.  
  1363.   return result;
  1364. }
  1365.  
  1366. Matrix
  1367. operator * (const Matrix& m, const DiagMatrix& a)
  1368. {
  1369.   int nr = m.rows ();
  1370.   int nc = m.cols ();
  1371.   int a_nr = a.rows ();
  1372.   int a_nc = a.cols ();
  1373.   if (nc != a_nr)
  1374.     {
  1375.       (*current_liboctave_error_handler)
  1376.     ("nonconformant matrix multiplication attempted");
  1377.       return Matrix ();
  1378.     }
  1379.  
  1380.   if (nr == 0 || nc == 0 || a_nc == 0)
  1381.     return Matrix (nr, a_nc, 0.0);
  1382.  
  1383.   double *c = new double [nr*a_nc];
  1384.   double *ctmp = 0;
  1385.  
  1386.   int a_len = a.length ();
  1387.   for (int j = 0; j < a_len; j++)
  1388.     {
  1389.       int idx = j * nr;
  1390.       ctmp = c + idx;
  1391.       if (a.elem (j, j) == 1.0)
  1392.     {
  1393.       for (int i = 0; i < nr; i++)
  1394.         ctmp[i] = m.elem (i, j);
  1395.     }
  1396.       else if (a.elem (j, j) == 0.0)
  1397.     {
  1398.       for (int i = 0; i < nr; i++)
  1399.         ctmp[i] = 0.0;
  1400.     }
  1401.       else
  1402.     {
  1403.       for (int i = 0; i < nr; i++)
  1404.         ctmp[i] = a.elem (j, j) * m.elem (i, j);
  1405.     }
  1406.     }
  1407.  
  1408.   if (a_nr < a_nc)
  1409.     {
  1410.       for (int i = nr * nc; i < nr * a_nc; i++)
  1411.     ctmp[i] = 0.0;
  1412.     }
  1413.  
  1414.   return Matrix (c, nr, a_nc);
  1415. }
  1416.  
  1417. // diagonal matrix by matrix -> matrix operations
  1418.  
  1419. Matrix
  1420. operator + (const DiagMatrix& m, const Matrix& a)
  1421. {
  1422.   int nr = m.rows ();
  1423.   int nc = m.cols ();
  1424.   if (nr != a.rows () || nc != a.cols ())
  1425.     {
  1426.       (*current_liboctave_error_handler)
  1427.     ("nonconformant matrix addition attempted");
  1428.       return Matrix ();
  1429.     }
  1430.  
  1431.   if (nr == 0 || nc == 0)
  1432.     return Matrix (nr, nc);
  1433.  
  1434.   Matrix result (a);
  1435.   for (int i = 0; i < m.length (); i++)
  1436.     result.elem (i, i) += m.elem (i, i);
  1437.  
  1438.   return result;
  1439. }
  1440.  
  1441. Matrix
  1442. operator - (const DiagMatrix& m, const Matrix& a)
  1443. {
  1444.   int nr = m.rows ();
  1445.   int nc = m.cols ();
  1446.   if (nr != a.rows () || nc != a.cols ())
  1447.     {
  1448.       (*current_liboctave_error_handler)
  1449.     ("nonconformant matrix subtraction attempted");
  1450.       return Matrix ();
  1451.     }
  1452.  
  1453.   if (nr == 0 || nc == 0)
  1454.     return Matrix (nr, nc);
  1455.  
  1456.   Matrix result (-a);
  1457.   for (int i = 0; i < m.length (); i++)
  1458.     result.elem (i, i) += m.elem (i, i);
  1459.  
  1460.   return result;
  1461. }
  1462.  
  1463. Matrix
  1464. operator * (const DiagMatrix& m, const Matrix& a)
  1465. {
  1466.   int nr = m.rows ();
  1467.   int nc = m.cols ();
  1468.   int a_nr = a.rows ();
  1469.   int a_nc = a.cols ();
  1470.   if (nc != a_nr)
  1471.     {
  1472.       (*current_liboctave_error_handler)
  1473.     ("nonconformant matrix multiplication attempted");
  1474.       return Matrix ();
  1475.     }
  1476.  
  1477.   if (nr == 0 || nc == 0 || a_nc == 0)
  1478.     return Matrix (nr, a_nc, 0.0);
  1479.  
  1480.   Matrix c (nr, a_nc);
  1481.  
  1482.   for (int i = 0; i < m.length (); i++)
  1483.     {
  1484.       if (m.elem (i, i) == 1.0)
  1485.     {
  1486.       for (int j = 0; j < a_nc; j++)
  1487.         c.elem (i, j) = a.elem (i, j);
  1488.     }
  1489.       else if (m.elem (i, i) == 0.0)
  1490.     {
  1491.       for (int j = 0; j < a_nc; j++)
  1492.         c.elem (i, j) = 0.0;
  1493.     }
  1494.       else
  1495.     {
  1496.       for (int j = 0; j < a_nc; j++)
  1497.         c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
  1498.     }
  1499.     }
  1500.  
  1501.   if (nr > nc)
  1502.     {
  1503.       for (int j = 0; j < a_nc; j++)
  1504.     for (int i = a_nr; i < nr; i++)
  1505.       c.elem (i, j) = 0.0;
  1506.     }
  1507.  
  1508.   return c;
  1509. }
  1510.  
  1511. // matrix by matrix -> matrix operations
  1512.  
  1513. Matrix
  1514. operator * (const Matrix& m, const Matrix& a)
  1515. {
  1516.   int nr = m.rows ();
  1517.   int nc = m.cols ();
  1518.   int a_nr = a.rows ();
  1519.   int a_nc = a.cols ();
  1520.   if (nc != a_nr)
  1521.     {
  1522.       (*current_liboctave_error_handler)
  1523.     ("nonconformant matrix multiplication attempted");
  1524.       return Matrix ();
  1525.     }
  1526.  
  1527.   if (nr == 0 || nc == 0 || a_nc == 0)
  1528.     return Matrix (nr, a_nc, 0.0);
  1529.  
  1530.   char trans  = 'N';
  1531.   char transa = 'N';
  1532.  
  1533.   int ld  = nr;
  1534.   int lda = a_nr;
  1535.  
  1536.   double alpha = 1.0;
  1537.   double beta  = 0.0;
  1538.  
  1539.   double *c = new double [nr*a_nc];
  1540.  
  1541.   F77_FCN (dgemm) (&trans, &transa, &nr, &a_nc, &nc, &alpha, m.data (),
  1542.            &ld, a.data (), &lda, &beta, c, &nr, 1L, 1L);
  1543.  
  1544.   return Matrix (c, nr, a_nc);
  1545. }
  1546.  
  1547. // other operations.
  1548.  
  1549. Matrix
  1550. map (d_d_Mapper f, const Matrix& a)
  1551. {
  1552.   Matrix b (a);
  1553.   b.map (f);
  1554.   return b;
  1555. }
  1556.  
  1557. Matrix
  1558. map (d_c_Mapper f, const ComplexMatrix& a)
  1559. {
  1560.   int a_nc = a.cols ();
  1561.   int a_nr = a.rows ();
  1562.   Matrix b (a_nr, a_nc);
  1563.   for (int j = 0; j < a_nc; j++)
  1564.     for (int i = 0; i < a_nr; i++)
  1565.       b.elem (i, j) = f (a.elem (i, j));
  1566.   return b;
  1567. }
  1568.  
  1569. void
  1570. Matrix::map (d_d_Mapper f)
  1571. {
  1572.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1573.  
  1574.   for (int i = 0; i < length (); i++)
  1575.     d[i] = f (d[i]);
  1576. }
  1577.  
  1578. // XXX FIXME XXX Do these really belong here?  They should maybe be
  1579. // cleaned up a bit, no?  What about corresponding functions for the
  1580. // Vectors?
  1581.  
  1582. Matrix
  1583. Matrix::all (void) const
  1584. {
  1585.   int nr = rows ();
  1586.   int nc = cols ();
  1587.   Matrix retval;
  1588.   if (nr > 0 && nc > 0)
  1589.     {
  1590.       if (nr == 1)
  1591.     {
  1592.       retval.resize (1, 1);
  1593.       retval.elem (0, 0) = 1.0;
  1594.       for (int j = 0; j < nc; j++)
  1595.         {
  1596.           if (elem (0, j) == 0.0)
  1597.         {
  1598.           retval.elem (0, 0) = 0.0;
  1599.           break;
  1600.         }
  1601.         }
  1602.     }
  1603.       else if (nc == 1)
  1604.     {
  1605.       retval.resize (1, 1);
  1606.       retval.elem (0, 0) = 1.0;
  1607.       for (int i = 0; i < nr; i++)
  1608.         {
  1609.           if (elem (i, 0) == 0.0)
  1610.         {
  1611.           retval.elem (0, 0) = 0.0;
  1612.           break;
  1613.         }
  1614.         }
  1615.     }
  1616.       else
  1617.     {
  1618.       retval.resize (1, nc);
  1619.       for (int j = 0; j < nc; j++)
  1620.         {
  1621.           retval.elem (0, j) = 1.0;
  1622.           for (int i = 0; i < nr; i++)
  1623.         {
  1624.           if (elem (i, j) == 0.0)
  1625.             {
  1626.               retval.elem (0, j) = 0.0;
  1627.               break;
  1628.             }
  1629.         }
  1630.         }
  1631.     }
  1632.     }
  1633.   return retval;
  1634. }
  1635.  
  1636. Matrix
  1637. Matrix::any (void) const
  1638. {
  1639.   int nr = rows ();
  1640.   int nc = cols ();
  1641.   Matrix retval;
  1642.   if (nr > 0 && nc > 0)
  1643.     {
  1644.       if (nr == 1)
  1645.     {
  1646.       retval.resize (1, 1);
  1647.       retval.elem (0, 0) = 0.0;
  1648.       for (int j = 0; j < nc; j++)
  1649.         {
  1650.           if (elem (0, j) != 0.0)
  1651.         {
  1652.           retval.elem (0, 0) = 1.0;
  1653.           break;
  1654.         }
  1655.         }
  1656.     }
  1657.       else if (nc == 1)
  1658.     {
  1659.       retval.resize (1, 1);
  1660.       retval.elem (0, 0) = 0.0;
  1661.       for (int i = 0; i < nr; i++)
  1662.         {
  1663.           if (elem (i, 0) != 0.0)
  1664.         {
  1665.           retval.elem (0, 0) = 1.0;
  1666.           break;
  1667.         }
  1668.         }
  1669.     }
  1670.       else
  1671.     {
  1672.       retval.resize (1, nc);
  1673.       for (int j = 0; j < nc; j++)
  1674.         {
  1675.           retval.elem (0, j) = 0.0;
  1676.           for (int i = 0; i < nr; i++)
  1677.         {
  1678.           if (elem (i, j) != 0.0)
  1679.             {
  1680.               retval.elem (0, j) = 1.0;
  1681.               break;
  1682.             }
  1683.         }
  1684.         }
  1685.     }
  1686.     }
  1687.   return retval;
  1688. }
  1689.  
  1690. Matrix
  1691. Matrix::cumprod (void) const
  1692. {
  1693.   Matrix retval;
  1694.  
  1695.   int nr = rows ();
  1696.   int nc = cols ();
  1697.  
  1698.   if (nr == 1)
  1699.     {
  1700.       retval.resize (1, nc);
  1701.       if (nc > 0)
  1702.     {
  1703.       double prod = elem (0, 0);
  1704.       for (int j = 0; j < nc; j++)
  1705.         {
  1706.           retval.elem (0, j) = prod;
  1707.           if (j < nc - 1)
  1708.         prod *= elem (0, j+1);
  1709.         }
  1710.     }
  1711.     }
  1712.   else if (nc == 1)
  1713.     {
  1714.       retval.resize (nr, 1);
  1715.       if (nr > 0)
  1716.     {
  1717.       double prod = elem (0, 0);
  1718.       for (int i = 0; i < nr; i++)
  1719.         {
  1720.           retval.elem (i, 0) = prod;
  1721.           if (i < nr - 1)
  1722.         prod *= elem (i+1, 0);
  1723.         }
  1724.     }
  1725.     }
  1726.   else
  1727.     {
  1728.       retval.resize (nr, nc);
  1729.       if (nr > 0 && nc > 0)
  1730.     {
  1731.       for (int j = 0; j < nc; j++)
  1732.         {
  1733.           double prod = elem (0, j);
  1734.           for (int i = 0; i < nr; i++)
  1735.         {
  1736.           retval.elem (i, j) = prod;
  1737.           if (i < nr - 1)
  1738.             prod *= elem (i+1, j);
  1739.         }
  1740.         }
  1741.     }
  1742.     }
  1743.   return retval;
  1744. }
  1745.  
  1746. Matrix
  1747. Matrix::cumsum (void) const
  1748. {
  1749.   Matrix retval;
  1750.  
  1751.   int nr = rows ();
  1752.   int nc = cols ();
  1753.  
  1754.   if (nr == 1)
  1755.     {
  1756.       retval.resize (1, nc);
  1757.       if (nc > 0)
  1758.     {
  1759.       double sum = elem (0, 0);
  1760.       for (int j = 0; j < nc; j++)
  1761.         {
  1762.           retval.elem (0, j) = sum;
  1763.           if (j < nc - 1)
  1764.         sum += elem (0, j+1);
  1765.         }
  1766.     }
  1767.     }
  1768.   else if (nc == 1)
  1769.     {
  1770.       retval.resize (nr, 1);
  1771.       if (nr > 0)
  1772.     {
  1773.       double sum = elem (0, 0);
  1774.       for (int i = 0; i < nr; i++)
  1775.         {
  1776.           retval.elem (i, 0) = sum;
  1777.           if (i < nr - 1)
  1778.         sum += elem (i+1, 0);
  1779.         }
  1780.     }
  1781.     }
  1782.   else
  1783.     {
  1784.       retval.resize (nr, nc);
  1785.       if (nr > 0 && nc > 0)
  1786.     {
  1787.       for (int j = 0; j < nc; j++)
  1788.         {
  1789.           double sum = elem (0, j);
  1790.           for (int i = 0; i < nr; i++)
  1791.         {
  1792.           retval.elem (i, j) = sum;
  1793.           if (i < nr - 1)
  1794.             sum += elem (i+1, j);
  1795.         }
  1796.         }
  1797.     }
  1798.     }
  1799.   return retval;
  1800. }
  1801.  
  1802. Matrix
  1803. Matrix::prod (void) const
  1804. {
  1805.   Matrix retval;
  1806.  
  1807.   int nr = rows ();
  1808.   int nc = cols ();
  1809.  
  1810.   if (nr == 1)
  1811.     {
  1812.       retval.resize (1, 1);
  1813.       retval.elem (0, 0) = 1.0;
  1814.       for (int j = 0; j < nc; j++)
  1815.     retval.elem (0, 0) *= elem (0, j);
  1816.     }
  1817.   else if (nc == 1)
  1818.     {
  1819.       retval.resize (1, 1);
  1820.       retval.elem (0, 0) = 1.0;
  1821.       for (int i = 0; i < nr; i++)
  1822.     retval.elem (0, 0) *= elem (i, 0);
  1823.     }
  1824.   else
  1825.     {
  1826.       if (nc == 0)
  1827.     {
  1828.       retval.resize (1, 1);
  1829.       retval.elem (0, 0) = 1.0;
  1830.     }
  1831.       else
  1832.     retval.resize (1, nc);
  1833.  
  1834.       for (int j = 0; j < nc; j++)
  1835.     {
  1836.       retval.elem (0, j) = 1.0;
  1837.       for (int i = 0; i < nr; i++)
  1838.         retval.elem (0, j) *= elem (i, j);
  1839.     }
  1840.     }
  1841.   return retval;
  1842. }
  1843.  
  1844. Matrix
  1845. Matrix::sum (void) const
  1846. {
  1847.   Matrix retval;
  1848.  
  1849.   int nr = rows ();
  1850.   int nc = cols ();
  1851.  
  1852.   if (nr == 1)
  1853.     {
  1854.       retval.resize (1, 1);
  1855.       retval.elem (0, 0) = 0.0;
  1856.       for (int j = 0; j < nc; j++)
  1857.     retval.elem (0, 0) += elem (0, j);
  1858.     }
  1859.   else if (nc == 1)
  1860.     {
  1861.       retval.resize (1, 1);
  1862.       retval.elem (0, 0) = 0.0;
  1863.       for (int i = 0; i < nr; i++)
  1864.     retval.elem (0, 0) += elem (i, 0);
  1865.     }
  1866.   else
  1867.     {
  1868.       if (nc == 0)
  1869.     {
  1870.       retval.resize (1, 1);
  1871.       retval.elem (0, 0) = 0.0;
  1872.     }
  1873.       else
  1874.     retval.resize (1, nc);
  1875.  
  1876.       for (int j = 0; j < nc; j++)
  1877.     {
  1878.       retval.elem (0, j) = 0.0;
  1879.       for (int i = 0; i < nr; i++)
  1880.         retval.elem (0, j) += elem (i, j);
  1881.     }
  1882.     }
  1883.   return retval;
  1884. }
  1885.  
  1886. Matrix
  1887. Matrix::sumsq (void) const
  1888. {
  1889.   Matrix retval;
  1890.  
  1891.   int nr = rows ();
  1892.   int nc = cols ();
  1893.  
  1894.   if (nr == 1)
  1895.     {
  1896.       retval.resize (1, 1);
  1897.       retval.elem (0, 0) = 0.0;
  1898.       for (int j = 0; j < nc; j++)
  1899.     {
  1900.       double d = elem (0, j);
  1901.       retval.elem (0, 0) += d * d;
  1902.     }
  1903.     }
  1904.   else if (nc == 1)
  1905.     {
  1906.       retval.resize (1, 1);
  1907.       retval.elem (0, 0) = 0.0;
  1908.       for (int i = 0; i < nr; i++)
  1909.     {
  1910.       double d = elem (i, 0);
  1911.       retval.elem (0, 0) += d * d;
  1912.     }
  1913.     }
  1914.   else
  1915.     {
  1916.       retval.resize (1, nc);
  1917.       for (int j = 0; j < nc; j++)
  1918.     {
  1919.       retval.elem (0, j) = 0.0;
  1920.       for (int i = 0; i < nr; i++)
  1921.         {
  1922.           double d = elem (i, j);
  1923.           retval.elem (0, j) += d * d;
  1924.         }
  1925.     }
  1926.     }
  1927.   return retval;
  1928. }
  1929.  
  1930. ColumnVector
  1931. Matrix::diag (void) const
  1932. {
  1933.   return diag (0);
  1934. }
  1935.  
  1936. ColumnVector
  1937. Matrix::diag (int k) const
  1938. {
  1939.   int nnr = rows ();
  1940.   int nnc = cols ();
  1941.   if (k > 0)
  1942.     nnc -= k;
  1943.   else if (k < 0)
  1944.     nnr += k;
  1945.  
  1946.   ColumnVector d;
  1947.  
  1948.   if (nnr > 0 && nnc > 0)
  1949.     {
  1950.       int ndiag = (nnr < nnc) ? nnr : nnc;
  1951.  
  1952.       d.resize (ndiag);
  1953.  
  1954.       if (k > 0)
  1955.     {
  1956.       for (int i = 0; i < ndiag; i++)
  1957.         d.elem (i) = elem (i, i+k);
  1958.     }
  1959.       else if ( k < 0)
  1960.     {
  1961.       for (int i = 0; i < ndiag; i++)
  1962.         d.elem (i) = elem (i-k, i);
  1963.     }
  1964.       else
  1965.     {
  1966.       for (int i = 0; i < ndiag; i++)
  1967.         d.elem (i) = elem (i, i);
  1968.     }
  1969.     }
  1970.   else
  1971.     cerr << "diag: requested diagonal out of range\n";
  1972.  
  1973.   return d;
  1974. }
  1975.  
  1976. ColumnVector
  1977. Matrix::row_min (void) const
  1978. {
  1979.   ColumnVector result;
  1980.  
  1981.   int nr = rows ();
  1982.   int nc = cols ();
  1983.  
  1984.   if (nr > 0 && nc > 0)
  1985.     {
  1986.       result.resize (nr);
  1987.  
  1988.       for (int i = 0; i < nr; i++)
  1989.     {
  1990.       double res = elem (i, 0);
  1991.       for (int j = 1; j < nc; j++)
  1992.         if (elem (i, j) < res)
  1993.           res = elem (i, j);
  1994.       result.elem (i) = res;
  1995.     }
  1996.     }
  1997.  
  1998.   return result;
  1999. }
  2000.  
  2001. ColumnVector
  2002. Matrix::row_min_loc (void) const
  2003. {
  2004.   ColumnVector result;
  2005.  
  2006.   int nr = rows ();
  2007.   int nc = cols ();
  2008.  
  2009.   if (nr > 0 && nc > 0)
  2010.     {
  2011.       result.resize (nr);
  2012.  
  2013.       for (int i = 0; i < nr; i++)
  2014.         {
  2015.           int res = 0;
  2016.           for (int j = 0; j < nc; j++)
  2017.             if (elem (i, j) < elem (i, res))
  2018.               res = j;
  2019.           result.elem (i) = (double) (res + 1);
  2020.         }
  2021.     }
  2022.  
  2023.   return result;
  2024. }
  2025.  
  2026. ColumnVector
  2027. Matrix::row_max (void) const
  2028. {
  2029.   ColumnVector result;
  2030.  
  2031.   int nr = rows ();
  2032.   int nc = cols ();
  2033.  
  2034.   if (nr > 0 && nc > 0)
  2035.     {
  2036.       result.resize (nr);
  2037.  
  2038.       for (int i = 0; i < nr; i++)
  2039.     {
  2040.       double res = elem (i, 0);
  2041.       for (int j = 1; j < nc; j++)
  2042.         if (elem (i, j) > res)
  2043.           res = elem (i, j);
  2044.       result.elem (i) = res;
  2045.     }
  2046.     }
  2047.  
  2048.   return result;
  2049. }
  2050.  
  2051. ColumnVector
  2052. Matrix::row_max_loc (void) const
  2053. {
  2054.   ColumnVector result;
  2055.  
  2056.   int nr = rows ();
  2057.   int nc = cols ();
  2058.  
  2059.   if (nr > 0 && nc > 0)
  2060.     {
  2061.       result.resize (nr);
  2062.  
  2063.       for (int i = 0; i < nr; i++)
  2064.         {
  2065.           int res = 0;
  2066.           for (int j = 0; j < nc; j++)
  2067.             if (elem (i, j) > elem (i, res))
  2068.               res = j;
  2069.           result.elem (i) = (double) (res + 1);
  2070.         }
  2071.     }
  2072.  
  2073.   return result;
  2074. }
  2075.  
  2076. RowVector
  2077. Matrix::column_min (void) const
  2078. {
  2079.   RowVector result;
  2080.  
  2081.   int nr = rows ();
  2082.   int nc = cols ();
  2083.  
  2084.   if (nr > 0 && nc > 0)
  2085.     {
  2086.       result.resize (nc);
  2087.  
  2088.       for (int j = 0; j < nc; j++)
  2089.     {
  2090.       double res = elem (0, j);
  2091.       for (int i = 1; i < nr; i++)
  2092.         if (elem (i, j) < res)
  2093.           res = elem (i, j);
  2094.       result.elem (j) = res;
  2095.     }
  2096.     }
  2097.  
  2098.   return result;
  2099. }
  2100. RowVector
  2101. Matrix::column_min_loc (void) const
  2102. {
  2103.   RowVector result;
  2104.  
  2105.   int nr = rows ();
  2106.   int nc = cols ();
  2107.  
  2108.   if (nr > 0 && nc > 0)
  2109.     {
  2110.       result.resize (nc);
  2111.  
  2112.       for (int j = 0; j < nc; j++)
  2113.         {
  2114.           int res = 0;
  2115.           for (int i = 0; i < nr; i++)
  2116.             if (elem (i, j) < elem (res, j))
  2117.               res = i;
  2118.           result.elem (j) = (double) (res + 1);
  2119.         }
  2120.     }
  2121.  
  2122.   return result;
  2123. }
  2124.  
  2125.  
  2126. RowVector
  2127. Matrix::column_max (void) const
  2128. {
  2129.   RowVector result;
  2130.  
  2131.   int nr = rows ();
  2132.   int nc = cols ();
  2133.  
  2134.   if (nr > 0 && nc > 0)
  2135.     {
  2136.       result.resize (nc);
  2137.  
  2138.       for (int j = 0; j < nc; j++)
  2139.     {
  2140.       double res = elem (0, j);
  2141.       for (int i = 1; i < nr; i++)
  2142.         if (elem (i, j) > res)
  2143.           res = elem (i, j);
  2144.       result.elem (j) = res;
  2145.     }
  2146.     }
  2147.  
  2148.   return result;
  2149. }
  2150.  
  2151. RowVector
  2152. Matrix::column_max_loc (void) const
  2153. {
  2154.   RowVector result;
  2155.  
  2156.   int nr = rows ();
  2157.   int nc = cols ();
  2158.  
  2159.   if (nr > 0 && nc > 0)
  2160.     {
  2161.       result.resize (nc);
  2162.  
  2163.       for (int j = 0; j < nc; j++)
  2164.         {
  2165.           int res = 0;
  2166.           for (int i = 0; i < nr; i++)
  2167.             if (elem (i, j) > elem (res, j))
  2168.               res = i;
  2169.           result.elem (j) = (double) (res + 1);
  2170.         }
  2171.     }
  2172.  
  2173.   return result;
  2174. }
  2175.  
  2176. ostream&
  2177. operator << (ostream& os, const Matrix& a)
  2178. {
  2179. //  int field_width = os.precision () + 7;
  2180.   for (int i = 0; i < a.rows (); i++)
  2181.     {
  2182.       for (int j = 0; j < a.cols (); j++)
  2183.     os << " " /* setw (field_width) */ << a.elem (i, j);
  2184.       os << "\n";
  2185.     }
  2186.   return os;
  2187. }
  2188.  
  2189. istream&
  2190. operator >> (istream& is, Matrix& a)
  2191. {
  2192.   int nr = a.rows ();
  2193.   int nc = a.cols ();
  2194.  
  2195.   if (nr < 1 || nc < 1)
  2196.     is.clear (ios::badbit);
  2197.   else
  2198.     {
  2199.       double tmp;
  2200.       for (int i = 0; i < nr; i++)
  2201.     for (int j = 0; j < nc; j++)
  2202.       {
  2203.         is >> tmp;
  2204.         if (is)
  2205.           a.elem (i, j) = tmp;
  2206.         else
  2207.           break;
  2208.       }
  2209.     }
  2210.  
  2211.   return is;
  2212. }
  2213.  
  2214. /*
  2215.  * Read an array of data froma file in binary format.
  2216.  */
  2217. int
  2218. Matrix::read (FILE *fptr, char *type)
  2219. {
  2220. // Allocate buffer pointers.
  2221.  
  2222.   union
  2223.     {
  2224.       void *vd;
  2225.       char *ch;
  2226.       u_char *uc;
  2227.       short *sh;
  2228.       u_short *us;
  2229.       int *in;
  2230.       u_int *ui;
  2231.       long *ln;
  2232.       u_long *ul;
  2233.       float *fl;
  2234.       double *db;
  2235.     }
  2236.   buf;
  2237.  
  2238. // Convert data to double.
  2239.  
  2240.   if (! type)
  2241.     {
  2242.       (*current_liboctave_error_handler)
  2243.     ("fread: invalid NULL type parameter");
  2244.       return 0;
  2245.     }    
  2246.  
  2247.   int count;
  2248.   int nitems = length ();
  2249.  
  2250.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  2251.  
  2252. #define DO_FREAD(TYPE,ELEM) \
  2253.   do \
  2254.     { \
  2255.       size_t size = sizeof (TYPE); \
  2256.       buf.ch = new char [size * nitems]; \
  2257.       count = fread (buf.ch, size, nitems, fptr); \
  2258.       for (int k = 0; k < count; k++) \
  2259.     d[k] = buf.ELEM[k]; \
  2260.       delete [] buf.ch; \
  2261.     } \
  2262.   while (0)
  2263.  
  2264.   if (strcasecmp (type, "double") == 0)
  2265.     DO_FREAD (double, db);
  2266.   else if (strcasecmp (type, "char") == 0)
  2267.     DO_FREAD (char, ch);
  2268.   else if (strcasecmp (type, "uchar") == 0)
  2269.     DO_FREAD (u_char, uc);
  2270.   else if (strcasecmp (type, "short") == 0)
  2271.     DO_FREAD (short, sh);
  2272.   else if (strcasecmp (type, "ushort") == 0)
  2273.     DO_FREAD (u_short, us);
  2274.   else if (strcasecmp (type, "int") == 0)
  2275.     DO_FREAD (int, in);
  2276.   else if (strcasecmp (type, "uint") == 0)
  2277.     DO_FREAD (u_int, ui);
  2278.   else if (strcasecmp (type, "long") == 0)
  2279.     DO_FREAD (long, ul);
  2280.   else if (strcasecmp (type, "float") == 0)
  2281.     DO_FREAD (float, fl);
  2282.   else
  2283.     {
  2284.       (*current_liboctave_error_handler)
  2285.     ("fread: invalid NULL type parameter");
  2286.       return 0;
  2287.     }
  2288.  
  2289.   return count;
  2290. }
  2291.  
  2292. /*
  2293.  * Write the data array to a file in binary format.
  2294.  */
  2295. int
  2296. Matrix::write (FILE *fptr, char *type)
  2297. {
  2298. // Allocate buffer pointers.
  2299.  
  2300.   union
  2301.     {
  2302.       void *vd;
  2303.       char *ch;
  2304.       u_char *uc;
  2305.       short *sh;
  2306.       u_short *us;
  2307.       int *in;
  2308.       u_int *ui;
  2309.       long *ln;
  2310.       u_long *ul;
  2311.       float *fl;
  2312.       double *db;
  2313.     }
  2314.   buf;
  2315.  
  2316.   int nitems = length ();
  2317.  
  2318.   double *d = fortran_vec ();
  2319.  
  2320. // Convert from double to correct size.
  2321.  
  2322.   if (! type)
  2323.     {
  2324.       (*current_liboctave_error_handler)
  2325.     ("fwrite: invalid NULL type parameter");
  2326.       return 0;
  2327.     }    
  2328.  
  2329.   size_t size;
  2330.   int count;
  2331.  
  2332. #define DO_FWRITE(TYPE,ELEM) \
  2333.   do \
  2334.     { \
  2335.       size = sizeof (TYPE); \
  2336.       buf.ELEM = new TYPE [nitems]; \
  2337.       for (int k = 0; k < nitems; k++) \
  2338.     buf.ELEM[k] = (TYPE) d[k]; \
  2339.       count = fwrite (buf.ELEM, size, nitems, fptr); \
  2340.       delete [] buf.ELEM; \
  2341.     } \
  2342.   while (0)
  2343.  
  2344.   if (strcasecmp (type, "double") == 0)
  2345.     DO_FWRITE (double, db);
  2346.   else if (strcasecmp (type, "char") == 0)
  2347.     DO_FWRITE (char, ch);
  2348.   else if (strcasecmp (type, "uchar") == 0)
  2349.     DO_FWRITE (u_char, uc);
  2350.   else if (strcasecmp (type, "short") == 0)
  2351.     DO_FWRITE (short, sh);
  2352.   else if (strcasecmp (type, "ushort") == 0)
  2353.     DO_FWRITE (u_short, us);
  2354.   else if (strcasecmp (type, "int") == 0)
  2355.     DO_FWRITE (int, in);
  2356.   else if (strcasecmp (type, "uint") == 0)
  2357.     DO_FWRITE (u_int, ui);
  2358.   else if (strcasecmp (type, "long") == 0)
  2359.     DO_FWRITE (long, ln);
  2360.   else if (strcasecmp (type, "ulong") == 0)
  2361.     DO_FWRITE (u_long, ul);
  2362.   else if (strcasecmp (type, "float") == 0)
  2363.     DO_FWRITE (float, fl);
  2364.   else
  2365.     {
  2366.       (*current_liboctave_error_handler)
  2367.     ("fwrite: unrecognized type parameter %s", type);
  2368.       return 0;
  2369.     }
  2370.  
  2371.   return count;
  2372. }
  2373.  
  2374. /*
  2375. ;;; Local Variables: ***
  2376. ;;; mode: C++ ***
  2377. ;;; page-delimiter: "^/\\*" ***
  2378. ;;; End: ***
  2379. */
  2380.