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