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 / CRowVector.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  16KB  |  752 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. int
  67. ComplexRowVector::operator == (const ComplexRowVector& 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. ComplexRowVector::operator != (const ComplexRowVector& a) const
  77. {
  78.   return !(*this == a);
  79. }
  80.  
  81. // destructive insert/delete/reorder operations
  82.  
  83. ComplexRowVector&
  84. ComplexRowVector::insert (const RowVector& a, int c)
  85. {
  86.   int a_len = a.length ();
  87.   if (c < 0 || c + 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 (c+i) = a.elem (i);
  95.  
  96.   return *this;
  97. }
  98.  
  99. ComplexRowVector&
  100. ComplexRowVector::insert (const ComplexRowVector& a, int c)
  101. {
  102.   int a_len = a.length ();
  103.   if (c < 0 || c + 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 (c+i) = a.elem (i);
  111.  
  112.   return *this;
  113. }
  114.  
  115. ComplexRowVector&
  116. ComplexRowVector::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. ComplexRowVector&
  126. ComplexRowVector::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. ComplexRowVector&
  136. ComplexRowVector::fill (double val, int c1, int c2)
  137. {
  138.   int len = length ();
  139.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  140.     {
  141.       (*current_liboctave_error_handler) ("range error for fill");
  142.       return *this;
  143.     }
  144.  
  145.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  146.  
  147.   for (int i = c1; i <= c2; i++)
  148.     elem (i) = val;
  149.  
  150.   return *this;
  151. }
  152.  
  153. ComplexRowVector&
  154. ComplexRowVector::fill (const Complex& val, int c1, int c2)
  155. {
  156.   int len = length ();
  157.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  158.     {
  159.       (*current_liboctave_error_handler) ("range error for fill");
  160.       return *this;
  161.     }
  162.  
  163.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  164.  
  165.   for (int i = c1; i <= c2; i++)
  166.     elem (i) = val;
  167.  
  168.   return *this;
  169. }
  170.  
  171. ComplexRowVector
  172. ComplexRowVector::append (const RowVector& a) const
  173. {
  174.   int len = length ();
  175.   int nc_insert = len;
  176.   ComplexRowVector retval (len + a.length ());
  177.   retval.insert (*this, 0);
  178.   retval.insert (a, nc_insert);
  179.   return retval;
  180. }
  181.  
  182. ComplexRowVector
  183. ComplexRowVector::append (const ComplexRowVector& a) const
  184. {
  185.   int len = length ();
  186.   int nc_insert = len;
  187.   ComplexRowVector retval (len + a.length ());
  188.   retval.insert (*this, 0);
  189.   retval.insert (a, nc_insert);
  190.   return retval;
  191. }
  192.  
  193. ComplexColumnVector
  194. ComplexRowVector::hermitian (void) const
  195. {
  196.   int len = length ();
  197.   return ComplexColumnVector (conj_dup (data (), len), len);
  198. }
  199.  
  200. ComplexColumnVector
  201. ComplexRowVector::transpose (void) const
  202. {
  203.   int len = length ();
  204.   return ComplexColumnVector (dup (data (), len), len);
  205. }
  206.  
  207. ComplexRowVector
  208. conj (const ComplexRowVector& a)
  209. {
  210.   int a_len = a.length ();
  211.   ComplexRowVector retval;
  212.   if (a_len > 0)
  213.     retval = ComplexRowVector (conj_dup (a.data (), a_len), a_len);
  214.   return retval;
  215. }
  216.  
  217. // resize is the destructive equivalent for this one
  218.  
  219. ComplexRowVector
  220. ComplexRowVector::extract (int c1, int c2) const
  221. {
  222.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  223.  
  224.   int new_c = c2 - c1 + 1;
  225.  
  226.   ComplexRowVector result (new_c);
  227.  
  228.   for (int i = 0; i < new_c; i++)
  229.     result.elem (i) = elem (c1+i);
  230.  
  231.   return result;
  232. }
  233.  
  234. // row vector by row vector -> row vector operations
  235.  
  236. ComplexRowVector&
  237. ComplexRowVector::operator += (const RowVector& 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. ComplexRowVector&
  257. ComplexRowVector::operator -= (const RowVector& 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. ComplexRowVector&
  277. ComplexRowVector::operator += (const ComplexRowVector& a)
  278. {
  279.   int len = length ();
  280.   if (len != a.length ())
  281.     {
  282.       (*current_liboctave_error_handler)
  283.     ("nonconformant vector += operation attempted");
  284.       return *this;
  285.     }
  286.  
  287.   if (len == 0)
  288.     return *this;
  289.  
  290.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  291.  
  292.   add2 (d, a.data (), len);
  293.   return *this;
  294. }
  295.  
  296. ComplexRowVector&
  297. ComplexRowVector::operator -= (const ComplexRowVector& a)
  298. {
  299.   int len = length ();
  300.   if (len != a.length ())
  301.     {
  302.       (*current_liboctave_error_handler)
  303.     ("nonconformant vector -= operation attempted");
  304.       return *this;
  305.     }
  306.  
  307.   if (len == 0)
  308.     return *this;
  309.  
  310.   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
  311.  
  312.   subtract2 (d, a.data (), len);
  313.   return *this;
  314. }
  315.  
  316. // row vector by scalar -> row vector operations
  317.  
  318. ComplexRowVector
  319. operator + (const ComplexRowVector& v, double s)
  320. {
  321.   int len = v.length ();
  322.   return ComplexRowVector (add (v.data (), len, s), len);
  323. }
  324.  
  325. ComplexRowVector
  326. operator - (const ComplexRowVector& v, double s)
  327. {
  328.   int len = v.length ();
  329.   return ComplexRowVector (subtract (v.data (), len, s), len);
  330. }
  331.  
  332. ComplexRowVector
  333. operator * (const ComplexRowVector& v, double s)
  334. {
  335.   int len = v.length ();
  336.   return ComplexRowVector (multiply (v.data (), len, s), len);
  337. }
  338.  
  339. ComplexRowVector
  340. operator / (const ComplexRowVector& v, double s)
  341. {
  342.   int len = v.length ();
  343.   return ComplexRowVector (divide (v.data (), len, s), len);
  344. }
  345.  
  346. ComplexRowVector
  347. operator + (const RowVector& v, const Complex& s)
  348. {
  349.   int len = v.length ();
  350.   return ComplexRowVector (add (v.data (), len, s), len);
  351. }
  352.  
  353. ComplexRowVector
  354. operator - (const RowVector& v, const Complex& s)
  355. {
  356.   int len = v.length ();
  357.   return ComplexRowVector (subtract (v.data (), len, s), len);
  358. }
  359.  
  360. ComplexRowVector
  361. operator * (const RowVector& v, const Complex& s)
  362. {
  363.   int len = v.length ();
  364.   return ComplexRowVector (multiply (v.data (), len, s), len);
  365. }
  366.  
  367. ComplexRowVector
  368. operator / (const RowVector& v, const Complex& s)
  369. {
  370.   int len = v.length ();
  371.   return ComplexRowVector (divide (v.data (), len, s), len);
  372. }
  373.  
  374. // scalar by row vector -> row vector operations
  375.  
  376. ComplexRowVector
  377. operator + (double s, const ComplexRowVector& a)
  378. {
  379.   int a_len = a.length ();
  380.   return ComplexRowVector (add (a.data (), a_len, s), a_len);
  381. }
  382.  
  383. ComplexRowVector
  384. operator - (double s, const ComplexRowVector& a)
  385. {
  386.   int a_len = a.length ();
  387.   return ComplexRowVector (subtract (s, a.data (), a_len), a_len);
  388. }
  389.  
  390. ComplexRowVector
  391. operator * (double s, const ComplexRowVector& a)
  392. {
  393.   int a_len = a.length ();
  394.   return ComplexRowVector (multiply (a.data (), a_len, s), a_len);
  395. }
  396.  
  397. ComplexRowVector
  398. operator / (double s, const ComplexRowVector& a)
  399. {
  400.   int a_len = a.length ();
  401.   return ComplexRowVector (divide (s, a.data (), a_len), a_len);
  402. }
  403.  
  404. ComplexRowVector
  405. operator + (const Complex& s, const RowVector& a)
  406. {
  407.   return ComplexRowVector ();
  408. }
  409.  
  410. ComplexRowVector
  411. operator - (const Complex& s, const RowVector& a)
  412. {
  413.   return ComplexRowVector ();
  414. }
  415.  
  416. ComplexRowVector
  417. operator * (const Complex& s, const RowVector& a)
  418. {
  419.   return ComplexRowVector ();
  420. }
  421.  
  422. ComplexRowVector
  423. operator / (const Complex& s, const RowVector& a)
  424. {
  425.   return ComplexRowVector ();
  426. }
  427.  
  428. // row vector by matrix -> row vector
  429.  
  430. ComplexRowVector
  431. operator * (const ComplexRowVector& v, const ComplexMatrix& a)
  432. {
  433.   int len = v.length ();
  434.   if (a.rows () != len)
  435.     {
  436.       (*current_liboctave_error_handler)
  437.     ("nonconformant vector multiplication attempted");
  438.       return ComplexRowVector ();
  439.     }
  440.  
  441.   if (len == 0 || a.cols () == 0)
  442.     return ComplexRowVector (0);
  443.  
  444. // Transpose A to form A'*x == (x'*A)'
  445.  
  446.   int a_nr = a.rows ();
  447.   int a_nc = a.cols ();
  448.  
  449.   char trans = 'T';
  450.   int ld = a_nr;
  451.   Complex alpha (1.0);
  452.   Complex beta (0.0);
  453.   int i_one = 1;
  454.  
  455.   Complex *y = new Complex [len];
  456.  
  457.   F77_FCN (zgemv) (&trans, &a_nc, &a_nr, &alpha, a.data (), &ld,
  458.            v.data (), &i_one, &beta, y, &i_one, 1L); 
  459.  
  460.   return ComplexRowVector (y, len);
  461. }
  462.  
  463. ComplexRowVector
  464. operator * (const RowVector& v, const ComplexMatrix& a)
  465. {
  466.   ComplexRowVector tmp (v);
  467.   return tmp * a;
  468. }
  469.  
  470. // row vector by row vector -> row vector operations
  471.  
  472. ComplexRowVector
  473. operator + (const ComplexRowVector& v, const RowVector& a)
  474. {
  475.   int len = v.length ();
  476.   if (len != a.length ())
  477.     {
  478.       (*current_liboctave_error_handler)
  479.     ("nonconformant vector addition attempted");
  480.       return ComplexRowVector ();
  481.     }
  482.  
  483.   if (len == 0)
  484.     return ComplexRowVector (0);
  485.  
  486.   return ComplexRowVector (add (v.data (), a.data (), len), len);
  487. }
  488.  
  489. ComplexRowVector
  490. operator - (const ComplexRowVector& v, const RowVector& a)
  491. {
  492.   int len = v.length ();
  493.   if (len != a.length ())
  494.     {
  495.       (*current_liboctave_error_handler)
  496.     ("nonconformant vector subtraction attempted");
  497.       return ComplexRowVector ();
  498.     }
  499.  
  500.   if (len == 0)
  501.     return ComplexRowVector (0);
  502.  
  503.   return ComplexRowVector (subtract (v.data (), a.data (), len), len);
  504. }
  505.  
  506. ComplexRowVector
  507. operator + (const RowVector& v, const ComplexRowVector& a)
  508. {
  509.   int len = v.length ();
  510.   if (len != a.length ())
  511.     {
  512.       (*current_liboctave_error_handler)
  513.     ("nonconformant vector addition attempted");
  514.       return ComplexRowVector ();
  515.     }
  516.  
  517.   if (len == 0)
  518.     return ComplexRowVector (0);
  519.  
  520.   return ComplexRowVector (add (v.data (), a.data (), len), len);
  521. }
  522.  
  523. ComplexRowVector
  524. operator - (const RowVector& v, const ComplexRowVector& a)
  525. {
  526.   int len = v.length ();
  527.   if (len != a.length ())
  528.     {
  529.       (*current_liboctave_error_handler)
  530.     ("nonconformant vector subtraction attempted");
  531.       return ComplexRowVector ();
  532.     }
  533.  
  534.   if (len == 0)
  535.     return ComplexRowVector (0);
  536.  
  537.   return ComplexRowVector (subtract (v.data (), a.data (), len), len);
  538. }
  539.  
  540. ComplexRowVector
  541. product (const ComplexRowVector& v, const RowVector& a)
  542. {
  543.   int len = v.length ();
  544.   if (len != a.length ())
  545.     {
  546.       (*current_liboctave_error_handler)
  547.     ("nonconformant vector product attempted");
  548.       return ComplexRowVector ();
  549.     }
  550.  
  551.   if (len == 0)
  552.     return ComplexRowVector (0);
  553.  
  554.   return ComplexRowVector (multiply (v.data (), a.data (), len), len);
  555. }
  556.  
  557. ComplexRowVector
  558. quotient (const ComplexRowVector& v, const RowVector& a)
  559. {
  560.   int len = v.length ();
  561.   if (len != a.length ())
  562.     {
  563.       (*current_liboctave_error_handler)
  564.     ("nonconformant vector quotient attempted");
  565.       return ComplexRowVector ();
  566.     }
  567.  
  568.   if (len == 0)
  569.     return ComplexRowVector (0);
  570.  
  571.   return ComplexRowVector (divide (v.data (), a.data (), len), len);
  572. }
  573.  
  574. ComplexRowVector
  575. product (const RowVector& v, const ComplexRowVector& a)
  576. {
  577.   int len = v.length ();
  578.   if (len != a.length ())
  579.     {
  580.       (*current_liboctave_error_handler)
  581.     ("nonconformant vector product attempted");
  582.       return ComplexRowVector ();
  583.     }
  584.  
  585.   if (len == 0)
  586.     return ComplexRowVector (0);
  587.  
  588.   return ComplexRowVector (multiply (v.data (), a.data (), len), len);
  589. }
  590.  
  591. ComplexRowVector
  592. quotient (const RowVector& v, const ComplexRowVector& a)
  593. {
  594.   int len = v.length ();
  595.   if (len != a.length ())
  596.     {
  597.       (*current_liboctave_error_handler)
  598.     ("nonconformant vector quotient attempted");
  599.       return ComplexRowVector ();
  600.     }
  601.  
  602.   if (len == 0)
  603.     return ComplexRowVector (0);
  604.  
  605.   return ComplexRowVector (divide (v.data (), a.data (), len), len);
  606. }
  607.  
  608. // other operations
  609.  
  610. ComplexRowVector
  611. map (c_c_Mapper f, const ComplexRowVector& a)
  612. {
  613.   ComplexRowVector b (a);
  614.   b.map (f);
  615.   return b;
  616. }
  617.  
  618. void
  619. ComplexRowVector::map (c_c_Mapper f)
  620. {
  621.   for (int i = 0; i < length (); i++)
  622.     elem (i) = f (elem (i));
  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. // row vector by column vector -> scalar
  699.  
  700. Complex
  701. operator * (const ComplexRowVector& v, const ColumnVector& a)
  702. {
  703.   ComplexColumnVector tmp (a);
  704.   return v * tmp;
  705. }
  706.  
  707. Complex
  708. operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
  709. {
  710.   int len = v.length ();
  711.   if (len != a.length ())
  712.     {
  713.       (*current_liboctave_error_handler)
  714.     ("nonconformant vector multiplication attempted");
  715.       return 0.0;
  716.     }
  717.  
  718.   Complex retval (0.0, 0.0);
  719.  
  720.   for (int i = 0; i < len; i++)
  721.     retval += v.elem (i) * a.elem (i);
  722.  
  723.   return retval;
  724. }
  725.  
  726. // other operations
  727.  
  728. ComplexRowVector
  729. linspace (const Complex& x1, const Complex& x2, int n)
  730. {
  731.   ComplexRowVector retval;
  732.  
  733.   if (n > 0)
  734.     {
  735.       retval.resize (n);
  736.       Complex delta = (x2 - x1) / (n - 1);
  737.       retval.elem (0) = x1;
  738.       for (int i = 1; i < n-1; i++)
  739.     retval.elem (i) = x1 + i * delta;
  740.       retval.elem (n-1) = x2;
  741.     }
  742.  
  743.   return retval;
  744. }
  745.  
  746. /*
  747. ;;; Local Variables: ***
  748. ;;; mode: C++ ***
  749. ;;; page-delimiter: "^/\\*" ***
  750. ;;; End: ***
  751. */
  752.