CRowVector.cc

Go to the documentation of this file.
00001 // RowVector manipulations.
00002 /*
00003 
00004 Copyright (C) 1994-2012 John W. Eaton
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 3 of the License, or (at your
00011 option) any later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, see
00020 <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #ifdef HAVE_CONFIG_H
00025 #include <config.h>
00026 #endif
00027 
00028 #include <iostream>
00029 
00030 #include "Array-util.h"
00031 #include "f77-fcn.h"
00032 #include "functor.h"
00033 #include "lo-error.h"
00034 #include "mx-base.h"
00035 #include "mx-inlines.cc"
00036 #include "oct-cmplx.h"
00037 
00038 // Fortran functions we call.
00039 
00040 extern "C"
00041 {
00042   F77_RET_T
00043   F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
00044                            const octave_idx_type&, const octave_idx_type&,
00045                            const Complex&, const Complex*,
00046                            const octave_idx_type&, const Complex*,
00047                            const octave_idx_type&, const Complex&, Complex*,
00048                            const octave_idx_type&
00049                            F77_CHAR_ARG_LEN_DECL);
00050 
00051   F77_RET_T
00052   F77_FUNC (xzdotu, XZDOTU) (const octave_idx_type&, const Complex*,
00053                              const octave_idx_type&, const Complex*,
00054                              const octave_idx_type&, Complex&);
00055 }
00056 
00057 // Complex Row Vector class
00058 
00059 bool
00060 ComplexRowVector::operator == (const ComplexRowVector& a) const
00061 {
00062   octave_idx_type len = length ();
00063   if (len != a.length ())
00064     return 0;
00065   return mx_inline_equal (len, data (), a.data ());
00066 }
00067 
00068 bool
00069 ComplexRowVector::operator != (const ComplexRowVector& a) const
00070 {
00071   return !(*this == a);
00072 }
00073 
00074 // destructive insert/delete/reorder operations
00075 
00076 ComplexRowVector&
00077 ComplexRowVector::insert (const RowVector& a, octave_idx_type c)
00078 {
00079   octave_idx_type a_len = a.length ();
00080 
00081   if (c < 0 || c + a_len > length ())
00082     {
00083       (*current_liboctave_error_handler) ("range error for insert");
00084       return *this;
00085     }
00086 
00087   if (a_len > 0)
00088     {
00089       make_unique ();
00090 
00091       for (octave_idx_type i = 0; i < a_len; i++)
00092         xelem (c+i) = a.elem (i);
00093     }
00094 
00095   return *this;
00096 }
00097 
00098 ComplexRowVector&
00099 ComplexRowVector::insert (const ComplexRowVector& a, octave_idx_type c)
00100 {
00101   octave_idx_type a_len = a.length ();
00102 
00103   if (c < 0 || c + a_len > length ())
00104     {
00105       (*current_liboctave_error_handler) ("range error for insert");
00106       return *this;
00107     }
00108 
00109   if (a_len > 0)
00110     {
00111       make_unique ();
00112 
00113       for (octave_idx_type i = 0; i < a_len; i++)
00114         xelem (c+i) = a.elem (i);
00115     }
00116 
00117   return *this;
00118 }
00119 
00120 ComplexRowVector&
00121 ComplexRowVector::fill (double val)
00122 {
00123   octave_idx_type len = length ();
00124 
00125   if (len > 0)
00126     {
00127       make_unique ();
00128 
00129       for (octave_idx_type i = 0; i < len; i++)
00130         xelem (i) = val;
00131     }
00132 
00133   return *this;
00134 }
00135 
00136 ComplexRowVector&
00137 ComplexRowVector::fill (const Complex& val)
00138 {
00139   octave_idx_type len = length ();
00140 
00141   if (len > 0)
00142     {
00143       make_unique ();
00144 
00145       for (octave_idx_type i = 0; i < len; i++)
00146         xelem (i) = val;
00147     }
00148 
00149   return *this;
00150 }
00151 
00152 ComplexRowVector&
00153 ComplexRowVector::fill (double val, octave_idx_type c1, octave_idx_type c2)
00154 {
00155   octave_idx_type len = length ();
00156 
00157   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
00158     {
00159       (*current_liboctave_error_handler) ("range error for fill");
00160       return *this;
00161     }
00162 
00163   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00164 
00165   if (c2 >= c1)
00166     {
00167       make_unique ();
00168 
00169       for (octave_idx_type i = c1; i <= c2; i++)
00170         xelem (i) = val;
00171     }
00172 
00173   return *this;
00174 }
00175 
00176 ComplexRowVector&
00177 ComplexRowVector::fill (const Complex& val, octave_idx_type c1, octave_idx_type c2)
00178 {
00179   octave_idx_type len = length ();
00180 
00181   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
00182     {
00183       (*current_liboctave_error_handler) ("range error for fill");
00184       return *this;
00185     }
00186 
00187   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00188 
00189   if (c2 >= c1)
00190     {
00191       make_unique ();
00192 
00193       for (octave_idx_type i = c1; i <= c2; i++)
00194         xelem (i) = val;
00195     }
00196 
00197   return *this;
00198 }
00199 
00200 ComplexRowVector
00201 ComplexRowVector::append (const RowVector& a) const
00202 {
00203   octave_idx_type len = length ();
00204   octave_idx_type nc_insert = len;
00205   ComplexRowVector retval (len + a.length ());
00206   retval.insert (*this, 0);
00207   retval.insert (a, nc_insert);
00208   return retval;
00209 }
00210 
00211 ComplexRowVector
00212 ComplexRowVector::append (const ComplexRowVector& a) const
00213 {
00214   octave_idx_type len = length ();
00215   octave_idx_type nc_insert = len;
00216   ComplexRowVector retval (len + a.length ());
00217   retval.insert (*this, 0);
00218   retval.insert (a, nc_insert);
00219   return retval;
00220 }
00221 
00222 ComplexColumnVector
00223 ComplexRowVector::hermitian (void) const
00224 {
00225   return MArray<Complex>::hermitian (std::conj);
00226 }
00227 
00228 ComplexColumnVector
00229 ComplexRowVector::transpose (void) const
00230 {
00231   return MArray<Complex>::transpose ();
00232 }
00233 
00234 ComplexRowVector
00235 conj (const ComplexRowVector& a)
00236 {
00237   return do_mx_unary_map<Complex, Complex, std::conj<double> > (a);
00238 }
00239 
00240 // resize is the destructive equivalent for this one
00241 
00242 ComplexRowVector
00243 ComplexRowVector::extract (octave_idx_type c1, octave_idx_type c2) const
00244 {
00245   if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
00246 
00247   octave_idx_type new_c = c2 - c1 + 1;
00248 
00249   ComplexRowVector result (new_c);
00250 
00251   for (octave_idx_type i = 0; i < new_c; i++)
00252     result.elem (i) = elem (c1+i);
00253 
00254   return result;
00255 }
00256 
00257 ComplexRowVector
00258 ComplexRowVector::extract_n (octave_idx_type r1, octave_idx_type n) const
00259 {
00260   ComplexRowVector result (n);
00261 
00262   for (octave_idx_type i = 0; i < n; i++)
00263     result.elem (i) = elem (r1+i);
00264 
00265   return result;
00266 }
00267 
00268 // row vector by row vector -> row vector operations
00269 
00270 ComplexRowVector&
00271 ComplexRowVector::operator += (const RowVector& a)
00272 {
00273   octave_idx_type len = length ();
00274 
00275   octave_idx_type a_len = a.length ();
00276 
00277   if (len != a_len)
00278     {
00279       gripe_nonconformant ("operator +=", len, a_len);
00280       return *this;
00281     }
00282 
00283   if (len == 0)
00284     return *this;
00285 
00286   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
00287 
00288   mx_inline_add2 (len, d, a.data ());
00289   return *this;
00290 }
00291 
00292 ComplexRowVector&
00293 ComplexRowVector::operator -= (const RowVector& a)
00294 {
00295   octave_idx_type len = length ();
00296 
00297   octave_idx_type a_len = a.length ();
00298 
00299   if (len != a_len)
00300     {
00301       gripe_nonconformant ("operator -=", len, a_len);
00302       return *this;
00303     }
00304 
00305   if (len == 0)
00306     return *this;
00307 
00308   Complex *d = fortran_vec (); // Ensures only one reference to my privates!
00309 
00310   mx_inline_sub2 (len, d, a.data ());
00311   return *this;
00312 }
00313 
00314 // row vector by matrix -> row vector
00315 
00316 ComplexRowVector
00317 operator * (const ComplexRowVector& v, const ComplexMatrix& a)
00318 {
00319   ComplexRowVector retval;
00320 
00321   octave_idx_type len = v.length ();
00322 
00323   octave_idx_type a_nr = a.rows ();
00324   octave_idx_type a_nc = a.cols ();
00325 
00326   if (a_nr != len)
00327     gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
00328   else
00329     {
00330       if (len == 0)
00331         retval.resize (a_nc, 0.0);
00332       else
00333         {
00334           // Transpose A to form A'*x == (x'*A)'
00335 
00336           octave_idx_type ld = a_nr;
00337 
00338           retval.resize (a_nc);
00339           Complex *y = retval.fortran_vec ();
00340 
00341           F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
00342                                    a_nr, a_nc, 1.0, a.data (),
00343                                    ld, v.data (), 1, 0.0, y, 1
00344                                    F77_CHAR_ARG_LEN (1)));
00345         }
00346     }
00347 
00348   return retval;
00349 }
00350 
00351 ComplexRowVector
00352 operator * (const RowVector& v, const ComplexMatrix& a)
00353 {
00354   ComplexRowVector tmp (v);
00355   return tmp * a;
00356 }
00357 
00358 // other operations
00359 
00360 Complex
00361 ComplexRowVector::min (void) const
00362 {
00363   octave_idx_type len = length ();
00364   if (len == 0)
00365     return Complex (0.0);
00366 
00367   Complex res = elem (0);
00368   double absres = std::abs (res);
00369 
00370   for (octave_idx_type i = 1; i < len; i++)
00371     if (std::abs (elem (i)) < absres)
00372       {
00373         res = elem (i);
00374         absres = std::abs (res);
00375       }
00376 
00377   return res;
00378 }
00379 
00380 Complex
00381 ComplexRowVector::max (void) const
00382 {
00383   octave_idx_type len = length ();
00384   if (len == 0)
00385     return Complex (0.0);
00386 
00387   Complex res = elem (0);
00388   double absres = std::abs (res);
00389 
00390   for (octave_idx_type i = 1; i < len; i++)
00391     if (std::abs (elem (i)) > absres)
00392       {
00393         res = elem (i);
00394         absres = std::abs (res);
00395       }
00396 
00397   return res;
00398 }
00399 
00400 // i/o
00401 
00402 std::ostream&
00403 operator << (std::ostream& os, const ComplexRowVector& a)
00404 {
00405 //  int field_width = os.precision () + 7;
00406   for (octave_idx_type i = 0; i < a.length (); i++)
00407     os << " " /* setw (field_width) */ << a.elem (i);
00408   return os;
00409 }
00410 
00411 std::istream&
00412 operator >> (std::istream& is, ComplexRowVector& a)
00413 {
00414   octave_idx_type len = a.length();
00415 
00416   if (len > 0)
00417     {
00418       Complex tmp;
00419       for (octave_idx_type i = 0; i < len; i++)
00420         {
00421           is >> tmp;
00422           if (is)
00423             a.elem (i) = tmp;
00424           else
00425             break;
00426         }
00427     }
00428   return is;
00429 }
00430 
00431 // row vector by column vector -> scalar
00432 
00433 // row vector by column vector -> scalar
00434 
00435 Complex
00436 operator * (const ComplexRowVector& v, const ColumnVector& a)
00437 {
00438   ComplexColumnVector tmp (a);
00439   return v * tmp;
00440 }
00441 
00442 Complex
00443 operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
00444 {
00445   Complex retval (0.0, 0.0);
00446 
00447   octave_idx_type len = v.length ();
00448 
00449   octave_idx_type a_len = a.length ();
00450 
00451   if (len != a_len)
00452     gripe_nonconformant ("operator *", len, a_len);
00453   else if (len != 0)
00454     F77_FUNC (xzdotu, XZDOTU) (len, v.data (), 1, a.data (), 1, retval);
00455 
00456   return retval;
00457 }
00458 
00459 // other operations
00460 
00461 ComplexRowVector
00462 linspace (const Complex& x1, const Complex& x2, octave_idx_type n)
00463 {
00464   if (n < 1) n = 1;
00465 
00466   NoAlias<ComplexRowVector> retval (n);
00467 
00468   Complex delta = (x2 - x1) / (n - 1.0);
00469   retval(0) = x1;
00470   for (octave_idx_type i = 1; i < n-1; i++)
00471     retval(i) = x1 + static_cast<double> (i)*delta;
00472   retval(n-1) = x2;
00473 
00474   return retval;
00475 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines