home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fs.zip / octave / octave-2.1.23 / liboctave / CMatrix.cc < prev    next >
C/C++ Source or Header  |  2000-01-15  |  59KB  |  2,970 lines
  1. // Matrix manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996, 1997 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <cfloat>
  33.  
  34. #include <iostream.h>
  35.  
  36. // XXX FIXME XXX
  37. #ifdef HAVE_SYS_TYPES_H
  38. #include <sys/types.h>
  39. #endif
  40.  
  41. #include "CMatrix.h"
  42. #include "CmplxAEPBAL.h"
  43. #include "CmplxDET.h"
  44. #include "CmplxSCHUR.h"
  45. #include "CmplxSVD.h"
  46. #include "f77-fcn.h"
  47. #include "lo-error.h"
  48. #include "lo-ieee.h"
  49. #include "lo-mappers.h"
  50. #include "lo-utils.h"
  51. #include "mx-base.h"
  52. #include "mx-cm-dm.h"
  53. #include "mx-dm-cm.h"
  54. #include "mx-cm-s.h"
  55. #include "mx-inlines.cc"
  56. #include "oct-cmplx.h"
  57.  
  58. // Fortran functions we call.
  59.  
  60. extern "C"
  61. {
  62.   int F77_FCN (zgebal, ZGEBAL) (const char*, const int&, Complex*,
  63.                                 const int&, int&, int&, double*, int&,
  64.                                 long, long);
  65.  
  66.   int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
  67.                                 const int&, const int&, double*,
  68.                                 const int&, double*, const int&,
  69.                                 int&, long, long);
  70.  
  71.   int F77_FCN (zgemm, ZGEMM) (const char*, const char*, const int&,
  72.                   const int&, const int&, const Complex&,
  73.                   const Complex*, const int&,
  74.                   const Complex*, const int&,
  75.                   const Complex&, Complex*, const int&, 
  76.                   long, long);
  77.  
  78.   int F77_FCN (zgeco, ZGECO) (Complex*, const int&, const int&, int*,
  79.                   double&, Complex*);
  80.  
  81.   int F77_FCN (zgedi, ZGEDI) (Complex*, const int&, const int&, int*,
  82.                   Complex*, Complex*, const int&);
  83.  
  84.   int F77_FCN (zgesl, ZGESL) (Complex*, const int&, const int&, int*,
  85.                   Complex*, const int&);
  86.  
  87.   int F77_FCN (zgelss, ZGELSS) (const int&, const int&, const int&,
  88.                 Complex*, const int&, Complex*,
  89.                 const int&, double*, double&, int&,
  90.                 Complex*, const int&, double*, int&);
  91.  
  92.   // Note that the original complex fft routines were not written for
  93.   // double complex arguments.  They have been modified by adding an
  94.   // implicit double precision (a-h,o-z) statement at the beginning of
  95.   // each subroutine.
  96.  
  97.   int F77_FCN (cffti, CFFTI) (const int&, Complex*);
  98.  
  99.   int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);
  100.  
  101.   int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);
  102.  
  103.   int F77_FCN (zlartg, ZLARTG) (const Complex&, const Complex&,
  104.                 double&, Complex&, Complex&);
  105.  
  106.   int F77_FCN (ztrsyl, ZTRSYL) (const char*, const char*, const int&,
  107.                 const int&, const int&,
  108.                 const Complex*, const int&,
  109.                 const Complex*, const int&, 
  110.                 const Complex*, const int&, double&,
  111.                 int&, long, long);
  112.  
  113.   int F77_FCN (xzlange, XZLANGE) (const char*, const int&,
  114.                   const int&, const Complex*,
  115.                   const int&, double*, double&); 
  116. }
  117.  
  118. static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
  119.  
  120. // Complex Matrix class
  121.  
  122. ComplexMatrix::ComplexMatrix (const Matrix& a)
  123.   : MArray2<Complex> (a.rows (), a.cols ())
  124. {
  125.   for (int j = 0; j < cols (); j++)
  126.     for (int i = 0; i < rows (); i++)
  127.       elem (i, j) = a.elem (i, j);
  128. }
  129.  
  130. ComplexMatrix::ComplexMatrix (const RowVector& rv)
  131.   : MArray2<Complex> (1, rv.length (), 0.0)
  132. {
  133.   for (int i = 0; i < rv.length (); i++)
  134.     elem (0, i) = rv.elem (i);
  135. }
  136.  
  137. ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
  138.   : MArray2<Complex> (cv.length (), 1, 0.0)
  139. {
  140.   for (int i = 0; i < cv.length (); i++)
  141.     elem (i, 0) = cv.elem (i);
  142. }
  143.  
  144. ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
  145.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  146. {
  147.   for (int i = 0; i < a.length (); i++)
  148.     elem (i, i) = a.elem (i, i);
  149. }
  150.  
  151. ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
  152.   : MArray2<Complex> (1, rv.length (), 0.0)
  153. {
  154.   for (int i = 0; i < rv.length (); i++)
  155.     elem (0, i) = rv.elem (i);
  156. }
  157.  
  158. ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
  159.   : MArray2<Complex> (cv.length (), 1, 0.0)
  160. {
  161.   for (int i = 0; i < cv.length (); i++)
  162.     elem (i, 0) = cv.elem (i);
  163. }
  164.  
  165. ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
  166.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  167. {
  168.   for (int i = 0; i < a.length (); i++)
  169.     elem (i, i) = a.elem (i, i);
  170. }
  171.  
  172. // XXX FIXME XXX -- could we use a templated mixed-type copy function
  173. // here?
  174.  
  175. ComplexMatrix::ComplexMatrix (const boolMatrix& a)
  176.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  177. {
  178.   for (int i = 0; i < a.cols (); i++)
  179.     for (int j = 0; j < a.rows (); j++)
  180.       elem (i, j) = a.elem (i, j);
  181. }
  182.  
  183. ComplexMatrix::ComplexMatrix (const charMatrix& a)
  184.   : MArray2<Complex> (a.rows (), a.cols (), 0.0)
  185. {
  186.   for (int i = 0; i < a.cols (); i++)
  187.     for (int j = 0; j < a.rows (); j++)
  188.       elem (i, j) = a.elem (i, j);
  189. }
  190.  
  191. bool
  192. ComplexMatrix::operator == (const ComplexMatrix& a) const
  193. {
  194.   if (rows () != a.rows () || cols () != a.cols ())
  195.     return false;
  196.  
  197.   return equal (data (), a.data (), length ());
  198. }
  199.  
  200. bool
  201. ComplexMatrix::operator != (const ComplexMatrix& a) const
  202. {
  203.   return !(*this == a);
  204. }
  205.  
  206. bool
  207. ComplexMatrix::is_hermitian (void) const
  208. {
  209.   int nr = rows ();
  210.   int nc = cols ();
  211.  
  212.   if (is_square () && nr > 0)
  213.     {
  214.       for (int i = 0; i < nr; i++)
  215.     for (int j = i; j < nc; j++)
  216.       if (elem (i, j) != conj (elem (j, i)))
  217.         return false;
  218.  
  219.       return true;
  220.     }
  221.  
  222.   return false;
  223. }
  224.  
  225. // destructive insert/delete/reorder operations
  226.  
  227. ComplexMatrix&
  228. ComplexMatrix::insert (const Matrix& a, int r, int c)
  229. {
  230.   int a_nr = a.rows ();
  231.   int a_nc = a.cols ();
  232.  
  233.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  234.     {
  235.       (*current_liboctave_error_handler) ("range error for insert");
  236.       return *this;
  237.     }
  238.  
  239.   for (int j = 0; j < a_nc; j++)
  240.     for (int i = 0; i < a_nr; i++)
  241.       elem (r+i, c+j) = a.elem (i, j);
  242.  
  243.   return *this;
  244. }
  245.  
  246. ComplexMatrix&
  247. ComplexMatrix::insert (const RowVector& a, int r, int c)
  248. {
  249.   int a_len = a.length ();
  250.   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
  251.     {
  252.       (*current_liboctave_error_handler) ("range error for insert");
  253.       return *this;
  254.     }
  255.  
  256.   for (int i = 0; i < a_len; i++)
  257.     elem (r, c+i) = a.elem (i);
  258.  
  259.   return *this;
  260. }
  261.  
  262. ComplexMatrix&
  263. ComplexMatrix::insert (const ColumnVector& a, int r, int c)
  264. {
  265.   int a_len = a.length ();
  266.   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
  267.     {
  268.       (*current_liboctave_error_handler) ("range error for insert");
  269.       return *this;
  270.     }
  271.  
  272.   for (int i = 0; i < a_len; i++)
  273.     elem (r+i, c) = a.elem (i);
  274.  
  275.   return *this;
  276. }
  277.  
  278. ComplexMatrix&
  279. ComplexMatrix::insert (const DiagMatrix& a, int r, int c)
  280. {
  281.   int a_nr = a.rows ();
  282.   int a_nc = a.cols ();
  283.  
  284.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  285.     {
  286.       (*current_liboctave_error_handler) ("range error for insert");
  287.       return *this;
  288.     }
  289.  
  290.   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
  291.  
  292.   for (int i = 0; i < a.length (); i++)
  293.     elem (r+i, c+i) = a.elem (i, i);
  294.  
  295.   return *this;
  296. }
  297.  
  298. ComplexMatrix&
  299. ComplexMatrix::insert (const ComplexMatrix& a, int r, int c)
  300. {
  301.   Array2<Complex>::insert (a, r, c);
  302.   return *this;
  303. }
  304.  
  305. ComplexMatrix&
  306. ComplexMatrix::insert (const ComplexRowVector& a, int r, int c)
  307. {
  308.   int a_len = a.length ();
  309.   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
  310.     {
  311.       (*current_liboctave_error_handler) ("range error for insert");
  312.       return *this;
  313.     }
  314.  
  315.   for (int i = 0; i < a_len; i++)
  316.     elem (r, c+i) = a.elem (i);
  317.  
  318.   return *this;
  319. }
  320.  
  321. ComplexMatrix&
  322. ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c)
  323. {
  324.   int a_len = a.length ();
  325.   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
  326.     {
  327.       (*current_liboctave_error_handler) ("range error for insert");
  328.       return *this;
  329.     }
  330.  
  331.   for (int i = 0; i < a_len; i++)
  332.     elem (r+i, c) = a.elem (i);
  333.  
  334.   return *this;
  335. }
  336.  
  337. ComplexMatrix&
  338. ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c)
  339. {
  340.   int a_nr = a.rows ();
  341.   int a_nc = a.cols ();
  342.  
  343.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  344.     {
  345.       (*current_liboctave_error_handler) ("range error for insert");
  346.       return *this;
  347.     }
  348.  
  349.   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
  350.  
  351.   for (int i = 0; i < a.length (); i++)
  352.     elem (r+i, c+i) = a.elem (i, i);
  353.  
  354.   return *this;
  355. }
  356.  
  357. ComplexMatrix&
  358. ComplexMatrix::fill (double val)
  359. {
  360.   int nr = rows ();
  361.   int nc = cols ();
  362.   if (nr > 0 && nc > 0)
  363.     for (int j = 0; j < nc; j++)
  364.       for (int i = 0; i < nr; i++)
  365.     elem (i, j) = val;
  366.  
  367.   return *this;
  368. }
  369.  
  370. ComplexMatrix&
  371. ComplexMatrix::fill (const Complex& val)
  372. {
  373.   int nr = rows ();
  374.   int nc = cols ();
  375.   if (nr > 0 && nc > 0)
  376.     for (int j = 0; j < nc; j++)
  377.       for (int i = 0; i < nr; i++)
  378.     elem (i, j) = val;
  379.  
  380.   return *this;
  381. }
  382.  
  383. ComplexMatrix&
  384. ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2)
  385. {
  386.   int nr = rows ();
  387.   int nc = cols ();
  388.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  389.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  390.     {
  391.       (*current_liboctave_error_handler) ("range error for fill");
  392.       return *this;
  393.     }
  394.  
  395.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  396.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  397.  
  398.   for (int j = c1; j <= c2; j++)
  399.     for (int i = r1; i <= r2; i++)
  400.       elem (i, j) = val;
  401.  
  402.   return *this;
  403. }
  404.  
  405. ComplexMatrix&
  406. ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2)
  407. {
  408.   int nr = rows ();
  409.   int nc = cols ();
  410.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  411.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  412.     {
  413.       (*current_liboctave_error_handler) ("range error for fill");
  414.       return *this;
  415.     }
  416.  
  417.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  418.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  419.  
  420.   for (int j = c1; j <= c2; j++)
  421.     for (int i = r1; i <= r2; i++)
  422.       elem (i, j) = val;
  423.  
  424.   return *this;
  425. }
  426.  
  427. ComplexMatrix
  428. ComplexMatrix::append (const Matrix& a) const
  429. {
  430.   int nr = rows ();
  431.   int nc = cols ();
  432.   if (nr != a.rows ())
  433.     {
  434.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  435.       return *this;
  436.     }
  437.  
  438.   int nc_insert = nc;
  439.   ComplexMatrix retval (nr, nc + a.cols ());
  440.   retval.insert (*this, 0, 0);
  441.   retval.insert (a, 0, nc_insert);
  442.   return retval;
  443. }
  444.  
  445. ComplexMatrix
  446. ComplexMatrix::append (const RowVector& a) const
  447. {
  448.   int nr = rows ();
  449.   int nc = cols ();
  450.   if (nr != 1)
  451.     {
  452.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  453.       return *this;
  454.     }
  455.  
  456.   int nc_insert = nc;
  457.   ComplexMatrix retval (nr, nc + a.length ());
  458.   retval.insert (*this, 0, 0);
  459.   retval.insert (a, 0, nc_insert);
  460.   return retval;
  461. }
  462.  
  463. ComplexMatrix
  464. ComplexMatrix::append (const ColumnVector& a) const
  465. {
  466.   int nr = rows ();
  467.   int nc = cols ();
  468.   if (nr != a.length ())
  469.     {
  470.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  471.       return *this;
  472.     }
  473.  
  474.   int nc_insert = nc;
  475.   ComplexMatrix retval (nr, nc + 1);
  476.   retval.insert (*this, 0, 0);
  477.   retval.insert (a, 0, nc_insert);
  478.   return retval;
  479. }
  480.  
  481. ComplexMatrix
  482. ComplexMatrix::append (const DiagMatrix& a) const
  483. {
  484.   int nr = rows ();
  485.   int nc = cols ();
  486.   if (nr != a.rows ())
  487.     {
  488.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  489.       return *this;
  490.     }
  491.  
  492.   int nc_insert = nc;
  493.   ComplexMatrix retval (nr, nc + a.cols ());
  494.   retval.insert (*this, 0, 0);
  495.   retval.insert (a, 0, nc_insert);
  496.   return retval;
  497. }
  498.  
  499. ComplexMatrix
  500. ComplexMatrix::append (const ComplexMatrix& a) const
  501. {
  502.   int nr = rows ();
  503.   int nc = cols ();
  504.   if (nr != a.rows ())
  505.     {
  506.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  507.       return *this;
  508.     }
  509.  
  510.   int nc_insert = nc;
  511.   ComplexMatrix retval (nr, nc + a.cols ());
  512.   retval.insert (*this, 0, 0);
  513.   retval.insert (a, 0, nc_insert);
  514.   return retval;
  515. }
  516.  
  517. ComplexMatrix
  518. ComplexMatrix::append (const ComplexRowVector& a) const
  519. {
  520.   int nr = rows ();
  521.   int nc = cols ();
  522.   if (nr != 1)
  523.     {
  524.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  525.       return *this;
  526.     }
  527.  
  528.   int nc_insert = nc;
  529.   ComplexMatrix retval (nr, nc + a.length ());
  530.   retval.insert (*this, 0, 0);
  531.   retval.insert (a, 0, nc_insert);
  532.   return retval;
  533. }
  534.  
  535. ComplexMatrix
  536. ComplexMatrix::append (const ComplexColumnVector& a) const
  537. {
  538.   int nr = rows ();
  539.   int nc = cols ();
  540.   if (nr != a.length ())
  541.     {
  542.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  543.       return *this;
  544.     }
  545.  
  546.   int nc_insert = nc;
  547.   ComplexMatrix retval (nr, nc + 1);
  548.   retval.insert (*this, 0, 0);
  549.   retval.insert (a, 0, nc_insert);
  550.   return retval;
  551. }
  552.  
  553. ComplexMatrix
  554. ComplexMatrix::append (const ComplexDiagMatrix& a) const
  555. {
  556.   int nr = rows ();
  557.   int nc = cols ();
  558.   if (nr != a.rows ())
  559.     {
  560.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  561.       return *this;
  562.     }
  563.  
  564.   int nc_insert = nc;
  565.   ComplexMatrix retval (nr, nc + a.cols ());
  566.   retval.insert (*this, 0, 0);
  567.   retval.insert (a, 0, nc_insert);
  568.   return retval;
  569. }
  570.  
  571. ComplexMatrix
  572. ComplexMatrix::stack (const Matrix& a) const
  573. {
  574.   int nr = rows ();
  575.   int nc = cols ();
  576.   if (nc != a.cols ())
  577.     {
  578.       (*current_liboctave_error_handler)
  579.     ("column dimension mismatch for stack");
  580.       return *this;
  581.     }
  582.  
  583.   int nr_insert = nr;
  584.   ComplexMatrix retval (nr + a.rows (), nc);
  585.   retval.insert (*this, 0, 0);
  586.   retval.insert (a, nr_insert, 0);
  587.   return retval;
  588. }
  589.  
  590. ComplexMatrix
  591. ComplexMatrix::stack (const RowVector& a) const
  592. {
  593.   int nr = rows ();
  594.   int nc = cols ();
  595.   if (nc != a.length ())
  596.     {
  597.       (*current_liboctave_error_handler)
  598.     ("column dimension mismatch for stack");
  599.       return *this;
  600.     }
  601.  
  602.   int nr_insert = nr;
  603.   ComplexMatrix retval (nr + 1, nc);
  604.   retval.insert (*this, 0, 0);
  605.   retval.insert (a, nr_insert, 0);
  606.   return retval;
  607. }
  608.  
  609. ComplexMatrix
  610. ComplexMatrix::stack (const ColumnVector& a) const
  611. {
  612.   int nr = rows ();
  613.   int nc = cols ();
  614.   if (nc != 1)
  615.     {
  616.       (*current_liboctave_error_handler)
  617.     ("column dimension mismatch for stack");
  618.       return *this;
  619.     }
  620.  
  621.   int nr_insert = nr;
  622.   ComplexMatrix retval (nr + a.length (), nc);
  623.   retval.insert (*this, 0, 0);
  624.   retval.insert (a, nr_insert, 0);
  625.   return retval;
  626. }
  627.  
  628. ComplexMatrix
  629. ComplexMatrix::stack (const DiagMatrix& a) const
  630. {
  631.   int nr = rows ();
  632.   int nc = cols ();
  633.   if (nc != a.cols ())
  634.     {
  635.       (*current_liboctave_error_handler)
  636.     ("column dimension mismatch for stack");
  637.       return *this;
  638.     }
  639.  
  640.   int nr_insert = nr;
  641.   ComplexMatrix retval (nr + a.rows (), nc);
  642.   retval.insert (*this, 0, 0);
  643.   retval.insert (a, nr_insert, 0);
  644.   return retval;
  645. }
  646.  
  647. ComplexMatrix
  648. ComplexMatrix::stack (const ComplexMatrix& a) const
  649. {
  650.   int nr = rows ();
  651.   int nc = cols ();
  652.   if (nc != a.cols ())
  653.     {
  654.       (*current_liboctave_error_handler)
  655.     ("column dimension mismatch for stack");
  656.       return *this;
  657.     }
  658.  
  659.   int nr_insert = nr;
  660.   ComplexMatrix retval (nr + a.rows (), nc);
  661.   retval.insert (*this, 0, 0);
  662.   retval.insert (a, nr_insert, 0);
  663.   return retval;
  664. }
  665.  
  666. ComplexMatrix
  667. ComplexMatrix::stack (const ComplexRowVector& a) const
  668. {
  669.   int nr = rows ();
  670.   int nc = cols ();
  671.   if (nc != a.length ())
  672.     {
  673.       (*current_liboctave_error_handler)
  674.     ("column dimension mismatch for stack");
  675.       return *this;
  676.     }
  677.  
  678.   int nr_insert = nr;
  679.   ComplexMatrix retval (nr + 1, nc);
  680.   retval.insert (*this, 0, 0);
  681.   retval.insert (a, nr_insert, 0);
  682.   return retval;
  683. }
  684.  
  685. ComplexMatrix
  686. ComplexMatrix::stack (const ComplexColumnVector& a) const
  687. {
  688.   int nr = rows ();
  689.   int nc = cols ();
  690.   if (nc != 1)
  691.     {
  692.       (*current_liboctave_error_handler)
  693.     ("column dimension mismatch for stack");
  694.       return *this;
  695.     }
  696.  
  697.   int nr_insert = nr;
  698.   ComplexMatrix retval (nr + a.length (), nc);
  699.   retval.insert (*this, 0, 0);
  700.   retval.insert (a, nr_insert, 0);
  701.   return retval;
  702. }
  703.  
  704. ComplexMatrix
  705. ComplexMatrix::stack (const ComplexDiagMatrix& a) const
  706. {
  707.   int nr = rows ();
  708.   int nc = cols ();
  709.   if (nc != a.cols ())
  710.     {
  711.       (*current_liboctave_error_handler)
  712.     ("column dimension mismatch for stack");
  713.       return *this;
  714.     }
  715.  
  716.   int nr_insert = nr;
  717.   ComplexMatrix retval (nr + a.rows (), nc);
  718.   retval.insert (*this, 0, 0);
  719.   retval.insert (a, nr_insert, 0);
  720.   return retval;
  721. }
  722.  
  723. ComplexMatrix
  724. ComplexMatrix::hermitian (void) const
  725. {
  726.   int nr = rows ();
  727.   int nc = cols ();
  728.   ComplexMatrix result;
  729.   if (length () > 0)
  730.     {
  731.       result.resize (nc, nr);
  732.       for (int j = 0; j < nc; j++)
  733.     for (int i = 0; i < nr; i++)
  734.       result.elem (j, i) = conj (elem (i, j));
  735.     }
  736.   return result;
  737. }
  738.  
  739. ComplexMatrix
  740. conj (const ComplexMatrix& a)
  741. {
  742.   int a_len = a.length ();
  743.   ComplexMatrix retval;
  744.   if (a_len > 0)
  745.     retval = ComplexMatrix (conj_dup (a.data (), a_len), a.rows (),
  746.                 a.cols ());
  747.   return retval;
  748. }
  749.  
  750. // resize is the destructive equivalent for this one
  751.  
  752. ComplexMatrix
  753. ComplexMatrix::extract (int r1, int c1, int r2, int c2) const
  754. {
  755.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  756.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  757.  
  758.   int new_r = r2 - r1 + 1;
  759.   int new_c = c2 - c1 + 1;
  760.  
  761.   ComplexMatrix result (new_r, new_c);
  762.  
  763.   for (int j = 0; j < new_c; j++)
  764.     for (int i = 0; i < new_r; i++)
  765.       result.elem (i, j) = elem (r1+i, c1+j);
  766.  
  767.   return result;
  768. }
  769.  
  770. // extract row or column i.
  771.  
  772. ComplexRowVector
  773. ComplexMatrix::row (int i) const
  774. {
  775.   int nc = cols ();
  776.   if (i < 0 || i >= rows ())
  777.     {
  778.       (*current_liboctave_error_handler) ("invalid row selection");
  779.       return ComplexRowVector ();
  780.     }
  781.  
  782.   ComplexRowVector retval (nc);
  783.   for (int j = 0; j < cols (); j++)
  784.     retval.elem (j) = elem (i, j);
  785.  
  786.   return retval;
  787. }
  788.  
  789. ComplexRowVector
  790. ComplexMatrix::row (char *s) const
  791. {
  792.   if (! s)
  793.     {
  794.       (*current_liboctave_error_handler) ("invalid row selection");
  795.       return ComplexRowVector ();
  796.     }
  797.  
  798.   char c = *s;
  799.   if (c == 'f' || c == 'F')
  800.     return row (0);
  801.   else if (c == 'l' || c == 'L')
  802.     return row (rows () - 1);
  803.   else
  804.     {
  805.       (*current_liboctave_error_handler) ("invalid row selection");
  806.       return ComplexRowVector ();
  807.     }
  808. }
  809.  
  810. ComplexColumnVector
  811. ComplexMatrix::column (int i) const
  812. {
  813.   int nr = rows ();
  814.   if (i < 0 || i >= cols ())
  815.     {
  816.       (*current_liboctave_error_handler) ("invalid column selection");
  817.       return ComplexColumnVector ();
  818.     }
  819.  
  820.   ComplexColumnVector retval (nr);
  821.   for (int j = 0; j < nr; j++)
  822.     retval.elem (j) = elem (j, i);
  823.  
  824.   return retval;
  825. }
  826.  
  827. ComplexColumnVector
  828. ComplexMatrix::column (char *s) const
  829. {
  830.   if (! s)
  831.     {
  832.       (*current_liboctave_error_handler) ("invalid column selection");
  833.       return ComplexColumnVector ();
  834.     }
  835.  
  836.   char c = *s;
  837.   if (c == 'f' || c == 'F')
  838.     return column (0);
  839.   else if (c == 'l' || c == 'L')
  840.     return column (cols () - 1);
  841.   else
  842.     {
  843.       (*current_liboctave_error_handler) ("invalid column selection");
  844.       return ComplexColumnVector ();
  845.     }
  846. }
  847.  
  848. ComplexMatrix
  849. ComplexMatrix::inverse (void) const
  850. {
  851.   int info;
  852.   double rcond;
  853.   return inverse (info, rcond);
  854. }
  855.  
  856. ComplexMatrix
  857. ComplexMatrix::inverse (int& info) const
  858. {
  859.   double rcond;
  860.   return inverse (info, rcond);
  861. }
  862.  
  863. ComplexMatrix
  864. ComplexMatrix::inverse (int& info, double& rcond, int force) const
  865. {
  866.   ComplexMatrix retval;
  867.  
  868.   int nr = rows ();
  869.   int nc = cols ();
  870.  
  871.   if (nr != nc)
  872.     (*current_liboctave_error_handler) ("inverse requires square matrix");
  873.   else
  874.     {
  875.       info = 0;
  876.  
  877.       Array<int> ipvt (nr);
  878.       int *pipvt = ipvt.fortran_vec ();
  879.  
  880.       Array<Complex> z (nr);
  881.       Complex *pz = z.fortran_vec ();
  882.  
  883.       retval = *this;
  884.       Complex *tmp_data = retval.fortran_vec ();
  885.  
  886.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz));
  887.  
  888.       if (f77_exception_encountered)
  889.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  890.       else
  891.     {
  892.       volatile double rcond_plus_one = rcond + 1.0;
  893.  
  894.       if (rcond_plus_one == 1.0)
  895.         info = -1;
  896.  
  897.       if (info == -1 && ! force)
  898.         retval = *this;  // Restore contents.
  899.       else
  900.         {
  901.           Complex *dummy = 0;
  902.  
  903.           F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy,
  904.                        pz, 1));
  905.  
  906.           if (f77_exception_encountered)
  907.         (*current_liboctave_error_handler)
  908.           ("unrecoverable error in zgedi");
  909.         }
  910.     }
  911.     }
  912.  
  913.   return retval;
  914. }
  915.  
  916. ComplexMatrix
  917. ComplexMatrix::pseudo_inverse (double tol)
  918. {
  919.   ComplexMatrix retval;
  920.  
  921.   ComplexSVD result (*this);
  922.  
  923.   DiagMatrix S = result.singular_values ();
  924.   ComplexMatrix U = result.left_singular_matrix ();
  925.   ComplexMatrix V = result.right_singular_matrix ();
  926.  
  927.   ColumnVector sigma = S.diag ();
  928.  
  929.   int r = sigma.length () - 1;
  930.   int nr = rows ();
  931.   int nc = cols ();
  932.  
  933.   if (tol <= 0.0)
  934.     {
  935.       if (nr > nc)
  936.     tol = nr * sigma.elem (0) * DBL_EPSILON;
  937.       else
  938.     tol = nc * sigma.elem (0) * DBL_EPSILON;
  939.     }
  940.  
  941.   while (r >= 0 && sigma.elem (r) < tol)
  942.     r--;
  943.  
  944.   if (r < 0)
  945.     retval = ComplexMatrix (nc, nr, 0.0);
  946.   else
  947.     {
  948.       ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
  949.       DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
  950.       ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
  951.       retval = Vr * D * Ur.hermitian ();
  952.     }
  953.  
  954.   return retval;
  955. }
  956.  
  957. ComplexMatrix
  958. ComplexMatrix::fourier (void) const
  959. {
  960.   ComplexMatrix retval;
  961.  
  962.   int nr = rows ();
  963.   int nc = cols ();
  964.  
  965.   int npts, nsamples;
  966.  
  967.   if (nr == 1 || nc == 1)
  968.     {
  969.       npts = nr > nc ? nr : nc;
  970.       nsamples = 1;
  971.     }
  972.   else
  973.     {
  974.       npts = nr;
  975.       nsamples = nc;
  976.     }
  977.  
  978.   int nn = 4*npts+15;
  979.  
  980.   Array<Complex> wsave (nn);
  981.   Complex *pwsave = wsave.fortran_vec ();
  982.  
  983.   retval = *this;
  984.   Complex *tmp_data = retval.fortran_vec ();
  985.  
  986.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  987.  
  988.   for (int j = 0; j < nsamples; j++)
  989.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  990.  
  991.   return retval;
  992. }
  993.  
  994. ComplexMatrix
  995. ComplexMatrix::ifourier (void) const
  996. {
  997.   ComplexMatrix retval;
  998.  
  999.   int nr = rows ();
  1000.   int nc = cols ();
  1001.  
  1002.   int npts, nsamples;
  1003.  
  1004.   if (nr == 1 || nc == 1)
  1005.     {
  1006.       npts = nr > nc ? nr : nc;
  1007.       nsamples = 1;
  1008.     }
  1009.   else
  1010.     {
  1011.       npts = nr;
  1012.       nsamples = nc;
  1013.     }
  1014.  
  1015.   int nn = 4*npts+15;
  1016.  
  1017.   Array<Complex> wsave (nn);
  1018.   Complex *pwsave = wsave.fortran_vec ();
  1019.  
  1020.   retval = *this;
  1021.   Complex *tmp_data = retval.fortran_vec ();
  1022.  
  1023.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1024.  
  1025.   for (int j = 0; j < nsamples; j++)
  1026.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  1027.  
  1028.   for (int j = 0; j < npts*nsamples; j++)
  1029.     tmp_data[j] = tmp_data[j] / npts;
  1030.  
  1031.   return retval;
  1032. }
  1033.  
  1034. ComplexMatrix
  1035. ComplexMatrix::fourier2d (void) const
  1036. {
  1037.   ComplexMatrix retval;
  1038.  
  1039.   int nr = rows ();
  1040.   int nc = cols ();
  1041.  
  1042.   int npts, nsamples;
  1043.  
  1044.   if (nr == 1 || nc == 1)
  1045.     {
  1046.       npts = nr > nc ? nr : nc;
  1047.       nsamples = 1;
  1048.     }
  1049.   else
  1050.     {
  1051.       npts = nr;
  1052.       nsamples = nc;
  1053.     }
  1054.  
  1055.   int nn = 4*npts+15;
  1056.  
  1057.   Array<Complex> wsave (nn);
  1058.   Complex *pwsave = wsave.fortran_vec ();
  1059.  
  1060.   retval = *this;
  1061.   Complex *tmp_data = retval.fortran_vec ();
  1062.  
  1063.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1064.  
  1065.   for (int j = 0; j < nsamples; j++)
  1066.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  1067.  
  1068.   npts = nc;
  1069.   nsamples = nr;
  1070.   nn = 4*npts+15;
  1071.  
  1072.   wsave.resize (nn);
  1073.   pwsave = wsave.fortran_vec ();
  1074.  
  1075.   Array<Complex> row (npts);
  1076.   Complex *prow = row.fortran_vec ();
  1077.  
  1078.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1079.  
  1080.   for (int j = 0; j < nsamples; j++)
  1081.     {
  1082.       for (int i = 0; i < npts; i++)
  1083.     prow[i] = tmp_data[i*nr + j];
  1084.  
  1085.       F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);
  1086.  
  1087.       for (int i = 0; i < npts; i++)
  1088.     tmp_data[i*nr + j] = prow[i];
  1089.     }
  1090.  
  1091.   return retval;
  1092. }
  1093.  
  1094. ComplexMatrix
  1095. ComplexMatrix::ifourier2d (void) const
  1096. {
  1097.   ComplexMatrix retval;
  1098.  
  1099.   int nr = rows ();
  1100.   int nc = cols ();
  1101.  
  1102.   int npts, nsamples;
  1103.  
  1104.   if (nr == 1 || nc == 1)
  1105.     {
  1106.       npts = nr > nc ? nr : nc;
  1107.       nsamples = 1;
  1108.     }
  1109.   else
  1110.     {
  1111.       npts = nr;
  1112.       nsamples = nc;
  1113.     }
  1114.  
  1115.   int nn = 4*npts+15;
  1116.  
  1117.   Array<Complex> wsave (nn);
  1118.   Complex *pwsave = wsave.fortran_vec ();
  1119.  
  1120.   retval = *this;
  1121.   Complex *tmp_data = retval.fortran_vec ();
  1122.  
  1123.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1124.  
  1125.   for (int j = 0; j < nsamples; j++)
  1126.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  1127.  
  1128.   for (int j = 0; j < npts*nsamples; j++)
  1129.     tmp_data[j] = tmp_data[j] / npts;
  1130.  
  1131.   npts = nc;
  1132.   nsamples = nr;
  1133.   nn = 4*npts+15;
  1134.  
  1135.   wsave.resize (nn);
  1136.   pwsave = wsave.fortran_vec ();
  1137.  
  1138.   Array<Complex> row (npts);
  1139.   Complex *prow = row.fortran_vec ();
  1140.  
  1141.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  1142.  
  1143.   for (int j = 0; j < nsamples; j++)
  1144.     {
  1145.       for (int i = 0; i < npts; i++)
  1146.     prow[i] = tmp_data[i*nr + j];
  1147.  
  1148.       F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);
  1149.  
  1150.       for (int i = 0; i < npts; i++)
  1151.     tmp_data[i*nr + j] = prow[i] / npts;
  1152.     }
  1153.  
  1154.   return retval;
  1155. }
  1156.  
  1157. ComplexDET
  1158. ComplexMatrix::determinant (void) const
  1159. {
  1160.   int info;
  1161.   double rcond;
  1162.   return determinant (info, rcond);
  1163. }
  1164.  
  1165. ComplexDET
  1166. ComplexMatrix::determinant (int& info) const
  1167. {
  1168.   double rcond;
  1169.   return determinant (info, rcond);
  1170. }
  1171.  
  1172. ComplexDET
  1173. ComplexMatrix::determinant (int& info, double& rcond) const
  1174. {
  1175.   ComplexDET retval;
  1176.  
  1177.   int nr = rows ();
  1178.   int nc = cols ();
  1179.  
  1180.   if (nr == 0 || nc == 0)
  1181.     {
  1182.       Complex d[2];
  1183.       d[0] = 1.0;
  1184.       d[1] = 0.0;
  1185.       retval = ComplexDET (d);
  1186.     }
  1187.   else
  1188.     {
  1189.       info = 0;
  1190.  
  1191.       Array<int> ipvt (nr);
  1192.       int *pipvt = ipvt.fortran_vec ();
  1193.  
  1194.       Array<Complex> z (nr);
  1195.       Complex *pz = z.fortran_vec ();
  1196.  
  1197.       ComplexMatrix atmp = *this;
  1198.       Complex *tmp_data = atmp.fortran_vec ();
  1199.  
  1200.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1201.  
  1202.       if (f77_exception_encountered)
  1203.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  1204.       else
  1205.     {
  1206.       volatile double rcond_plus_one = rcond + 1.0;
  1207.  
  1208.       if (rcond_plus_one == 1.0)
  1209.         {
  1210.           info = -1;
  1211.           retval = ComplexDET ();
  1212.         }
  1213.       else
  1214.         {
  1215.           Complex d[2];
  1216.  
  1217.           F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));
  1218.  
  1219.           if (f77_exception_encountered)
  1220.         (*current_liboctave_error_handler)
  1221.           ("unrecoverable error in dgedi");
  1222.           else
  1223.         retval = ComplexDET (d);
  1224.         }
  1225.     }
  1226.     }
  1227.  
  1228.   return retval;
  1229. }
  1230.  
  1231. ComplexMatrix
  1232. ComplexMatrix::solve (const Matrix& b) const
  1233. {
  1234.   int info;
  1235.   double rcond;
  1236.   return solve (b, info, rcond);
  1237. }
  1238.  
  1239. ComplexMatrix
  1240. ComplexMatrix::solve (const Matrix& b, int& info) const
  1241. {
  1242.   double rcond;
  1243.   return solve (b, info, rcond);
  1244. }
  1245.  
  1246. ComplexMatrix
  1247. ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const
  1248. {
  1249.   ComplexMatrix tmp (b);
  1250.   return solve (tmp, info, rcond);
  1251. }
  1252.  
  1253. ComplexMatrix
  1254. ComplexMatrix::solve (const ComplexMatrix& b) const
  1255. {
  1256.   int info;
  1257.   double rcond;
  1258.   return solve (b, info, rcond);
  1259. }
  1260.  
  1261. ComplexMatrix
  1262. ComplexMatrix::solve (const ComplexMatrix& b, int& info) const
  1263. {
  1264.   double rcond;
  1265.   return solve (b, info, rcond);
  1266. }
  1267. ComplexMatrix
  1268. ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
  1269. {
  1270.   ComplexMatrix retval;
  1271.  
  1272.   int nr = rows ();
  1273.   int nc = cols ();
  1274.  
  1275.   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
  1276.     (*current_liboctave_error_handler)
  1277.       ("matrix dimension mismatch in solution of linear equations");
  1278.   else
  1279.     {
  1280.       info = 0;
  1281.  
  1282.       Array<int> ipvt (nr);
  1283.       int *pipvt = ipvt.fortran_vec ();
  1284.  
  1285.       Array<Complex> z (nr);
  1286.       Complex *pz = z.fortran_vec ();
  1287.  
  1288.       ComplexMatrix atmp = *this;
  1289.       Complex *tmp_data = atmp.fortran_vec ();
  1290.  
  1291.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1292.  
  1293.       if (f77_exception_encountered)
  1294.     (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
  1295.       else
  1296.     {
  1297.       volatile double rcond_plus_one = rcond + 1.0;
  1298.  
  1299.       if (rcond_plus_one == 1.0)
  1300.         {
  1301.           info = -2;
  1302.         }
  1303.       else
  1304.         {
  1305.           retval = b;
  1306.           Complex *result = retval.fortran_vec ();
  1307.  
  1308.           int b_nc = b.cols ();
  1309.  
  1310.           for (volatile int j = 0; j < b_nc; j++)
  1311.         {
  1312.           F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt,
  1313.                        &result[nr*j], 0));
  1314.  
  1315.           if (f77_exception_encountered)
  1316.             {
  1317.               (*current_liboctave_error_handler)
  1318.             ("unrecoverable error in dgesl");
  1319.  
  1320.               break;
  1321.             }
  1322.         }
  1323.         }
  1324.     }
  1325.     }
  1326.  
  1327.   return retval;
  1328. }
  1329.  
  1330. ComplexColumnVector
  1331. ComplexMatrix::solve (const ComplexColumnVector& b) const
  1332. {
  1333.   int info;
  1334.   double rcond;
  1335.   return solve (b, info, rcond);
  1336. }
  1337.  
  1338. ComplexColumnVector
  1339. ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
  1340. {
  1341.   double rcond;
  1342.   return solve (b, info, rcond);
  1343. }
  1344.  
  1345. ComplexColumnVector
  1346. ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
  1347.               double& rcond) const
  1348. {
  1349.   ComplexColumnVector retval;
  1350.  
  1351.   int nr = rows ();
  1352.   int nc = cols ();
  1353.  
  1354.   if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
  1355.     (*current_liboctave_error_handler)
  1356.       ("matrix dimension mismatch in solution of linear equations");
  1357.   else
  1358.     {
  1359.       info = 0;
  1360.  
  1361.       Array<int> ipvt (nr);
  1362.       int *pipvt = ipvt.fortran_vec ();
  1363.  
  1364.       Array<Complex> z (nr);
  1365.       Complex *pz = z.fortran_vec ();
  1366.  
  1367.       ComplexMatrix atmp = *this;
  1368.       Complex *tmp_data = atmp.fortran_vec ();
  1369.  
  1370.       F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1371.  
  1372.       if (f77_exception_encountered)
  1373.     (*current_liboctave_error_handler)
  1374.       ("unrecoverable error in dgeco");
  1375.       else
  1376.     {
  1377.       volatile double rcond_plus_one = rcond + 1.0;
  1378.  
  1379.       if (rcond_plus_one == 1.0)
  1380.         {
  1381.           info = -2;
  1382.         }
  1383.       else
  1384.         {
  1385.           retval = b;
  1386.           Complex *result = retval.fortran_vec ();
  1387.  
  1388.           F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0));
  1389.  
  1390.           if (f77_exception_encountered)
  1391.         (*current_liboctave_error_handler)
  1392.           ("unrecoverable error in dgesl");
  1393.         }
  1394.     }
  1395.     }
  1396.  
  1397.   return retval;
  1398. }
  1399.  
  1400. ComplexMatrix
  1401. ComplexMatrix::lssolve (const ComplexMatrix& b) const
  1402. {
  1403.   int info;
  1404.   int rank;
  1405.   return lssolve (b, info, rank);
  1406. }
  1407.  
  1408. ComplexMatrix
  1409. ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const
  1410. {
  1411.   int rank;
  1412.   return lssolve (b, info, rank);
  1413. }
  1414.  
  1415. ComplexMatrix
  1416. ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
  1417. {
  1418.   ComplexMatrix retval;
  1419.  
  1420.   int nrhs = b.cols ();
  1421.  
  1422.   int m = rows ();
  1423.   int n = cols ();
  1424.  
  1425.   if (m == 0 || n == 0 || m != b.rows ())
  1426.     (*current_liboctave_error_handler)
  1427.       ("matrix dimension mismatch solution of linear equations");
  1428.   else
  1429.     {
  1430.       ComplexMatrix atmp = *this;
  1431.       Complex *tmp_data = atmp.fortran_vec ();
  1432.  
  1433.       int nrr = m > n ? m : n;
  1434.       ComplexMatrix result (nrr, nrhs);
  1435.  
  1436.       for (int j = 0; j < nrhs; j++)
  1437.     for (int i = 0; i < m; i++)
  1438.       result.elem (i, j) = b.elem (i, j);
  1439.  
  1440.       Complex *presult = result.fortran_vec ();
  1441.  
  1442.       int len_s = m < n ? m : n;
  1443.       Array<double> s (len_s);
  1444.       double *ps = s.fortran_vec ();
  1445.  
  1446.       double rcond = -1.0;
  1447.  
  1448.       int lwork;
  1449.       if (m < n)
  1450.     lwork = 2*m + (nrhs > n ? nrhs : n);
  1451.       else
  1452.     lwork = 2*n + (nrhs > m ? nrhs : m);
  1453.  
  1454.       lwork *= 16;
  1455.  
  1456.       Array<Complex> work (lwork);
  1457.       Complex *pwork = work.fortran_vec ();
  1458.  
  1459.       int lrwork = (5 * (m < n ? m : n)) - 4;
  1460.       lrwork = lrwork > 1 ? lrwork : 1;
  1461.       Array<double> rwork (lrwork);
  1462.       double *prwork = rwork.fortran_vec ();
  1463.  
  1464.       F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
  1465.                  nrr, ps, rcond, rank, pwork, lwork,
  1466.                  prwork, info));
  1467.  
  1468.       if (f77_exception_encountered)
  1469.     (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
  1470.       else
  1471.     {
  1472.       retval.resize (n, nrhs);
  1473.       for (int j = 0; j < nrhs; j++)
  1474.         for (int i = 0; i < n; i++)
  1475.           retval.elem (i, j) = result.elem (i, j);
  1476.     }
  1477.     }
  1478.  
  1479.   return retval;
  1480. }
  1481.  
  1482. ComplexColumnVector
  1483. ComplexMatrix::lssolve (const ComplexColumnVector& b) const
  1484. {
  1485.   int info;
  1486.   int rank;
  1487.   return lssolve (b, info, rank);
  1488. }
  1489.  
  1490. ComplexColumnVector
  1491. ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const
  1492. {
  1493.   int rank;
  1494.   return lssolve (b, info, rank);
  1495. }
  1496.  
  1497. ComplexColumnVector
  1498. ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
  1499.             int& rank) const
  1500. {
  1501.   ComplexColumnVector retval;
  1502.  
  1503.   int nrhs = 1;
  1504.  
  1505.   int m = rows ();
  1506.   int n = cols ();
  1507.  
  1508.   if (m == 0 || n == 0 || m != b.length ())
  1509.     (*current_liboctave_error_handler)
  1510.       ("matrix dimension mismatch solution of least squares problem");
  1511.   else
  1512.     {
  1513.       ComplexMatrix atmp = *this;
  1514.       Complex *tmp_data = atmp.fortran_vec ();
  1515.  
  1516.       int nrr = m > n ? m : n;
  1517.       ComplexColumnVector result (nrr);
  1518.  
  1519.       for (int i = 0; i < m; i++)
  1520.     result.elem (i) = b.elem (i);
  1521.  
  1522.       Complex *presult = result.fortran_vec ();
  1523.  
  1524.       int len_s = m < n ? m : n;
  1525.       Array<double> s (len_s);
  1526.       double *ps = s.fortran_vec ();
  1527.  
  1528.       double rcond = -1.0;
  1529.  
  1530.       int lwork;
  1531.       if (m < n)
  1532.     lwork = 2*m + (nrhs > n ? nrhs : n);
  1533.       else
  1534.     lwork = 2*n + (nrhs > m ? nrhs : m);
  1535.  
  1536.       lwork *= 16;
  1537.  
  1538.       Array<Complex> work (lwork);
  1539.       Complex *pwork = work.fortran_vec ();
  1540.  
  1541.       int lrwork = (5 * (m < n ? m : n)) - 4;
  1542.       lrwork = lrwork > 1 ? lrwork : 1;
  1543.       Array<double> rwork (lrwork);
  1544.       double *prwork = rwork.fortran_vec ();
  1545.  
  1546.       F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
  1547.                  nrr, ps, rcond, rank, pwork, lwork,
  1548.                  prwork, info));
  1549.  
  1550.       if (f77_exception_encountered)
  1551.     (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
  1552.       else
  1553.     {
  1554.       retval.resize (n);
  1555.       for (int i = 0; i < n; i++)
  1556.         retval.elem (i) = result.elem (i);
  1557.     }
  1558.     }
  1559.  
  1560.   return retval;
  1561. }
  1562.  
  1563. // Constants for matrix exponential calculation.
  1564.  
  1565. static double padec [] =
  1566. {
  1567.   5.0000000000000000e-1,
  1568.   1.1666666666666667e-1,
  1569.   1.6666666666666667e-2,
  1570.   1.6025641025641026e-3,
  1571.   1.0683760683760684e-4,
  1572.   4.8562548562548563e-6,
  1573.   1.3875013875013875e-7,
  1574.   1.9270852604185938e-9,
  1575. };
  1576.  
  1577. ComplexMatrix
  1578. ComplexMatrix::expm (void) const
  1579. {
  1580.   ComplexMatrix retval;
  1581.  
  1582.   ComplexMatrix m = *this;
  1583.  
  1584.   int nc = columns ();
  1585.  
  1586.   // Preconditioning step 1: trace normalization to reduce dynamic
  1587.   // range of poles, but avoid making stable eigenvalues unstable.
  1588.  
  1589.   // trace shift value
  1590.   Complex trshift = 0.0;
  1591.  
  1592.   for (int i = 0; i < nc; i++)
  1593.     trshift += m.elem (i, i);
  1594.  
  1595.   trshift /= nc;
  1596.  
  1597.   if (trshift.real () < 0.0)
  1598.     trshift = trshift.imag ();
  1599.  
  1600.   for (int i = 0; i < nc; i++)
  1601.     m.elem (i, i) -= trshift;
  1602.  
  1603.   // Preconditioning step 2: eigenvalue balancing.
  1604.   // code follows development in AEPBAL
  1605.  
  1606.   Complex *mp = m.fortran_vec ();
  1607.   int ilo, ihi, info;
  1608.  
  1609.   // FIXME: should pass job as a parameter in expm
  1610.   char job = 'B';
  1611.  
  1612.   Array<double> scale(nc);
  1613.   double *pscale = scale.fortran_vec ();
  1614.  
  1615.   F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi, pscale, info,
  1616.                  1L, 1L));
  1617.  
  1618.   if (f77_exception_encountered)
  1619.     {
  1620.       (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
  1621.       return retval;
  1622.     }
  1623.  
  1624.   // construct balancing matrices dmat, dinv
  1625.  
  1626.   Matrix dmat = Matrix (nc, nc, 0.0);
  1627.   Matrix dinv = Matrix (nc, nc, 0.0);
  1628.  
  1629.   for (int i = 0; i < nc; i++)
  1630.     dmat(i,i) = dinv(i,i) = 1.0;
  1631.  
  1632.   // pscale contains real data, so just use dgebak, dside=R => dmat := D*dmat
  1633.   char dside = 'R';
  1634.   F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
  1635.                  dmat.fortran_vec(), nc, info, 1L, 1L));
  1636.  
  1637.   if (f77_exception_encountered)
  1638.     {
  1639.       (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
  1640.       return retval;
  1641.     }
  1642.  
  1643.   // dgebak, dside=L => dinv := dinv*D^{-1}
  1644.   dside = 'L';
  1645.   F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
  1646.                  dinv.fortran_vec(), nc, info, 1L, 1L));
  1647.  
  1648.   if (f77_exception_encountered)
  1649.     {
  1650.       (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
  1651.       return retval;
  1652.     }
  1653.  
  1654.   // Preconditioning step 3: scaling.
  1655.  
  1656.   ColumnVector work (nc);
  1657.   double inf_norm;
  1658.  
  1659.   F77_XFCN (xzlange, XZLANGE, ("I", nc, nc, m.fortran_vec (), nc,
  1660.                    work.fortran_vec (), inf_norm));
  1661.  
  1662.   if (f77_exception_encountered)
  1663.     {
  1664.       (*current_liboctave_error_handler) ("unrecoverable error in zlange");
  1665.       return retval;
  1666.     }
  1667.  
  1668.   int sqpow = (inf_norm > 0.0
  1669.            ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
  1670.  
  1671.   // Check whether we need to square at all.
  1672.  
  1673.   if (sqpow < 0)
  1674.     sqpow = 0;
  1675.  
  1676.   if (sqpow > 0)
  1677.     {
  1678.       double scale_factor = 1.0;
  1679.       for (int i = 0; i < sqpow; i++)
  1680.     scale_factor *= 2.0;
  1681.  
  1682.       m = m / scale_factor;
  1683.     }
  1684.  
  1685.   // npp, dpp: pade' approx polynomial matrices.
  1686.  
  1687.   ComplexMatrix npp (nc, nc, 0.0);
  1688.   ComplexMatrix dpp = npp;
  1689.  
  1690.   // Now powers a^8 ... a^1.
  1691.  
  1692.   int minus_one_j = -1;
  1693.   for (int j = 7; j >= 0; j--)
  1694.     {
  1695.       npp = m * npp + m * padec[j];
  1696.       dpp = m * dpp + m * (minus_one_j * padec[j]);
  1697.       minus_one_j *= -1;
  1698.     }
  1699.  
  1700.   // Zero power.
  1701.  
  1702.   dpp = -dpp;
  1703.   for (int j = 0; j < nc; j++)
  1704.     {
  1705.       npp.elem (j, j) += 1.0;
  1706.       dpp.elem (j, j) += 1.0;
  1707.     }
  1708.  
  1709.   // Compute pade approximation = inverse (dpp) * npp.
  1710.  
  1711.   retval = dpp.solve (npp);
  1712.     
  1713.   // Reverse preconditioning step 3: repeated squaring.
  1714.  
  1715.   while (sqpow)
  1716.     {
  1717.       retval = retval * retval;
  1718.       sqpow--;
  1719.     }
  1720.  
  1721.   // Reverse preconditioning step 2: inverse balancing.
  1722.   // FIXME: need to figure out how to get dgebak to do
  1723.   // dmat*retval*dinv instead of dinv*retval*dmat
  1724.   // This works for now
  1725.   retval = dmat*retval*dinv;
  1726.  
  1727.   // Reverse preconditioning step 1: fix trace normalization.
  1728.  
  1729.   return exp (trshift) * retval;
  1730. }
  1731.  
  1732. // column vector by row vector -> matrix operations
  1733.  
  1734. ComplexMatrix
  1735. operator * (const ColumnVector& v, const ComplexRowVector& a)
  1736. {
  1737.   ComplexColumnVector tmp (v);
  1738.   return tmp * a;
  1739. }
  1740.  
  1741. ComplexMatrix
  1742. operator * (const ComplexColumnVector& a, const RowVector& b)
  1743. {
  1744.   ComplexRowVector tmp (b);
  1745.   return a * tmp;
  1746. }
  1747.  
  1748. ComplexMatrix
  1749. operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
  1750. {
  1751.   ComplexMatrix retval;
  1752.  
  1753.   int len = v.length ();
  1754.  
  1755.   if (len != 0)
  1756.     {
  1757.       int a_len = a.length ();
  1758.  
  1759.       retval.resize (len, a_len);
  1760.       Complex *c = retval.fortran_vec ();
  1761.  
  1762.       F77_XFCN (zgemm, ZGEMM, ("N", "N", len, a_len, 1, 1.0,
  1763.                    v.data (), len, a.data (), 1, 0.0,
  1764.                    c, len, 1L, 1L)); 
  1765.  
  1766.       if (f77_exception_encountered)
  1767.     (*current_liboctave_error_handler)
  1768.       ("unrecoverable error in zgemm");
  1769.     }
  1770.  
  1771.   return retval;
  1772. }
  1773.  
  1774. // matrix by diagonal matrix -> matrix operations
  1775.  
  1776. ComplexMatrix&
  1777. ComplexMatrix::operator += (const DiagMatrix& a)
  1778. {
  1779.   int nr = rows ();
  1780.   int nc = cols ();
  1781.  
  1782.   int a_nr = rows ();
  1783.   int a_nc = cols ();
  1784.  
  1785.   if (nr != a_nr || nc != a_nc)
  1786.     {
  1787.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1788.       return *this;
  1789.     }
  1790.  
  1791.   for (int i = 0; i < a.length (); i++)
  1792.     elem (i, i) += a.elem (i, i);
  1793.  
  1794.   return *this;
  1795. }
  1796.  
  1797. ComplexMatrix&
  1798. ComplexMatrix::operator -= (const DiagMatrix& a)
  1799. {
  1800.   int nr = rows ();
  1801.   int nc = cols ();
  1802.  
  1803.   int a_nr = rows ();
  1804.   int a_nc = cols ();
  1805.  
  1806.   if (nr != a_nr || nc != a_nc)
  1807.     {
  1808.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1809.       return *this;
  1810.     }
  1811.  
  1812.   for (int i = 0; i < a.length (); i++)
  1813.     elem (i, i) -= a.elem (i, i);
  1814.  
  1815.   return *this;
  1816. }
  1817.  
  1818. ComplexMatrix&
  1819. ComplexMatrix::operator += (const ComplexDiagMatrix& a)
  1820. {
  1821.   int nr = rows ();
  1822.   int nc = cols ();
  1823.  
  1824.   int a_nr = rows ();
  1825.   int a_nc = cols ();
  1826.  
  1827.   if (nr != a_nr || nc != a_nc)
  1828.     {
  1829.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1830.       return *this;
  1831.     }
  1832.  
  1833.   for (int i = 0; i < a.length (); i++)
  1834.     elem (i, i) += a.elem (i, i);
  1835.  
  1836.   return *this;
  1837. }
  1838.  
  1839. ComplexMatrix&
  1840. ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
  1841. {
  1842.   int nr = rows ();
  1843.   int nc = cols ();
  1844.  
  1845.   int a_nr = rows ();
  1846.   int a_nc = cols ();
  1847.  
  1848.   if (nr != a_nr || nc != a_nc)
  1849.     {
  1850.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1851.       return *this;
  1852.     }
  1853.  
  1854.   for (int i = 0; i < a.length (); i++)
  1855.     elem (i, i) -= a.elem (i, i);
  1856.  
  1857.   return *this;
  1858. }
  1859.  
  1860. // matrix by matrix -> matrix operations
  1861.  
  1862. ComplexMatrix&
  1863. ComplexMatrix::operator += (const Matrix& a)
  1864. {
  1865.   int nr = rows ();
  1866.   int nc = cols ();
  1867.  
  1868.   int a_nr = a.rows ();
  1869.   int a_nc = a.cols ();
  1870.  
  1871.   if (nr != a_nr || nc != a_nc)
  1872.     {
  1873.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1874.       return *this;
  1875.     }
  1876.  
  1877.   if (nr == 0 || nc == 0)
  1878.     return *this;
  1879.  
  1880.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  1881.  
  1882.   add2 (d, a.data (), length ());
  1883.   return *this;
  1884. }
  1885.  
  1886. ComplexMatrix&
  1887. ComplexMatrix::operator -= (const Matrix& a)
  1888. {
  1889.   int nr = rows ();
  1890.   int nc = cols ();
  1891.  
  1892.   int a_nr = a.rows ();
  1893.   int a_nc = a.cols ();
  1894.  
  1895.   if (nr != a_nr || nc != a_nc)
  1896.     {
  1897.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1898.       return *this;
  1899.     }
  1900.  
  1901.   if (nr == 0 || nc == 0)
  1902.     return *this;
  1903.  
  1904.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  1905.  
  1906.   subtract2 (d, a.data (), length ());
  1907.   return *this;
  1908. }
  1909.  
  1910. ComplexMatrix&
  1911. ComplexMatrix::operator += (const ComplexMatrix& a)
  1912. {
  1913.   int nr = rows ();
  1914.   int nc = cols ();
  1915.  
  1916.   int a_nr = a.rows ();
  1917.   int a_nc = a.cols ();
  1918.  
  1919.   if (nr != a_nr || nc != a_nc)
  1920.     {
  1921.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1922.       return *this;
  1923.     }
  1924.  
  1925.   if (nr == 0 || nc == 0)
  1926.     return *this;
  1927.  
  1928.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  1929.  
  1930.   add2 (d, a.data (), length ());
  1931.   return *this;
  1932. }
  1933.  
  1934. ComplexMatrix&
  1935. ComplexMatrix::operator -= (const ComplexMatrix& a)
  1936. {
  1937.   int nr = rows ();
  1938.   int nc = cols ();
  1939.  
  1940.   int a_nr = a.rows ();
  1941.   int a_nc = a.cols ();
  1942.  
  1943.   if (nr != a_nr || nc != a_nc)
  1944.     {
  1945.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1946.       return *this;
  1947.     }
  1948.  
  1949.   if (nr == 0 || nc == 0)
  1950.     return *this;
  1951.  
  1952.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  1953.  
  1954.   subtract2 (d, a.data (), length ());
  1955.   return *this;
  1956. }
  1957.  
  1958. // unary operations
  1959.  
  1960. boolMatrix
  1961. ComplexMatrix::operator ! (void) const
  1962. {
  1963.   int nr = rows ();
  1964.   int nc = cols ();
  1965.  
  1966.   boolMatrix b (nr, nc);
  1967.  
  1968.   for (int j = 0; j < nc; j++)
  1969.     for (int i = 0; i < nr; i++)
  1970.       b.elem (i, j) = elem (i, j) != 0.0;
  1971.  
  1972.   return b;
  1973. }
  1974.  
  1975. // other operations
  1976.  
  1977. ComplexMatrix
  1978. ComplexMatrix::map (c_c_Mapper f) const
  1979. {
  1980.   ComplexMatrix b (*this);
  1981.   return b.apply (f);
  1982. }
  1983.  
  1984. Matrix
  1985. ComplexMatrix::map (d_c_Mapper f) const
  1986. {
  1987.   int nr = rows ();
  1988.   int nc = cols ();
  1989.  
  1990.   Matrix retval (nr, nc);
  1991.  
  1992.   for (int j = 0; j < nc; j++)
  1993.     for (int i = 0; i < nr; i++)
  1994.       retval(i,j) = f (elem(i,j));
  1995.  
  1996.   return retval;
  1997. }
  1998.  
  1999. boolMatrix
  2000. ComplexMatrix::map (b_c_Mapper f) const
  2001. {
  2002.   int nr = rows ();
  2003.   int nc = cols ();
  2004.  
  2005.   boolMatrix retval (nr, nc);
  2006.  
  2007.   for (int j = 0; j < nc; j++)
  2008.     for (int i = 0; i < nr; i++)
  2009.       retval(i,j) = f (elem(i,j));
  2010.  
  2011.   return retval;
  2012. }
  2013.  
  2014. ComplexMatrix&
  2015. ComplexMatrix::apply (c_c_Mapper f)
  2016. {
  2017.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  2018.  
  2019.   for (int i = 0; i < length (); i++)
  2020.     d[i] = f (d[i]);
  2021.  
  2022.   return *this;
  2023. }
  2024.  
  2025. bool
  2026. ComplexMatrix::any_element_is_inf_or_nan (void) const
  2027. {
  2028.   int nr = rows ();
  2029.   int nc = cols ();
  2030.  
  2031.   for (int j = 0; j < nc; j++)
  2032.     for (int i = 0; i < nr; i++)
  2033.       {
  2034.     Complex val = elem (i, j);
  2035.     if (xisinf (val) || xisnan (val))
  2036.       return true;
  2037.       }
  2038.  
  2039.   return false;
  2040. }
  2041.  
  2042. // Return true if no elements have imaginary components.
  2043.  
  2044. bool
  2045. ComplexMatrix::all_elements_are_real (void) const
  2046. {
  2047.   int nr = rows ();
  2048.   int nc = cols ();
  2049.  
  2050.   for (int j = 0; j < nc; j++)
  2051.     for (int i = 0; i < nr; i++)
  2052.       if (imag (elem (i, j)) != 0.0)
  2053.     return false;
  2054.  
  2055.   return true;
  2056. }
  2057.  
  2058. // Return nonzero if any element of CM has a non-integer real or
  2059. // imaginary part.  Also extract the largest and smallest (real or
  2060. // imaginary) values and return them in MAX_VAL and MIN_VAL. 
  2061.  
  2062. bool
  2063. ComplexMatrix::all_integers (double& max_val, double& min_val) const
  2064. {
  2065.   int nr = rows ();
  2066.   int nc = cols ();
  2067.  
  2068.   if (nr > 0 && nc > 0)
  2069.     {
  2070.       Complex val = elem (0, 0);
  2071.  
  2072.       double r_val = real (val);
  2073.       double i_val = imag (val);
  2074.  
  2075.       max_val = r_val;
  2076.       min_val = r_val;
  2077.  
  2078.       if (i_val > max_val)
  2079.     max_val = i_val;
  2080.  
  2081.       if (i_val < max_val)
  2082.     min_val = i_val;
  2083.     }
  2084.   else
  2085.     return false;
  2086.  
  2087.   for (int j = 0; j < nc; j++)
  2088.     for (int i = 0; i < nr; i++)
  2089.       {
  2090.     Complex val = elem (i, j);
  2091.  
  2092.     double r_val = real (val);
  2093.     double i_val = imag (val);
  2094.  
  2095.     if (r_val > max_val)
  2096.       max_val = r_val;
  2097.  
  2098.     if (i_val > max_val)
  2099.       max_val = i_val;
  2100.  
  2101.     if (r_val < min_val)
  2102.       min_val = r_val;
  2103.  
  2104.     if (i_val < min_val)
  2105.       min_val = i_val;
  2106.  
  2107.     if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
  2108.       return false;
  2109.       }
  2110.  
  2111.   return true;
  2112. }
  2113.  
  2114. bool
  2115. ComplexMatrix::too_large_for_float (void) const
  2116. {
  2117.   int nr = rows ();
  2118.   int nc = cols ();
  2119.  
  2120.   for (int j = 0; j < nc; j++)
  2121.     for (int i = 0; i < nr; i++)
  2122.       {
  2123.     Complex val = elem (i, j);
  2124.  
  2125.     double r_val = real (val);
  2126.     double i_val = imag (val);
  2127.  
  2128.     if (r_val > FLT_MAX
  2129.         || i_val > FLT_MAX
  2130.         || r_val < FLT_MIN
  2131.         || i_val < FLT_MIN)
  2132.       return true;
  2133.       }
  2134.  
  2135.   return false;
  2136. }
  2137.  
  2138. boolMatrix
  2139. ComplexMatrix::all (void) const
  2140. {
  2141.   int nr = rows ();
  2142.   int nc = cols ();
  2143.   boolMatrix retval;
  2144.   if (nr > 0 && nc > 0)
  2145.     {
  2146.       if (nr == 1)
  2147.     {
  2148.       retval.resize (1, 1);
  2149.       retval.elem (0, 0) = true;
  2150.       for (int j = 0; j < nc; j++)
  2151.         {
  2152.           if (elem (0, j) == 0.0)
  2153.         {
  2154.           retval.elem (0, 0) = false;
  2155.           break;
  2156.         }
  2157.         }
  2158.     }
  2159.       else if (nc == 1)
  2160.     {
  2161.       retval.resize (1, 1);
  2162.       retval.elem (0, 0) = true;
  2163.       for (int i = 0; i < nr; i++)
  2164.         {
  2165.           if (elem (i, 0) == 0.0)
  2166.         {
  2167.           retval.elem (0, 0) = false;
  2168.           break;
  2169.         }
  2170.         }
  2171.     }
  2172.       else
  2173.     {
  2174.       retval.resize (1, nc);
  2175.       for (int j = 0; j < nc; j++)
  2176.         {
  2177.           retval.elem (0, j) = true;
  2178.           for (int i = 0; i < nr; i++)
  2179.         {
  2180.           if (elem (i, j) == 0.0)
  2181.             {
  2182.               retval.elem (0, j) = false;
  2183.               break;
  2184.             }
  2185.         }
  2186.         }
  2187.     }
  2188.     }
  2189.   return retval;
  2190. }
  2191.  
  2192. boolMatrix
  2193. ComplexMatrix::any (void) const
  2194. {
  2195.   int nr = rows ();
  2196.   int nc = cols ();
  2197.   boolMatrix retval;
  2198.   if (nr > 0 && nc > 0)
  2199.     {
  2200.       if (nr == 1)
  2201.     {
  2202.       retval.resize (1, 1);
  2203.       retval.elem (0, 0) = false;
  2204.       for (int j = 0; j < nc; j++)
  2205.         {
  2206.           if (elem (0, j) != 0.0)
  2207.         {
  2208.           retval.elem (0, 0) = true;
  2209.           break;
  2210.         }
  2211.         }
  2212.     }
  2213.       else if (nc == 1)
  2214.     {
  2215.       retval.resize (1, 1);
  2216.       retval.elem (0, 0) = false;
  2217.       for (int i = 0; i < nr; i++)
  2218.         {
  2219.           if (elem (i, 0) != 0.0)
  2220.         {
  2221.           retval.elem (0, 0) = true;
  2222.           break;
  2223.         }
  2224.         }
  2225.     }
  2226.       else
  2227.     {
  2228.       retval.resize (1, nc);
  2229.       for (int j = 0; j < nc; j++)
  2230.         {
  2231.           retval.elem (0, j) = false;
  2232.           for (int i = 0; i < nr; i++)
  2233.         {
  2234.           if (elem (i, j) != 0.0)
  2235.             {
  2236.               retval.elem (0, j) = true;
  2237.               break;
  2238.             }
  2239.         }
  2240.         }
  2241.     }
  2242.     }
  2243.   return retval;
  2244. }
  2245.  
  2246. ComplexMatrix
  2247. ComplexMatrix::cumprod (void) const
  2248. {
  2249.   int nr = rows ();
  2250.   int nc = cols ();
  2251.   ComplexMatrix retval;
  2252.   if (nr > 0 && nc > 0)
  2253.     {
  2254.       if (nr == 1)
  2255.     {
  2256.       retval.resize (1, nc);
  2257.       Complex prod = elem (0, 0);
  2258.       for (int j = 0; j < nc; j++)
  2259.         {
  2260.           retval.elem (0, j) = prod;
  2261.           if (j < nc - 1)
  2262.         prod *= elem (0, j+1);
  2263.         }
  2264.     }
  2265.       else if (nc == 1)
  2266.     {
  2267.       retval.resize (nr, 1);
  2268.       Complex prod = elem (0, 0);
  2269.       for (int i = 0; i < nr; i++)
  2270.         {
  2271.           retval.elem (i, 0) = prod;
  2272.           if (i < nr - 1)
  2273.         prod *= elem (i+1, 0);
  2274.         }
  2275.     }
  2276.       else
  2277.     {
  2278.       retval.resize (nr, nc);
  2279.       for (int j = 0; j < nc; j++)
  2280.         {
  2281.           Complex prod = elem (0, j);
  2282.           for (int i = 0; i < nr; i++)
  2283.         {
  2284.           retval.elem (i, j) = prod;
  2285.           if (i < nr - 1)
  2286.             prod *= elem (i+1, j);
  2287.         }
  2288.         }
  2289.     }
  2290.     }
  2291.   return retval;
  2292. }
  2293.  
  2294. ComplexMatrix
  2295. ComplexMatrix::cumsum (void) const
  2296. {
  2297.   int nr = rows ();
  2298.   int nc = cols ();
  2299.   ComplexMatrix retval;
  2300.   if (nr > 0 && nc > 0)
  2301.     {
  2302.       if (nr == 1)
  2303.     {
  2304.       retval.resize (1, nc);
  2305.       Complex sum = elem (0, 0);
  2306.       for (int j = 0; j < nc; j++)
  2307.         {
  2308.           retval.elem (0, j) = sum;
  2309.           if (j < nc - 1)
  2310.         sum += elem (0, j+1);
  2311.         }
  2312.     }
  2313.       else if (nc == 1)
  2314.     {
  2315.       retval.resize (nr, 1);
  2316.       Complex sum = elem (0, 0);
  2317.       for (int i = 0; i < nr; i++)
  2318.         {
  2319.           retval.elem (i, 0) = sum;
  2320.           if (i < nr - 1)
  2321.         sum += elem (i+1, 0);
  2322.         }
  2323.     }
  2324.       else
  2325.     {
  2326.       retval.resize (nr, nc);
  2327.       for (int j = 0; j < nc; j++)
  2328.         {
  2329.           Complex sum = elem (0, j);
  2330.           for (int i = 0; i < nr; i++)
  2331.         {
  2332.           retval.elem (i, j) = sum;
  2333.           if (i < nr - 1)
  2334.             sum += elem (i+1, j);
  2335.         }
  2336.         }
  2337.     }
  2338.     }
  2339.   return retval;
  2340. }
  2341.  
  2342. ComplexMatrix
  2343. ComplexMatrix::prod (void) const
  2344. {
  2345.   int nr = rows ();
  2346.   int nc = cols ();
  2347.   ComplexMatrix retval;
  2348.   if (nr > 0 && nc > 0)
  2349.     {
  2350.       if (nr == 1)
  2351.     {
  2352.       retval.resize (1, 1);
  2353.       retval.elem (0, 0) = 1.0;
  2354.       for (int j = 0; j < nc; j++)
  2355.         retval.elem (0, 0) *= elem (0, j);
  2356.     }
  2357.       else if (nc == 1)
  2358.     {
  2359.       retval.resize (1, 1);
  2360.       retval.elem (0, 0) = 1.0;
  2361.       for (int i = 0; i < nr; i++)
  2362.         retval.elem (0, 0) *= elem (i, 0);
  2363.     }
  2364.       else
  2365.     {
  2366.       retval.resize (1, nc);
  2367.       for (int j = 0; j < nc; j++)
  2368.         {
  2369.           retval.elem (0, j) = 1.0;
  2370.           for (int i = 0; i < nr; i++)
  2371.         retval.elem (0, j) *= elem (i, j);
  2372.         }
  2373.     }
  2374.     }
  2375.   return retval;
  2376. }
  2377.  
  2378. ComplexMatrix
  2379. ComplexMatrix::sum (void) const
  2380. {
  2381.   int nr = rows ();
  2382.   int nc = cols ();
  2383.   ComplexMatrix retval;
  2384.   if (nr > 0 && nc > 0)
  2385.     {
  2386.       if (nr == 1)
  2387.     {
  2388.       retval.resize (1, 1);
  2389.       retval.elem (0, 0) = 0.0;
  2390.       for (int j = 0; j < nc; j++)
  2391.         retval.elem (0, 0) += elem (0, j);
  2392.     }
  2393.       else if (nc == 1)
  2394.     {
  2395.       retval.resize (1, 1);
  2396.       retval.elem (0, 0) = 0.0;
  2397.       for (int i = 0; i < nr; i++)
  2398.         retval.elem (0, 0) += elem (i, 0);
  2399.     }
  2400.       else
  2401.     {
  2402.       retval.resize (1, nc);
  2403.       for (int j = 0; j < nc; j++)
  2404.         {
  2405.           retval.elem (0, j) = 0.0;
  2406.           for (int i = 0; i < nr; i++)
  2407.         retval.elem (0, j) += elem (i, j);
  2408.         }
  2409.     }
  2410.     }
  2411.   return retval;
  2412. }
  2413.  
  2414. ComplexMatrix
  2415. ComplexMatrix::sumsq (void) const
  2416. {
  2417.   int nr = rows ();
  2418.   int nc = cols ();
  2419.   ComplexMatrix retval;
  2420.   if (nr > 0 && nc > 0)
  2421.     {
  2422.       if (nr == 1)
  2423.     {
  2424.       retval.resize (1, 1);
  2425.       retval.elem (0, 0) = 0.0;
  2426.       for (int j = 0; j < nc; j++)
  2427.         {
  2428.           Complex d = elem (0, j);
  2429.           retval.elem (0, 0) += d * conj (d);
  2430.         }
  2431.     }
  2432.       else if (nc == 1)
  2433.     {
  2434.       retval.resize (1, 1);
  2435.       retval.elem (0, 0) = 0.0;
  2436.       for (int i = 0; i < nr; i++)
  2437.         {
  2438.           Complex d = elem (i, 0);
  2439.           retval.elem (0, 0) += d * conj (d);
  2440.         }
  2441.     }
  2442.       else
  2443.     {
  2444.       retval.resize (1, nc);
  2445.       for (int j = 0; j < nc; j++)
  2446.         {
  2447.           retval.elem (0, j) = 0.0;
  2448.           for (int i = 0; i < nr; i++)
  2449.         {
  2450.           Complex d = elem (i, j);
  2451.           retval.elem (0, j) += d * conj (d);
  2452.         }
  2453.         }
  2454.     }
  2455.     }
  2456.   return retval;
  2457. }
  2458.  
  2459. ComplexColumnVector
  2460. ComplexMatrix::diag (void) const
  2461. {
  2462.   return diag (0);
  2463. }
  2464.  
  2465. ComplexColumnVector
  2466. ComplexMatrix::diag (int k) const
  2467. {
  2468.   int nnr = rows ();
  2469.   int nnc = cols ();
  2470.   if (k > 0)
  2471.     nnc -= k;
  2472.   else if (k < 0)
  2473.     nnr += k;
  2474.  
  2475.   ComplexColumnVector d;
  2476.  
  2477.   if (nnr > 0 && nnc > 0)
  2478.     {
  2479.       int ndiag = (nnr < nnc) ? nnr : nnc;
  2480.  
  2481.       d.resize (ndiag);
  2482.  
  2483.       if (k > 0)
  2484.     {
  2485.       for (int i = 0; i < ndiag; i++)
  2486.         d.elem (i) = elem (i, i+k);
  2487.     }
  2488.       else if ( k < 0)
  2489.     {
  2490.       for (int i = 0; i < ndiag; i++)
  2491.         d.elem (i) = elem (i-k, i);
  2492.     }
  2493.       else
  2494.     {
  2495.       for (int i = 0; i < ndiag; i++)
  2496.         d.elem (i) = elem (i, i);
  2497.     }
  2498.     }
  2499.   else
  2500.     cerr << "diag: requested diagonal out of range\n";
  2501.  
  2502.   return d;
  2503. }
  2504.  
  2505. bool
  2506. ComplexMatrix::row_is_real_only (int i) const
  2507. {
  2508.   bool retval = true;
  2509.  
  2510.   int nc = columns ();
  2511.  
  2512.   for (int j = 0; j < nc; j++)
  2513.     {
  2514.       if (imag (elem (i, j)) != 0.0)
  2515.     {
  2516.       retval = false;
  2517.       break;
  2518.     }
  2519.     }
  2520.  
  2521.   return retval;          
  2522. }
  2523.  
  2524. bool
  2525. ComplexMatrix::column_is_real_only (int j) const
  2526. {
  2527.   bool retval = true;
  2528.  
  2529.   int nr = rows ();
  2530.  
  2531.   for (int i = 0; i < nr; i++)
  2532.     {
  2533.       if (imag (elem (i, j)) != 0.0)
  2534.     {
  2535.       retval = false;
  2536.       break;
  2537.     }
  2538.     }
  2539.  
  2540.   return retval;          
  2541. }
  2542.  
  2543. ComplexColumnVector
  2544. ComplexMatrix::row_min (void) const
  2545. {
  2546.   Array<int> index;
  2547.   return row_min (index);
  2548. }
  2549.  
  2550. ComplexColumnVector
  2551. ComplexMatrix::row_min (Array<int>& index) const
  2552. {
  2553.   ComplexColumnVector result;
  2554.  
  2555.   int nr = rows ();
  2556.   int nc = cols ();
  2557.  
  2558.   if (nr > 0 && nc > 0)
  2559.     {
  2560.       result.resize (nr);
  2561.       index.resize (nr);
  2562.  
  2563.       for (int i = 0; i < nr; i++)
  2564.         {
  2565.       int idx = 0;
  2566.  
  2567.       Complex tmp_min = elem (i, idx);
  2568.  
  2569.       bool real_only = row_is_real_only (i);
  2570.  
  2571.       double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
  2572.  
  2573.       if (xisnan (tmp_min))
  2574.         idx = -1;
  2575.       else
  2576.         {
  2577.           for (int j = 1; j < nc; j++)
  2578.         {
  2579.           Complex tmp = elem (i, j);
  2580.  
  2581.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  2582.  
  2583.           if (xisnan (tmp))
  2584.             {
  2585.               idx = -1;
  2586.               break;
  2587.             }
  2588.           else if (abs_tmp < abs_min)
  2589.             {
  2590.               idx = j;
  2591.               tmp_min = tmp;
  2592.               abs_min = abs_tmp;
  2593.             }
  2594.         }
  2595.  
  2596.           result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_min;
  2597.           index.elem (i) = idx;
  2598.         }
  2599.         }
  2600.     }
  2601.  
  2602.   return result;
  2603. }
  2604.  
  2605. ComplexColumnVector
  2606. ComplexMatrix::row_max (void) const
  2607. {
  2608.   Array<int> index;
  2609.   return row_max (index);
  2610. }
  2611.  
  2612. ComplexColumnVector
  2613. ComplexMatrix::row_max (Array<int>& index) const
  2614. {
  2615.   ComplexColumnVector result;
  2616.  
  2617.   int nr = rows ();
  2618.   int nc = cols ();
  2619.  
  2620.   if (nr > 0 && nc > 0)
  2621.     {
  2622.       result.resize (nr);
  2623.       index.resize (nr);
  2624.  
  2625.       for (int i = 0; i < nr; i++)
  2626.         {
  2627.       int idx = 0;
  2628.  
  2629.       Complex tmp_max = elem (i, idx);
  2630.  
  2631.       bool real_only = row_is_real_only (i);
  2632.  
  2633.       double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
  2634.  
  2635.       if (xisnan (tmp_max))
  2636.         idx = -1;
  2637.       else
  2638.         {
  2639.           for (int j = 1; j < nc; j++)
  2640.         {
  2641.           Complex tmp = elem (i, j);
  2642.  
  2643.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  2644.  
  2645.           if (xisnan (tmp))
  2646.             {
  2647.               idx = -1;
  2648.               break;
  2649.             }
  2650.           else if (abs_tmp > abs_max)
  2651.             {
  2652.               idx = j;
  2653.               tmp_max = tmp;
  2654.               abs_max = abs_tmp;
  2655.             }
  2656.         }
  2657.  
  2658.           result.elem (i) = (idx < 0) ? Complex_NaN_result : tmp_max;
  2659.           index.elem (i) = idx;
  2660.         }
  2661.         }
  2662.     }
  2663.  
  2664.   return result;
  2665. }
  2666.  
  2667. ComplexRowVector
  2668. ComplexMatrix::column_min (void) const
  2669. {
  2670.   Array<int> index;
  2671.   return column_min (index);
  2672. }
  2673.  
  2674. ComplexRowVector
  2675. ComplexMatrix::column_min (Array<int>& index) const
  2676. {
  2677.   ComplexRowVector result;
  2678.  
  2679.   int nr = rows ();
  2680.   int nc = cols ();
  2681.  
  2682.   if (nr > 0 && nc > 0)
  2683.     {
  2684.       result.resize (nc);
  2685.       index.resize (nc);
  2686.  
  2687.       for (int j = 0; j < nc; j++)
  2688.         {
  2689.       int idx = 0;
  2690.  
  2691.       Complex tmp_min = elem (idx, j);
  2692.  
  2693.       bool real_only = column_is_real_only (j);
  2694.  
  2695.       double abs_min = real_only ? real (tmp_min) : abs (tmp_min);
  2696.  
  2697.       if (xisnan (tmp_min))
  2698.         idx = -1;
  2699.       else
  2700.         {
  2701.           for (int i = 1; i < nr; i++)
  2702.         {
  2703.           Complex tmp = elem (i, j);
  2704.  
  2705.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  2706.  
  2707.           if (xisnan (tmp))
  2708.             {
  2709.               idx = -1;
  2710.               break;
  2711.             }
  2712.           else if (abs_tmp < abs_min)
  2713.             {
  2714.               idx = i;
  2715.               tmp_min = tmp;
  2716.               abs_min = abs_tmp;
  2717.             }
  2718.         }
  2719.  
  2720.           result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_min;
  2721.           index.elem (j) = idx;
  2722.         }
  2723.         }
  2724.     }
  2725.  
  2726.   return result;
  2727. }
  2728.  
  2729. ComplexRowVector
  2730. ComplexMatrix::column_max (void) const
  2731. {
  2732.   Array<int> index;
  2733.   return column_max (index);
  2734. }
  2735.  
  2736. ComplexRowVector
  2737. ComplexMatrix::column_max (Array<int>& index) const
  2738. {
  2739.   ComplexRowVector result;
  2740.  
  2741.   int nr = rows ();
  2742.   int nc = cols ();
  2743.  
  2744.   if (nr > 0 && nc > 0)
  2745.     {
  2746.       result.resize (nc);
  2747.       index.resize (nc);
  2748.  
  2749.       for (int j = 0; j < nc; j++)
  2750.         {
  2751.       int idx = 0;
  2752.  
  2753.       Complex tmp_max = elem (idx, j);
  2754.  
  2755.       bool real_only = column_is_real_only (j);
  2756.  
  2757.       double abs_max = real_only ? real (tmp_max) : abs (tmp_max);
  2758.  
  2759.       if (xisnan (tmp_max))
  2760.         idx = -1;
  2761.       else
  2762.         {
  2763.           for (int i = 1; i < nr; i++)
  2764.         {
  2765.           Complex tmp = elem (i, j);
  2766.  
  2767.           double abs_tmp = real_only ? real (tmp) : abs (tmp);
  2768.  
  2769.           if (xisnan (tmp))
  2770.             {
  2771.               idx = -1;
  2772.               break;
  2773.             }
  2774.           else if (abs_tmp > abs_max)
  2775.             {
  2776.               idx = i;
  2777.               tmp_max = tmp;
  2778.               abs_max = abs_tmp;
  2779.             }
  2780.         }
  2781.  
  2782.           result.elem (j) = (idx < 0) ? Complex_NaN_result : tmp_max;
  2783.           index.elem (j) = idx;
  2784.         }
  2785.         }
  2786.     }
  2787.  
  2788.   return result;
  2789. }
  2790.  
  2791. // i/o
  2792.  
  2793. ostream&
  2794. operator << (ostream& os, const ComplexMatrix& a)
  2795. {
  2796. //  int field_width = os.precision () + 7;
  2797.   for (int i = 0; i < a.rows (); i++)
  2798.     {
  2799.       for (int j = 0; j < a.cols (); j++)
  2800.     os << " " /* setw (field_width) */ << a.elem (i, j);
  2801.       os << "\n";
  2802.     }
  2803.   return os;
  2804. }
  2805.  
  2806. istream&
  2807. operator >> (istream& is, ComplexMatrix& a)
  2808. {
  2809.   int nr = a.rows ();
  2810.   int nc = a.cols ();
  2811.  
  2812.   if (nr < 1 || nc < 1)
  2813.     is.clear (ios::badbit);
  2814.   else
  2815.     {
  2816.       Complex tmp;
  2817.       for (int i = 0; i < nr; i++)
  2818.     for (int j = 0; j < nc; j++)
  2819.       {
  2820.         is >> tmp;
  2821.         if (is)
  2822.           a.elem (i, j) = tmp;
  2823.         else
  2824.           goto done;
  2825.       }
  2826.     }
  2827.  
  2828. done:
  2829.  
  2830.   return is;
  2831. }
  2832.  
  2833. ComplexMatrix
  2834. Givens (const Complex& x, const Complex& y)
  2835. {
  2836.   double cc;
  2837.   Complex cs, temp_r;
  2838.  
  2839.   F77_FCN (zlartg, ZLARTG) (x, y, cc, cs, temp_r);
  2840.  
  2841.   ComplexMatrix g (2, 2);
  2842.  
  2843.   g.elem (0, 0) = cc;
  2844.   g.elem (1, 1) = cc;
  2845.   g.elem (0, 1) = cs;
  2846.   g.elem (1, 0) = -conj (cs);
  2847.  
  2848.   return g;
  2849. }
  2850.  
  2851. ComplexMatrix
  2852. Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
  2853.        const ComplexMatrix& c)
  2854. {
  2855.   ComplexMatrix retval;
  2856.  
  2857.   // XXX FIXME XXX -- need to check that a, b, and c are all the same
  2858.   // size.
  2859.  
  2860.   // Compute Schur decompositions
  2861.  
  2862.   ComplexSCHUR as (a, "U");
  2863.   ComplexSCHUR bs (b, "U");
  2864.   
  2865.   // Transform c to new coordinates.
  2866.  
  2867.   ComplexMatrix ua = as.unitary_matrix ();
  2868.   ComplexMatrix sch_a = as.schur_matrix ();
  2869.  
  2870.   ComplexMatrix ub = bs.unitary_matrix ();
  2871.   ComplexMatrix sch_b = bs.schur_matrix ();
  2872.   
  2873.   ComplexMatrix cx = ua.hermitian () * c * ub;
  2874.  
  2875.   // Solve the sylvester equation, back-transform, and return the
  2876.   // solution.
  2877.  
  2878.   int a_nr = a.rows ();
  2879.   int b_nr = b.rows ();
  2880.  
  2881.   double scale;
  2882.   int info;
  2883.  
  2884.   Complex *pa = sch_a.fortran_vec ();
  2885.   Complex *pb = sch_b.fortran_vec ();
  2886.   Complex *px = cx.fortran_vec ();
  2887.   
  2888.   F77_XFCN (ztrsyl, ZTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb,
  2889.                  b_nr, px, a_nr, scale,
  2890.                  info, 1L, 1L));
  2891.  
  2892.   if (f77_exception_encountered)
  2893.     (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl");
  2894.   else
  2895.     {
  2896.       // XXX FIXME XXX -- check info?
  2897.  
  2898.       retval = -ua * cx * ub.hermitian ();
  2899.     }
  2900.  
  2901.   return retval;
  2902. }
  2903.  
  2904. ComplexMatrix
  2905. operator * (const ComplexMatrix& m, const Matrix& a)
  2906. {
  2907.   ComplexMatrix tmp (a);
  2908.   return m * tmp;
  2909. }
  2910.  
  2911. ComplexMatrix
  2912. operator * (const Matrix& m, const ComplexMatrix& a)
  2913. {
  2914.   ComplexMatrix tmp (m);
  2915.   return tmp * a;
  2916. }
  2917.  
  2918. ComplexMatrix
  2919. operator * (const ComplexMatrix& m, const ComplexMatrix& a)
  2920. {
  2921.   ComplexMatrix retval;
  2922.  
  2923.   int nr = m.rows ();
  2924.   int nc = m.cols ();
  2925.  
  2926.   int a_nr = a.rows ();
  2927.   int a_nc = a.cols ();
  2928.  
  2929.   if (nc != a_nr)
  2930.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  2931.   else
  2932.     {
  2933.       if (nr == 0 || nc == 0 || a_nc == 0)
  2934.     retval.resize (nr, nc, 0.0);
  2935.       else
  2936.     {
  2937.       int ld  = nr;
  2938.       int lda = a.rows ();
  2939.  
  2940.       retval.resize (nr, a_nc);
  2941.       Complex *c = retval.fortran_vec ();
  2942.  
  2943.       F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
  2944.                    m.data (), ld, a.data (), lda, 0.0,
  2945.                    c, nr, 1L, 1L));
  2946.  
  2947.       if (f77_exception_encountered)
  2948.         (*current_liboctave_error_handler)
  2949.           ("unrecoverable error in zgemm");
  2950.     }
  2951.     }
  2952.  
  2953.   return retval;
  2954. }
  2955.  
  2956. MS_CMP_OPS(ComplexMatrix, real, Complex, real)
  2957. MS_BOOL_OPS(ComplexMatrix, Complex)
  2958.  
  2959. SM_CMP_OPS(Complex, real, ComplexMatrix, real)
  2960. SM_BOOL_OPS(Complex, ComplexMatrix)
  2961.  
  2962. MM_CMP_OPS(ComplexMatrix, real, ComplexMatrix, real)
  2963. MM_BOOL_OPS(ComplexMatrix, ComplexMatrix)
  2964.  
  2965. /*
  2966. ;;; Local Variables: ***
  2967. ;;; mode: C++ ***
  2968. ;;; End: ***
  2969. */
  2970.