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