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