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