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 / dRowVector.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  7KB  |  391 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. int
  63. RowVector::operator == (const RowVector& a) const
  64. {
  65.   int len = length ();
  66.   if (len != a.length ())
  67.     return 0;
  68.   return equal (data (), a.data (), len);
  69. }
  70.  
  71. int
  72. RowVector::operator != (const RowVector& a) const
  73. {
  74.   return !(*this == a);
  75. }
  76.  
  77. RowVector&
  78. RowVector::insert (const RowVector& a, int c)
  79. {
  80.   int a_len = a.length ();
  81.   if (c < 0 || c + a_len - 1 > length ())
  82.     {
  83.       (*current_liboctave_error_handler) ("range error for insert");
  84.       return *this;
  85.     }
  86.  
  87.   for (int i = 0; i < a_len; i++)
  88.     elem (c+i) = a.elem (i);
  89.  
  90.   return *this;
  91. }
  92.  
  93. RowVector&
  94. RowVector::fill (double val)
  95. {
  96.   int len = length ();
  97.   if (len > 0)
  98.     for (int i = 0; i < len; i++)
  99.       elem (i) = val;
  100.   return *this;
  101. }
  102.  
  103. RowVector&
  104. RowVector::fill (double val, int c1, int c2)
  105. {
  106.   int len = length ();
  107.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  108.     {
  109.       (*current_liboctave_error_handler) ("range error for fill");
  110.       return *this;
  111.     }
  112.  
  113.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  114.  
  115.   for (int i = c1; i <= c2; i++)
  116.     elem (i) = val;
  117.  
  118.   return *this;
  119. }
  120.  
  121. RowVector
  122. RowVector::append (const RowVector& a) const
  123. {
  124.   int len = length ();
  125.   int nc_insert = len;
  126.   RowVector retval (len + a.length ());
  127.   retval.insert (*this, 0);
  128.   retval.insert (a, nc_insert);
  129.   return retval;
  130. }
  131.  
  132. ColumnVector
  133. RowVector::transpose (void) const
  134. {
  135.   int len = length ();
  136.   return ColumnVector (dup (data (), len), len);
  137. }
  138.  
  139. RowVector
  140. real (const ComplexRowVector& a)
  141. {
  142.   int a_len = a.length ();
  143.   RowVector retval;
  144.   if (a_len > 0)
  145.     retval = RowVector (real_dup (a.data (), a_len), a_len);
  146.   return retval;
  147. }
  148.  
  149. RowVector
  150. imag (const ComplexRowVector& a)
  151. {
  152.   int a_len = a.length ();
  153.   RowVector retval;
  154.   if (a_len > 0)
  155.     retval = RowVector (imag_dup (a.data (), a_len), a_len);
  156.   return retval;
  157. }
  158.  
  159. RowVector
  160. RowVector::extract (int c1, int c2) const
  161. {
  162.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  163.  
  164.   int new_c = c2 - c1 + 1;
  165.  
  166.   RowVector result (new_c);
  167.  
  168.   for (int i = 0; i < new_c; i++)
  169.     result.elem (i) = elem (c1+i);
  170.  
  171.   return result;
  172. }
  173.  
  174. // row vector by row vector -> row vector operations
  175.  
  176. RowVector&
  177. RowVector::operator += (const RowVector& a)
  178. {
  179.   int len = length ();
  180.   if (len != a.length ())
  181.     {
  182.       (*current_liboctave_error_handler)
  183.     ("nonconformant vector += operation attempted");
  184.       return *this;
  185.     }
  186.  
  187.   if (len == 0)
  188.     return *this;
  189.  
  190.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  191.  
  192.   add2 (d, a.data (), len);
  193.   return *this;
  194. }
  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.   subtract2 (d, a.data (), len);
  213.   return *this;
  214. }
  215.  
  216. // row vector by matrix -> row vector
  217.  
  218. RowVector
  219. operator * (const RowVector& v, const Matrix& a)
  220. {
  221.   int len = v.length ();
  222.   if (a.rows () != len)
  223.     {
  224.       (*current_liboctave_error_handler)
  225.     ("nonconformant vector multiplication attempted");
  226.       return RowVector ();
  227.     }
  228.  
  229.   if (len == 0 || a.cols () == 0)
  230.     return RowVector (0);
  231.  
  232. // Transpose A to form A'*x == (x'*A)'
  233.  
  234.   int a_nr = a.rows ();
  235.   int a_nc = a.cols ();
  236.  
  237.   char trans = 'T';
  238.   int ld = a_nr;
  239.   double alpha = 1.0;
  240.   double beta  = 0.0;
  241.   int i_one = 1;
  242.  
  243.   double *y = new double [len];
  244.  
  245.   F77_FCN (dgemv) (&trans, &a_nc, &a_nr, &alpha, a.data (), &ld,
  246.            v.data (), &i_one, &beta, y, &i_one, 1L); 
  247.  
  248.   return RowVector (y, len);
  249. }
  250.  
  251. // other operations
  252.  
  253. RowVector
  254. map (d_d_Mapper f, const RowVector& a)
  255. {
  256.   RowVector b (a);
  257.   b.map (f);
  258.   return b;
  259. }
  260.  
  261. RowVector
  262. map (d_c_Mapper f, const ComplexRowVector& a)
  263. {
  264.   int a_len = a.length ();
  265.   RowVector b (a_len);
  266.   for (int i = 0; i < a_len; i++)
  267.     b.elem (i) = f (a.elem (i));
  268.   return b;
  269. }
  270.  
  271. void
  272. RowVector::map (d_d_Mapper f)
  273. {
  274.   for (int i = 0; i < length (); i++)
  275.     elem (i) = f (elem (i));
  276. }
  277.  
  278. double
  279. RowVector::min (void) const
  280. {
  281.   int len = length ();
  282.   if (len == 0)
  283.     return 0;
  284.  
  285.   double res = elem (0);
  286.  
  287.   for (int i = 1; i < len; i++)
  288.     if (elem (i) < res)
  289.       res = elem (i);
  290.  
  291.   return res;
  292. }
  293.  
  294. double
  295. RowVector::max (void) const
  296. {
  297.   int len = length ();
  298.   if (len == 0)
  299.     return 0;
  300.  
  301.   double res = elem (0);
  302.  
  303.   for (int i = 1; i < len; i++)
  304.     if (elem (i) > res)
  305.       res = elem (i);
  306.  
  307.   return res;
  308. }
  309.  
  310. ostream&
  311. operator << (ostream& os, const RowVector& a)
  312. {
  313. //  int field_width = os.precision () + 7;
  314.   for (int i = 0; i < a.length (); i++)
  315.     os << " " /* setw (field_width) */ << a.elem (i);
  316.   return os;
  317. }
  318.  
  319. istream&
  320. operator >> (istream& is, RowVector& a)
  321. {
  322.   int len = a.length();
  323.  
  324.   if (len < 1)
  325.     is.clear (ios::badbit);
  326.   else
  327.     {
  328.       double tmp;
  329.       for (int i = 0; i < len; i++)
  330.         {
  331.           is >> tmp;
  332.           if (is)
  333.             a.elem (i) = tmp;
  334.           else
  335.             break;
  336.         }
  337.     }
  338.   return is;
  339. }
  340.  
  341. // other operations
  342.  
  343. RowVector
  344. linspace (double x1, double x2, int n)
  345. {
  346.   RowVector retval;
  347.  
  348.   if (n > 0)
  349.     {
  350.       retval.resize (n);
  351.       double delta = (x2 - x1) / (n - 1);
  352.       retval.elem (0) = x1;
  353.       for (int i = 1; i < n-1; i++)
  354.     retval.elem (i) = x1 + i * delta;
  355.       retval.elem (n-1) = x2;
  356.     }
  357.  
  358.   return retval;
  359. }
  360.  
  361. // row vector by column vector -> scalar
  362.  
  363. double
  364. operator * (const RowVector& v, const ColumnVector& a)
  365. {
  366.   int len = v.length ();
  367.   if (len != a.length ())
  368.     {
  369.       (*current_liboctave_error_handler)
  370.     ("nonconformant vector multiplication attempted");
  371.       return 0.0;
  372.     }
  373.  
  374.   int i_one = 1;
  375.   return F77_FCN (ddot) (&len, v.data (), &i_one, a.data (), &i_one);
  376. }
  377.  
  378. Complex
  379. operator * (const RowVector& v, const ComplexColumnVector& a)
  380. {
  381.   ComplexRowVector tmp (v);
  382.   return (tmp * a);
  383. }
  384.  
  385. /*
  386. ;;; Local Variables: ***
  387. ;;; mode: C++ ***
  388. ;;; page-delimiter: "^/\\*" ***
  389. ;;; End: ***
  390. */
  391.