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