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 / CDiagMatrix.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  19KB  |  856 lines
  1. // DiagMatrix 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 <iostream.h>
  29.  
  30. #include <Complex.h>
  31.  
  32. #include "mx-base.h"
  33. #include "mx-inlines.cc"
  34. #include "lo-error.h"
  35.  
  36. /*
  37.  * Complex Diagonal Matrix class
  38.  */
  39.  
  40. #define KLUDGE_DIAG_MATRICES
  41. #define TYPE Complex
  42. #define KL_DMAT_TYPE ComplexDiagMatrix
  43. #include "mx-kludge.cc"
  44. #undef KLUDGE_DIAG_MATRICES
  45. #undef TYPE
  46. #undef KL_DMAT_TYPE
  47.  
  48. ComplexDiagMatrix::ComplexDiagMatrix (const RowVector& a)
  49.   : DiagArray<Complex> (a.length ())
  50. {
  51.   for (int i = 0; i < length (); i++)
  52.     elem (i, i) = a.elem (i);
  53. }
  54.  
  55. ComplexDiagMatrix::ComplexDiagMatrix (const ColumnVector& a)
  56.   : DiagArray<Complex> (a.length ())
  57. {
  58.   for (int i = 0; i < length (); i++)
  59.     elem (i, i) = a.elem (i);
  60. }
  61.  
  62. ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a)
  63.   : DiagArray<Complex> (a.rows (), a.cols ())
  64. {
  65.   for (int i = 0; i < length (); i++)
  66.     elem (i, i) = a.elem (i, i);
  67. }
  68.  
  69. int
  70. ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const
  71. {
  72.   if (rows () != a.rows () || cols () != a.cols ())
  73.     return 0;
  74.  
  75.   return equal (data (), a.data (), length ());
  76. }
  77.  
  78. int
  79. ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const
  80. {
  81.   return !(*this == a);
  82. }
  83.  
  84. ComplexDiagMatrix&
  85. ComplexDiagMatrix::fill (double val)
  86. {
  87.   for (int i = 0; i < length (); i++)
  88.     elem (i, i) = val;
  89.   return *this;
  90. }
  91.  
  92. ComplexDiagMatrix&
  93. ComplexDiagMatrix::fill (const Complex& val)
  94. {
  95.   for (int i = 0; i < length (); i++)
  96.     elem (i, i) = val;
  97.   return *this;
  98. }
  99.  
  100. ComplexDiagMatrix&
  101. ComplexDiagMatrix::fill (double val, int beg, int end)
  102. {
  103.   if (beg < 0 || end >= length () || end < beg)
  104.     {
  105.       (*current_liboctave_error_handler) ("range error for fill");
  106.       return *this;
  107.     }
  108.  
  109.   for (int i = beg; i < end; i++)
  110.     elem (i, i) = val;
  111.  
  112.   return *this;
  113. }
  114.  
  115. ComplexDiagMatrix&
  116. ComplexDiagMatrix::fill (const Complex& val, int beg, int end)
  117. {
  118.   if (beg < 0 || end >= length () || end < beg)
  119.     {
  120.       (*current_liboctave_error_handler) ("range error for fill");
  121.       return *this;
  122.     }
  123.  
  124.   for (int i = beg; i < end; i++)
  125.     elem (i, i) = val;
  126.  
  127.   return *this;
  128. }
  129.  
  130. ComplexDiagMatrix&
  131. ComplexDiagMatrix::fill (const ColumnVector& a)
  132. {
  133.   int len = length ();
  134.   if (a.length () != len)
  135.     {
  136.       (*current_liboctave_error_handler) ("range error for fill");
  137.       return *this;
  138.     }
  139.  
  140.   for (int i = 0; i < len; i++)
  141.     elem (i, i) = a.elem (i);
  142.  
  143.   return *this;
  144. }
  145.  
  146. ComplexDiagMatrix&
  147. ComplexDiagMatrix::fill (const ComplexColumnVector& a)
  148. {
  149.   int len = length ();
  150.   if (a.length () != len)
  151.     {
  152.       (*current_liboctave_error_handler) ("range error for fill");
  153.       return *this;
  154.     }
  155.  
  156.   for (int i = 0; i < len; i++)
  157.     elem (i, i) = a.elem (i);
  158.  
  159.   return *this;
  160. }
  161.  
  162. ComplexDiagMatrix&
  163. ComplexDiagMatrix::fill (const RowVector& a)
  164. {
  165.   int len = length ();
  166.   if (a.length () != len)
  167.     {
  168.       (*current_liboctave_error_handler) ("range error for fill");
  169.       return *this;
  170.     }
  171.  
  172.   for (int i = 0; i < len; i++)
  173.     elem (i, i) = a.elem (i);
  174.  
  175.   return *this;
  176. }
  177.  
  178. ComplexDiagMatrix&
  179. ComplexDiagMatrix::fill (const ComplexRowVector& a)
  180. {
  181.   int len = length ();
  182.   if (a.length () != len)
  183.     {
  184.       (*current_liboctave_error_handler) ("range error for fill");
  185.       return *this;
  186.     }
  187.  
  188.   for (int i = 0; i < len; i++)
  189.     elem (i, i) = a.elem (i);
  190.  
  191.   return *this;
  192. }
  193.  
  194. ComplexDiagMatrix&
  195. ComplexDiagMatrix::fill (const ColumnVector& a, int beg)
  196. {
  197.   int a_len = a.length ();
  198.   if (beg < 0 || beg + a_len >= length ())
  199.     {
  200.       (*current_liboctave_error_handler) ("range error for fill");
  201.       return *this;
  202.     }
  203.  
  204.   for (int i = 0; i < a_len; i++)
  205.     elem (i+beg, i+beg) = a.elem (i);
  206.  
  207.   return *this;
  208. }
  209.  
  210. ComplexDiagMatrix&
  211. ComplexDiagMatrix::fill (const ComplexColumnVector& a, int beg)
  212. {
  213.   int a_len = a.length ();
  214.   if (beg < 0 || beg + a_len >= length ())
  215.     {
  216.       (*current_liboctave_error_handler) ("range error for fill");
  217.       return *this;
  218.     }
  219.  
  220.   for (int i = 0; i < a_len; i++)
  221.     elem (i+beg, i+beg) = a.elem (i);
  222.  
  223.   return *this;
  224. }
  225.  
  226. ComplexDiagMatrix&
  227. ComplexDiagMatrix::fill (const RowVector& a, int beg)
  228. {
  229.   int a_len = a.length ();
  230.   if (beg < 0 || beg + a_len >= length ())
  231.     {
  232.       (*current_liboctave_error_handler) ("range error for fill");
  233.       return *this;
  234.     }
  235.  
  236.   for (int i = 0; i < a_len; i++)
  237.     elem (i+beg, i+beg) = a.elem (i);
  238.  
  239.   return *this;
  240. }
  241.  
  242. ComplexDiagMatrix&
  243. ComplexDiagMatrix::fill (const ComplexRowVector& a, int beg)
  244. {
  245.   int a_len = a.length ();
  246.   if (beg < 0 || beg + a_len >= length ())
  247.     {
  248.       (*current_liboctave_error_handler) ("range error for fill");
  249.       return *this;
  250.     }
  251.  
  252.   for (int i = 0; i < a_len; i++)
  253.     elem (i+beg, i+beg) = a.elem (i);
  254.  
  255.   return *this;
  256. }
  257.  
  258. ComplexDiagMatrix
  259. ComplexDiagMatrix::hermitian (void) const
  260. {
  261.   return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ());
  262. }
  263.  
  264. ComplexDiagMatrix
  265. ComplexDiagMatrix::transpose (void) const
  266. {
  267.   return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ());
  268. }
  269.  
  270. ComplexDiagMatrix
  271. conj (const ComplexDiagMatrix& a)
  272. {
  273.   ComplexDiagMatrix retval;
  274.   int a_len = a.length ();
  275.   if (a_len > 0)
  276.     retval = ComplexDiagMatrix (conj_dup (a.data (), a_len),
  277.                 a.rows (), a.cols ());
  278.   return retval;
  279. }
  280.  
  281. // resize is the destructive analog for this one
  282.  
  283. ComplexMatrix
  284. ComplexDiagMatrix::extract (int r1, int c1, int r2, int c2) const
  285. {
  286.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  287.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  288.  
  289.   int new_r = r2 - r1 + 1;
  290.   int new_c = c2 - c1 + 1;
  291.  
  292.   ComplexMatrix result (new_r, new_c);
  293.  
  294.   for (int j = 0; j < new_c; j++)
  295.     for (int i = 0; i < new_r; i++)
  296.       result.elem (i, j) = elem (r1+i, c1+j);
  297.  
  298.   return result;
  299. }
  300.  
  301. // extract row or column i.
  302.  
  303. ComplexRowVector
  304. ComplexDiagMatrix::row (int i) const
  305. {
  306.   int nr = rows ();
  307.   int nc = cols ();
  308.   if (i < 0 || i >= nr)
  309.     {
  310.       (*current_liboctave_error_handler) ("invalid row selection");
  311.       return RowVector (); 
  312.     }
  313.  
  314.   ComplexRowVector retval (nc, 0.0);
  315.   if (nr <= nc || (nr > nc && i < nc))
  316.     retval.elem (i) = elem (i, i);
  317.  
  318.   return retval;
  319. }
  320.  
  321. ComplexRowVector
  322. ComplexDiagMatrix::row (char *s) const
  323. {
  324.   if (! s)
  325.     {
  326.       (*current_liboctave_error_handler) ("invalid row selection");
  327.       return ComplexRowVector (); 
  328.     }
  329.  
  330.   char c = *s;
  331.   if (c == 'f' || c == 'F')
  332.     return row (0);
  333.   else if (c == 'l' || c == 'L')
  334.     return row (rows () - 1);
  335.   else
  336.     {
  337.       (*current_liboctave_error_handler) ("invalid row selection");
  338.       return ComplexRowVector ();
  339.     }
  340. }
  341.  
  342. ComplexColumnVector
  343. ComplexDiagMatrix::column (int i) const
  344. {
  345.   int nr = rows ();
  346.   int nc = cols ();
  347.   if (i < 0 || i >= nc)
  348.     {
  349.       (*current_liboctave_error_handler) ("invalid column selection");
  350.       return ColumnVector (); 
  351.     }
  352.  
  353.   ComplexColumnVector retval (nr, 0.0);
  354.   if (nr >= nc || (nr < nc && i < nr))
  355.     retval.elem (i) = elem (i, i);
  356.  
  357.   return retval;
  358. }
  359.  
  360. ComplexColumnVector
  361. ComplexDiagMatrix::column (char *s) const
  362. {
  363.   if (! s)
  364.     {
  365.       (*current_liboctave_error_handler) ("invalid column selection");
  366.       return ColumnVector (); 
  367.     }
  368.  
  369.   char c = *s;
  370.   if (c == 'f' || c == 'F')
  371.     return column (0);
  372.   else if (c == 'l' || c == 'L')
  373.     return column (cols () - 1);
  374.   else
  375.     {
  376.       (*current_liboctave_error_handler) ("invalid column selection");
  377.       return ColumnVector (); 
  378.     }
  379. }
  380.  
  381. ComplexDiagMatrix
  382. ComplexDiagMatrix::inverse (void) const
  383. {
  384.   int info;
  385.   return inverse (info);
  386. }
  387.  
  388. ComplexDiagMatrix
  389. ComplexDiagMatrix::inverse (int& info) const
  390. {
  391.   int nr = rows ();
  392.   int nc = cols ();
  393.   if (nr != nc)
  394.     {
  395.       (*current_liboctave_error_handler) ("inverse requires square matrix");
  396.       return DiagMatrix ();
  397.     }
  398.  
  399.   ComplexDiagMatrix retval (nr, nc);
  400.  
  401.   info = 0;
  402.   for (int i = 0; i < length (); i++)
  403.     {
  404.       if (elem (i, i) == 0.0)
  405.     {
  406.       info = -1;
  407.       return *this;
  408.     }
  409.       else
  410.     retval.elem (i, i) = 1.0 / elem (i, i);
  411.     }
  412.  
  413.   return *this;
  414. }
  415.  
  416. // diagonal matrix by diagonal matrix -> diagonal matrix operations
  417.  
  418. ComplexDiagMatrix&
  419. ComplexDiagMatrix::operator += (const DiagMatrix& a)
  420. {
  421.   int nr = rows ();
  422.   int nc = cols ();
  423.   if (nr != a.rows () || nc != a.cols ())
  424.     {
  425.       (*current_liboctave_error_handler)
  426.     ("nonconformant matrix += operation attempted");
  427.       return *this;
  428.     }
  429.  
  430.   if (nr == 0 || nc == 0)
  431.     return *this;
  432.  
  433.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  434.  
  435.   add2 (d, a.data (), length ());
  436.   return *this;
  437. }
  438.  
  439. ComplexDiagMatrix&
  440. ComplexDiagMatrix::operator -= (const DiagMatrix& a)
  441. {
  442.   int nr = rows ();
  443.   int nc = cols ();
  444.   if (nr != a.rows () || nc != a.cols ())
  445.     {
  446.       (*current_liboctave_error_handler)
  447.     ("nonconformant matrix -= operation attempted");
  448.       return *this;
  449.     }
  450.  
  451.   if (nr == 0 || nc == 0)
  452.     return *this;
  453.  
  454.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  455.  
  456.   subtract2 (d, a.data (), length ());
  457.   return *this;
  458. }
  459.  
  460. ComplexDiagMatrix&
  461. ComplexDiagMatrix::operator += (const ComplexDiagMatrix& a)
  462. {
  463.   int nr = rows ();
  464.   int nc = cols ();
  465.   if (nr != a.rows () || nc != a.cols ())
  466.     {
  467.       (*current_liboctave_error_handler)
  468.     ("nonconformant matrix += operation attempted");
  469.       return *this;
  470.     }
  471.  
  472.   if (nr == 0 || nc == 0)
  473.     return *this;
  474.  
  475.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  476.  
  477.   add2 (d, a.data (), length ());
  478.   return *this;
  479. }
  480.  
  481. ComplexDiagMatrix&
  482. ComplexDiagMatrix::operator -= (const ComplexDiagMatrix& a)
  483. {
  484.   int nr = rows ();
  485.   int nc = cols ();
  486.   if (nr != a.rows () || nc != a.cols ())
  487.     {
  488.       (*current_liboctave_error_handler)
  489.     ("nonconformant matrix -= operation attempted");
  490.       return *this;
  491.     }
  492.  
  493.   if (nr == 0 || nc == 0)
  494.     return *this;
  495.  
  496.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  497.  
  498.   subtract2 (d, a.data (), length ());
  499.   return *this;
  500. }
  501.  
  502. // diagonal matrix by scalar -> diagonal matrix operations
  503.  
  504. ComplexDiagMatrix
  505. operator * (const ComplexDiagMatrix& a, double s)
  506. {
  507.   return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
  508.                 a.rows (), a.cols ());
  509. }
  510.  
  511. ComplexDiagMatrix
  512. operator / (const ComplexDiagMatrix& a, double s)
  513. {
  514.   return ComplexDiagMatrix (divide (a.data (), a.length (), s),
  515.                 a.rows (), a.cols ());
  516. }
  517.  
  518. ComplexDiagMatrix
  519. operator * (const DiagMatrix& a, const Complex& s)
  520. {
  521.   return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
  522.                 a.rows (), a.cols ());
  523. }
  524.  
  525. ComplexDiagMatrix
  526. operator / (const DiagMatrix& a, const Complex& s)
  527. {
  528.   return ComplexDiagMatrix (divide (a.data (), a.length (), s),
  529.                 a.rows (), a.cols ());
  530. }
  531.  
  532. // scalar by diagonal matrix -> diagonal matrix operations
  533.  
  534. ComplexDiagMatrix
  535. operator * (double s, const ComplexDiagMatrix& a)
  536. {
  537.   return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
  538.                 a.rows (), a.cols ());
  539. }
  540.  
  541. ComplexDiagMatrix
  542. operator * (const Complex& s, const DiagMatrix& a)
  543. {
  544.   return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
  545.                 a.rows (), a.cols ());
  546. }
  547.  
  548. // diagonal matrix by diagonal matrix -> diagonal matrix operations
  549.  
  550. ComplexDiagMatrix
  551. operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
  552. {
  553.   int nr_a = a.rows ();
  554.   int nc_a = a.cols ();
  555.   int nr_b = b.rows ();
  556.   int nc_b = b.cols ();
  557.   if (nc_a != nr_b)
  558.     {
  559.       (*current_liboctave_error_handler)
  560.         ("nonconformant matrix multiplication attempted");
  561.       return ComplexDiagMatrix ();
  562.     }
  563.  
  564.   if (nr_a == 0 || nc_a == 0 || nc_b == 0)
  565.     return ComplexDiagMatrix (nr_a, nc_a, 0.0);
  566.  
  567.   ComplexDiagMatrix c (nr_a, nc_b);
  568.  
  569.   int len = nr_a < nc_b ? nr_a : nc_b;
  570.  
  571.   for (int i = 0; i < len; i++)
  572.     {
  573.       Complex a_element = a.elem (i, i);
  574.       Complex b_element = b.elem (i, i);
  575.  
  576.       if (a_element == 0.0 || b_element == 0.0)
  577.         c.elem (i, i) = 0.0;
  578.       else if (a_element == 1.0)
  579.         c.elem (i, i) = b_element;
  580.       else if (b_element == 1.0)
  581.         c.elem (i, i) = a_element;
  582.       else
  583.         c.elem (i, i) = a_element * b_element;
  584.     }
  585.  
  586.   return c;
  587. }
  588.  
  589. ComplexDiagMatrix
  590. operator + (const ComplexDiagMatrix& m, const DiagMatrix& a)
  591. {
  592.   int nr = m.rows ();
  593.   int nc = m.cols ();
  594.   if (nr != a.rows () || nc != a.cols ())
  595.     {
  596.       (*current_liboctave_error_handler)
  597.     ("nonconformant matrix addition attempted");
  598.       return ComplexDiagMatrix ();
  599.     }
  600.  
  601.   if (nr == 0 || nc == 0)
  602.     return ComplexDiagMatrix (nr, nc);
  603.  
  604.   return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
  605. }
  606.  
  607. ComplexDiagMatrix
  608. operator - (const ComplexDiagMatrix& m, const DiagMatrix& a)
  609. {
  610.   int nr = m.rows ();
  611.   int nc = m.cols ();
  612.   if (nr != a.rows () || nc != a.cols ())
  613.     {
  614.       (*current_liboctave_error_handler)
  615.     ("nonconformant matrix subtraction attempted");
  616.       return ComplexDiagMatrix ();
  617.     }
  618.  
  619.   if (nr == 0 || nc == 0)
  620.     return ComplexDiagMatrix (nr, nc);
  621.  
  622.   return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
  623.                 nr, nc);
  624. }
  625.  
  626. ComplexDiagMatrix
  627. operator * (const ComplexDiagMatrix& a, const DiagMatrix& b)
  628. {
  629.   int nr_a = a.rows ();
  630.   int nc_a = a.cols ();
  631.   int nr_b = b.rows ();
  632.   int nc_b = b.cols ();
  633.   if (nc_a != nr_b)
  634.     {
  635.       (*current_liboctave_error_handler)
  636.         ("nonconformant matrix multiplication attempted");
  637.       return ComplexDiagMatrix ();
  638.     }
  639.  
  640.   if (nr_a == 0 || nc_a == 0 || nc_b == 0)
  641.     return ComplexDiagMatrix (nr_a, nc_a, 0.0);
  642.  
  643.   ComplexDiagMatrix c (nr_a, nc_b);
  644.  
  645.   int len = nr_a < nc_b ? nr_a : nc_b;
  646.  
  647.   for (int i = 0; i < len; i++)
  648.     {
  649.       Complex a_element = a.elem (i, i);
  650.       double b_element = b.elem (i, i);
  651.  
  652.       if (a_element == 0.0 || b_element == 0.0)
  653.         c.elem (i, i) = 0.0;
  654.       else if (a_element == 1.0)
  655.         c.elem (i, i) = b_element;
  656.       else if (b_element == 1.0)
  657.         c.elem (i, i) = a_element;
  658.       else
  659.         c.elem (i, i) = a_element * b_element;
  660.     }
  661.  
  662.   return c;
  663. }
  664.  
  665. ComplexDiagMatrix
  666. operator + (const DiagMatrix& m, const ComplexDiagMatrix& a)
  667. {
  668.   int nr = m.rows ();
  669.   int nc = m.cols ();
  670.   if (nr != a.rows () || nc != a.cols ())
  671.     {
  672.       (*current_liboctave_error_handler)
  673.     ("nonconformant matrix addition attempted");
  674.       return ComplexDiagMatrix ();
  675.     }
  676.  
  677.   if (nc == 0 || nr == 0)
  678.     return ComplexDiagMatrix (nr, nc);
  679.  
  680.   return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()),  nr, nc);
  681. }
  682.  
  683. ComplexDiagMatrix
  684. operator - (const DiagMatrix& m, const ComplexDiagMatrix& a)
  685. {
  686.   int nr = m.rows ();
  687.   int nc = m.cols ();
  688.   if (nr != a.rows () || nc != a.cols ())
  689.     {
  690.       (*current_liboctave_error_handler)
  691.     ("nonconformant matrix subtraction attempted");
  692.       return ComplexDiagMatrix ();
  693.     }
  694.  
  695.   if (nc == 0 || nr == 0)
  696.     return ComplexDiagMatrix (nr, nc);
  697.  
  698.   return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
  699.                 nr, nc);
  700. }
  701.  
  702. ComplexDiagMatrix
  703. operator * (const DiagMatrix& a, const ComplexDiagMatrix& b)
  704. {
  705.   int nr_a = a.rows ();
  706.   int nc_a = a.cols ();
  707.   int nr_b = b.rows ();
  708.   int nc_b = b.cols ();
  709.   if (nc_a != nr_b)
  710.     {
  711.       (*current_liboctave_error_handler)
  712.         ("nonconformant matrix multiplication attempted");
  713.       return ComplexDiagMatrix ();
  714.     }
  715.  
  716.   if (nr_a == 0 || nc_a == 0 || nc_b == 0)
  717.     return ComplexDiagMatrix (nr_a, nc_a, 0.0);
  718.  
  719.   ComplexDiagMatrix c (nr_a, nc_b);
  720.  
  721.   int len = nr_a < nc_b ? nr_a : nc_b;
  722.  
  723.   for (int i = 0; i < len; i++)
  724.     {
  725.       double a_element = a.elem (i, i);
  726.       Complex b_element = b.elem (i, i);
  727.  
  728.       if (a_element == 0.0 || b_element == 0.0)
  729.         c.elem (i, i) = 0.0;
  730.       else if (a_element == 1.0)
  731.         c.elem (i, i) = b_element;
  732.       else if (b_element == 1.0)
  733.         c.elem (i, i) = a_element;
  734.       else
  735.         c.elem (i, i) = a_element * b_element;
  736.     }
  737.  
  738.   return c;
  739. }
  740.  
  741. ComplexDiagMatrix
  742. product (const ComplexDiagMatrix& m, const DiagMatrix& a)
  743. {
  744.   int nr = m.rows ();
  745.   int nc = m.cols ();
  746.   if (nr != a.rows () || nc != a.cols ())
  747.     {
  748.       (*current_liboctave_error_handler)
  749.     ("nonconformant matrix product attempted");
  750.       return ComplexDiagMatrix ();
  751.     }
  752.  
  753.   if (nr == 0 || nc == 0)
  754.     return ComplexDiagMatrix (nr, nc);
  755.  
  756.   return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
  757.                 nr, nc);
  758. }
  759.  
  760. ComplexDiagMatrix
  761. product (const DiagMatrix& m, const ComplexDiagMatrix& a)
  762. {
  763.   int nr = m.rows ();
  764.   int nc = m.cols ();
  765.   if (nr != a.rows () || nc != a.cols ())
  766.     {
  767.       (*current_liboctave_error_handler)
  768.     ("nonconformant matrix product attempted");
  769.       return ComplexDiagMatrix ();
  770.     }
  771.  
  772.   if (nc == 0 || nr == 0)
  773.     return ComplexDiagMatrix (nr, nc);
  774.  
  775.   return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
  776.                 nr, nc);
  777. }
  778.  
  779. // other operations
  780.  
  781. ComplexColumnVector
  782. ComplexDiagMatrix::diag (void) const
  783. {
  784.   return diag (0);
  785. }
  786.  
  787. // Could be optimized...
  788.  
  789. ComplexColumnVector
  790. ComplexDiagMatrix::diag (int k) const
  791. {
  792.   int nnr = rows ();
  793.   int nnc = cols ();
  794.   if (k > 0)
  795.     nnc -= k;
  796.   else if (k < 0)
  797.     nnr += k;
  798.  
  799.   ComplexColumnVector d;
  800.  
  801.   if (nnr > 0 && nnc > 0)
  802.     {
  803.       int ndiag = (nnr < nnc) ? nnr : nnc;
  804.  
  805.       d.resize (ndiag);
  806.  
  807.       if (k > 0)
  808.     {
  809.       for (int i = 0; i < ndiag; i++)
  810.         d.elem (i) = elem (i, i+k);
  811.     }
  812.       else if ( k < 0)
  813.     {
  814.       for (int i = 0; i < ndiag; i++)
  815.         d.elem (i) = elem (i-k, i);
  816.     }
  817.       else
  818.     {
  819.       for (int i = 0; i < ndiag; i++)
  820.         d.elem (i) = elem (i, i);
  821.     }
  822.     }
  823.   else
  824.     cerr << "diag: requested diagonal out of range\n";
  825.  
  826.   return d;
  827. }
  828.  
  829. // i/o
  830.  
  831. ostream&
  832. operator << (ostream& os, const ComplexDiagMatrix& a)
  833. {
  834.   Complex ZERO (0.0);
  835. //  int field_width = os.precision () + 7;
  836.   for (int i = 0; i < a.rows (); i++)
  837.     {
  838.       for (int j = 0; j < a.cols (); j++)
  839.     {
  840.       if (i == j)
  841.         os << " " /* setw (field_width) */ << a.elem (i, i);
  842.       else
  843.         os << " " /* setw (field_width) */ << ZERO;
  844.     }
  845.       os << "\n";
  846.     }
  847.   return os;
  848. }
  849.  
  850. /*
  851. ;;; Local Variables: ***
  852. ;;; mode: C++ ***
  853. ;;; page-delimiter: "^/\\*" ***
  854. ;;; End: ***
  855. */
  856.