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 / dRowVector.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-14  |  10.0 KB  |  532 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 (dgemv) (const char*, const int*, const int*,
  42.                const double*, const double*, const int*,
  43.                const double*, const int*, const double*,
  44.                double*, const int*, long);
  45.  
  46.   double F77_FCN (ddot) (const int*, const double*, const int*,
  47.              const double*, const int*);
  48. }
  49.  
  50. /*
  51.  * Row Vector class.
  52.  */
  53.  
  54. #define KLUDGE_VECTORS
  55. #define TYPE double
  56. #define KL_VEC_TYPE RowVector
  57. #include "mx-kludge.cc"
  58. #undef KLUDGE_VECTORS
  59. #undef TYPE
  60. #undef KL_VEC_TYPE
  61.  
  62. #if 0
  63. RowVector&
  64. RowVector::resize (int n)
  65. {
  66.   if (n < 0)
  67.     {
  68.       (*current_liboctave_error_handler)
  69.     ("can't resize to negative dimension");
  70.       return *this;
  71.     }
  72.  
  73.   double *new_data = 0;
  74.   if (n > 0)
  75.     {
  76.       new_data = new double [n];
  77.       int min_len = len < n ? len : n;
  78.  
  79.       for (int i = 0; i < min_len; i++)
  80.     new_data[i] = data[i];
  81.     }
  82.  
  83.   delete [] data;
  84.   len = n;
  85.   data = new_data;
  86.  
  87.   return *this;
  88. }
  89.  
  90. RowVector&
  91. RowVector::resize (int n, double val)
  92. {
  93.   int old_len = len;
  94.   resize (n);
  95.   for (int i = old_len; i < len; i++)
  96.     data[i] = val;
  97.  
  98.   return *this;
  99. }
  100. #endif
  101.  
  102. int
  103. RowVector::operator == (const RowVector& a) const
  104. {
  105.   int len = length ();
  106.   if (len != a.length ())
  107.     return 0;
  108.   return equal (data (), a.data (), len);
  109. }
  110.  
  111. int
  112. RowVector::operator != (const RowVector& a) const
  113. {
  114.   return !(*this == a);
  115. }
  116.  
  117. RowVector&
  118. RowVector::insert (const RowVector& a, int c)
  119. {
  120.   int a_len = a.length ();
  121.   if (c < 0 || c + a_len - 1 > length ())
  122.     {
  123.       (*current_liboctave_error_handler) ("range error for insert");
  124.       return *this;
  125.     }
  126.  
  127.   for (int i = 0; i < a_len; i++)
  128.     elem (c+i) = a.elem (i);
  129.  
  130.   return *this;
  131. }
  132.  
  133. RowVector&
  134. RowVector::fill (double val)
  135. {
  136.   int len = length ();
  137.   if (len > 0)
  138.     for (int i = 0; i < len; i++)
  139.       elem (i) = val;
  140.   return *this;
  141. }
  142.  
  143. RowVector&
  144. RowVector::fill (double val, int c1, int c2)
  145. {
  146.   int len = length ();
  147.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  148.     {
  149.       (*current_liboctave_error_handler) ("range error for fill");
  150.       return *this;
  151.     }
  152.  
  153.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  154.  
  155.   for (int i = c1; i <= c2; i++)
  156.     elem (i) = val;
  157.  
  158.   return *this;
  159. }
  160.  
  161. RowVector
  162. RowVector::append (const RowVector& a) const
  163. {
  164.   int len = length ();
  165.   int nc_insert = len;
  166.   RowVector retval (len + a.length ());
  167.   retval.insert (*this, 0);
  168.   retval.insert (a, nc_insert);
  169.   return retval;
  170. }
  171.  
  172. ColumnVector
  173. RowVector::transpose (void) const
  174. {
  175.   int len = length ();
  176.   return ColumnVector (dup (data (), len), len);
  177. }
  178.  
  179. RowVector
  180. RowVector::extract (int c1, int c2) const
  181. {
  182.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  183.  
  184.   int new_c = c2 - c1 + 1;
  185.  
  186.   RowVector result (new_c);
  187.  
  188.   for (int i = 0; i < new_c; i++)
  189.     result.elem (i) = elem (c1+i);
  190.  
  191.   return result;
  192. }
  193.  
  194. // row vector by row vector -> row vector operations
  195.  
  196. RowVector&
  197. RowVector::operator += (const RowVector& a)
  198. {
  199.   int len = length ();
  200.   if (len != a.length ())
  201.     {
  202.       (*current_liboctave_error_handler)
  203.     ("nonconformant vector += operation attempted");
  204.       return *this;
  205.     }
  206.  
  207.   if (len == 0)
  208.     return *this;
  209.  
  210.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  211.  
  212.   add2 (d, a.data (), len);
  213.   return *this;
  214. }
  215.  
  216. RowVector&
  217. RowVector::operator -= (const RowVector& a)
  218. {
  219.   int len = length ();
  220.   if (len != a.length ())
  221.     {
  222.       (*current_liboctave_error_handler)
  223.     ("nonconformant vector -= operation attempted");
  224.       return *this;
  225.     }
  226.  
  227.   if (len == 0)
  228.     return *this;
  229.  
  230.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  231.  
  232.   subtract2 (d, a.data (), len);
  233.   return *this;
  234. }
  235.  
  236. // row vector by scalar -> row vector operations
  237.  
  238. ComplexRowVector
  239. operator + (const RowVector& v, const Complex& s)
  240. {
  241.   int len = v.length ();
  242.   return ComplexRowVector (add (v.data (), len, s), len);
  243. }
  244.  
  245. ComplexRowVector
  246. operator - (const RowVector& v, const Complex& s)
  247. {
  248.   int len = v.length ();
  249.   return ComplexRowVector (subtract (v.data (), len, s), len);
  250. }
  251.  
  252. ComplexRowVector
  253. operator * (const RowVector& v, const Complex& s)
  254. {
  255.   int len = v.length ();
  256.   return ComplexRowVector (multiply (v.data (), len, s), len);
  257. }
  258.  
  259. ComplexRowVector
  260. operator / (const RowVector& v, const Complex& s)
  261. {
  262.   int len = v.length ();
  263.   return ComplexRowVector (divide (v.data (), len, s), len);
  264. }
  265.  
  266. // scalar by row vector -> row vector operations
  267.  
  268. ComplexRowVector
  269. operator + (const Complex& s, const RowVector& a)
  270. {
  271.   return ComplexRowVector ();
  272. }
  273.  
  274. ComplexRowVector
  275. operator - (const Complex& s, const RowVector& a)
  276. {
  277.   return ComplexRowVector ();
  278. }
  279.  
  280. ComplexRowVector
  281. operator * (const Complex& s, const RowVector& a)
  282. {
  283.   return ComplexRowVector ();
  284. }
  285.  
  286. ComplexRowVector
  287. operator / (const Complex& s, const RowVector& a)
  288. {
  289.   return ComplexRowVector ();
  290. }
  291.  
  292. // row vector by column vector -> scalar
  293.  
  294. double
  295. operator * (const RowVector& v, const ColumnVector& a)
  296. {
  297.   int len = v.length ();
  298.   if (len != a.length ())
  299.     {
  300.       (*current_liboctave_error_handler)
  301.     ("nonconformant vector multiplication attempted");
  302.       return 0.0;
  303.     }
  304.  
  305.   int i_one = 1;
  306.   return F77_FCN (ddot) (&len, v.data (), &i_one, a.data (), &i_one);
  307. }
  308.  
  309. Complex
  310. operator * (const RowVector& v, const ComplexColumnVector& a)
  311. {
  312.   ComplexRowVector tmp (v);
  313.   return tmp * a;
  314. }
  315.  
  316. // row vector by matrix -> row vector
  317.  
  318. RowVector
  319. operator * (const RowVector& v, const Matrix& a)
  320. {
  321.   int len = v.length ();
  322.   if (a.rows () != len)
  323.     {
  324.       (*current_liboctave_error_handler)
  325.     ("nonconformant vector multiplication attempted");
  326.       return RowVector ();
  327.     }
  328.  
  329.   if (len == 0 || a.cols () == 0)
  330.     return RowVector (0);
  331.  
  332. // Transpose A to form A'*x == (x'*A)'
  333.  
  334.   int a_nr = a.rows ();
  335.   int a_nc = a.cols ();
  336.  
  337.   char trans = 'T';
  338.   int ld = a_nr;
  339.   double alpha = 1.0;
  340.   double beta  = 0.0;
  341.   int i_one = 1;
  342.  
  343.   double *y = new double [len];
  344.  
  345.   F77_FCN (dgemv) (&trans, &a_nc, &a_nr, &alpha, a.data (), &ld,
  346.            v.data (), &i_one, &beta, y, &i_one, 1L); 
  347.  
  348.   return RowVector (y, len);
  349. }
  350.  
  351. ComplexRowVector
  352. operator * (const RowVector& v, const ComplexMatrix& a)
  353. {
  354.   ComplexRowVector tmp (v);
  355.   return tmp * a;
  356. }
  357.  
  358. // row vector by row vector -> row vector operations
  359.  
  360. ComplexRowVector
  361. operator + (const RowVector& v, const ComplexRowVector& a)
  362. {
  363.   int len = v.length ();
  364.   if (len != a.length ())
  365.     {
  366.       (*current_liboctave_error_handler)
  367.     ("nonconformant vector addition attempted");
  368.       return ComplexRowVector ();
  369.     }
  370.  
  371.   if (len == 0)
  372.     return ComplexRowVector (0);
  373.  
  374.   return ComplexRowVector (add (v.data (), a.data (), len), len);
  375. }
  376.  
  377. ComplexRowVector
  378. operator - (const RowVector& v, const ComplexRowVector& a)
  379. {
  380.   int len = v.length ();
  381.   if (len != a.length ())
  382.     {
  383.       (*current_liboctave_error_handler)
  384.     ("nonconformant vector subtraction attempted");
  385.       return ComplexRowVector ();
  386.     }
  387.  
  388.   if (len == 0)
  389.     return ComplexRowVector (0);
  390.  
  391.   return ComplexRowVector (subtract (v.data (), a.data (), len), len);
  392. }
  393.  
  394. ComplexRowVector
  395. product (const RowVector& v, const ComplexRowVector& a)
  396. {
  397.   int len = v.length ();
  398.   if (len != a.length ())
  399.     {
  400.       (*current_liboctave_error_handler)
  401.     ("nonconformant vector product attempted");
  402.       return ComplexRowVector ();
  403.     }
  404.  
  405.   if (len == 0)
  406.     return ComplexRowVector (0);
  407.  
  408.   return ComplexRowVector (multiply (v.data (), a.data (), len), len);
  409. }
  410.  
  411. ComplexRowVector
  412. quotient (const RowVector& v, const ComplexRowVector& a)
  413. {
  414.   int len = v.length ();
  415.   if (len != a.length ())
  416.     {
  417.       (*current_liboctave_error_handler)
  418.     ("nonconformant vector quotient attempted");
  419.       return ComplexRowVector ();
  420.     }
  421.  
  422.   if (len == 0)
  423.     return ComplexRowVector (0);
  424.  
  425.   return ComplexRowVector (divide (v.data (), a.data (), len), len);
  426. }
  427.  
  428. // other operations
  429.  
  430. RowVector
  431. map (d_d_Mapper f, const RowVector& a)
  432. {
  433.   RowVector b (a);
  434.   b.map (f);
  435.   return b;
  436. }
  437.  
  438. void
  439. RowVector::map (d_d_Mapper f)
  440. {
  441.   for (int i = 0; i < length (); i++)
  442.     elem (i) = f (elem (i));
  443. }
  444.  
  445. RowVector
  446. linspace (double x1, double x2, int n)
  447. {
  448.   RowVector retval;
  449.  
  450.   if (n > 0)
  451.     {
  452.       retval.resize (n);
  453.       double delta = (x2 - x1) / (n - 1);
  454.       retval.elem (0) = x1;
  455.       for (int i = 1; i < n-1; i++)
  456.     retval.elem (i) = x1 + i * delta;
  457.       retval.elem (n-1) = x2;
  458.     }
  459.  
  460.   return retval;
  461. }
  462.  
  463. double
  464. RowVector::min (void) const
  465. {
  466.   int len = length ();
  467.   if (len == 0)
  468.     return 0;
  469.  
  470.   double res = elem (0);
  471.  
  472.   for (int i = 1; i < len; i++)
  473.     if (elem (i) < res)
  474.       res = elem (i);
  475.  
  476.   return res;
  477. }
  478.  
  479. double
  480. RowVector::max (void) const
  481. {
  482.   int len = length ();
  483.   if (len == 0)
  484.     return 0;
  485.  
  486.   double res = elem (0);
  487.  
  488.   for (int i = 1; i < len; i++)
  489.     if (elem (i) > res)
  490.       res = elem (i);
  491.  
  492.   return res;
  493. }
  494.  
  495. ostream&
  496. operator << (ostream& os, const RowVector& a)
  497. {
  498. //  int field_width = os.precision () + 7;
  499.   for (int i = 0; i < a.length (); i++)
  500.     os << " " /* setw (field_width) */ << a.elem (i);
  501.   return os;
  502. }
  503.  
  504. istream&
  505. operator >> (istream& is, RowVector& a)
  506. {
  507.   int len = a.length();
  508.  
  509.   if (len < 1)
  510.     is.clear (ios::badbit);
  511.   else
  512.     {
  513.       double tmp;
  514.       for (int i = 0; i < len; i++)
  515.         {
  516.           is >> tmp;
  517.           if (is)
  518.             a.elem (i) = tmp;
  519.           else
  520.             break;
  521.         }
  522.     }
  523.   return is;
  524. }
  525.  
  526. /*
  527. ;;; Local Variables: ***
  528. ;;; mode: C++ ***
  529. ;;; page-delimiter: "^/\\*" ***
  530. ;;; End: ***
  531. */
  532.