00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <iostream>
00033
00034 #include "Array-util.h"
00035 #include "f77-fcn.h"
00036 #include "lo-error.h"
00037 #include "mx-base.h"
00038 #include "mx-inlines.cc"
00039 #include "oct-cmplx.h"
00040
00041
00042
00043 extern "C"
00044 {
00045 F77_RET_T
00046 F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
00047 const int&, const int&, const double&,
00048 const double*, const int&, const double*,
00049 const int&, const double&, double*, const int&
00050 F77_CHAR_ARG_LEN_DECL);
00051
00052 double F77_FUNC (ddot, DDOT) (const int&, const double*, const int&,
00053 const double*, const int&);
00054 }
00055
00056
00057
00058 bool
00059 RowVector::operator == (const RowVector& a) const
00060 {
00061 int len = length ();
00062 if (len != a.length ())
00063 return 0;
00064 return mx_inline_equal (data (), a.data (), len);
00065 }
00066
00067 bool
00068 RowVector::operator != (const RowVector& a) const
00069 {
00070 return !(*this == a);
00071 }
00072
00073 RowVector&
00074 RowVector::insert (const RowVector& a, int c)
00075 {
00076 int a_len = a.length ();
00077
00078 if (c < 0 || c + a_len > length ())
00079 {
00080 (*current_liboctave_error_handler) ("range error for insert");
00081 return *this;
00082 }
00083
00084 if (a_len > 0)
00085 {
00086 make_unique ();
00087
00088 for (int i = 0; i < a_len; i++)
00089 xelem (c+i) = a.elem (i);
00090 }
00091
00092 return *this;
00093 }
00094
00095 RowVector&
00096 RowVector::fill (double val)
00097 {
00098 int len = length ();
00099
00100 if (len > 0)
00101 {
00102 make_unique ();
00103
00104 for (int i = 0; i < len; i++)
00105 xelem (i) = val;
00106 }
00107
00108 return *this;
00109 }
00110
00111 RowVector&
00112 RowVector::fill (double val, int c1, int c2)
00113 {
00114 int len = length ();
00115
00116 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
00117 {
00118 (*current_liboctave_error_handler) ("range error for fill");
00119 return *this;
00120 }
00121
00122 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00123
00124 if (c2 >= c1)
00125 {
00126 make_unique ();
00127
00128 for (int i = c1; i <= c2; i++)
00129 xelem (i) = val;
00130 }
00131
00132 return *this;
00133 }
00134
00135 RowVector
00136 RowVector::append (const RowVector& a) const
00137 {
00138 int len = length ();
00139 int nc_insert = len;
00140 RowVector retval (len + a.length ());
00141 retval.insert (*this, 0);
00142 retval.insert (a, nc_insert);
00143 return retval;
00144 }
00145
00146 ColumnVector
00147 RowVector::transpose (void) const
00148 {
00149 return ColumnVector (*this);
00150 }
00151
00152 RowVector
00153 real (const ComplexRowVector& a)
00154 {
00155 int a_len = a.length ();
00156 RowVector retval;
00157 if (a_len > 0)
00158 retval = RowVector (mx_inline_real_dup (a.data (), a_len), a_len);
00159 return retval;
00160 }
00161
00162 RowVector
00163 imag (const ComplexRowVector& a)
00164 {
00165 int a_len = a.length ();
00166 RowVector retval;
00167 if (a_len > 0)
00168 retval = RowVector (mx_inline_imag_dup (a.data (), a_len), a_len);
00169 return retval;
00170 }
00171
00172 RowVector
00173 RowVector::extract (int c1, int c2) const
00174 {
00175 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00176
00177 int new_c = c2 - c1 + 1;
00178
00179 RowVector result (new_c);
00180
00181 for (int i = 0; i < new_c; i++)
00182 result.xelem (i) = elem (c1+i);
00183
00184 return result;
00185 }
00186
00187 RowVector
00188 RowVector::extract_n (int r1, int n) const
00189 {
00190 RowVector result (n);
00191
00192 for (int i = 0; i < n; i++)
00193 result.xelem (i) = elem (r1+i);
00194
00195 return result;
00196 }
00197
00198
00199
00200 RowVector
00201 operator * (const RowVector& v, const Matrix& a)
00202 {
00203 RowVector retval;
00204
00205 int len = v.length ();
00206
00207 int a_nr = a.rows ();
00208 int a_nc = a.cols ();
00209
00210 if (a_nr != len)
00211 gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
00212 else
00213 {
00214 if (len == 0)
00215 retval.resize (a_nc, 0.0);
00216 else
00217 {
00218
00219
00220 int ld = a_nr;
00221
00222 retval.resize (a_nc);
00223 double *y = retval.fortran_vec ();
00224
00225 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
00226 a_nr, a_nc, 1.0, a.data (),
00227 ld, v.data (), 1, 0.0, y, 1
00228 F77_CHAR_ARG_LEN (1)));
00229
00230 if (f77_exception_encountered)
00231 (*current_liboctave_error_handler)
00232 ("unrecoverable error in dgemv");
00233 }
00234 }
00235
00236 return retval;
00237 }
00238
00239
00240
00241 RowVector
00242 RowVector::map (d_d_Mapper f) const
00243 {
00244 RowVector b (*this);
00245 return b.apply (f);
00246 }
00247
00248 RowVector&
00249 RowVector::apply (d_d_Mapper f)
00250 {
00251 double *d = fortran_vec ();
00252
00253 for (int i = 0; i < length (); i++)
00254 d[i] = f (d[i]);
00255
00256 return *this;
00257 }
00258
00259 double
00260 RowVector::min (void) const
00261 {
00262 int len = length ();
00263 if (len == 0)
00264 return 0;
00265
00266 double res = elem (0);
00267
00268 for (int i = 1; i < len; i++)
00269 if (elem (i) < res)
00270 res = elem (i);
00271
00272 return res;
00273 }
00274
00275 double
00276 RowVector::max (void) const
00277 {
00278 int len = length ();
00279 if (len == 0)
00280 return 0;
00281
00282 double res = elem (0);
00283
00284 for (int i = 1; i < len; i++)
00285 if (elem (i) > res)
00286 res = elem (i);
00287
00288 return res;
00289 }
00290
00291 std::ostream&
00292 operator << (std::ostream& os, const RowVector& a)
00293 {
00294
00295
00296 for (int i = 0; i < a.length (); i++)
00297 os << " " << a.elem (i);
00298 return os;
00299 }
00300
00301 std::istream&
00302 operator >> (std::istream& is, RowVector& a)
00303 {
00304 int len = a.length();
00305
00306 if (len < 1)
00307 is.clear (std::ios::badbit);
00308 else
00309 {
00310 double tmp;
00311 for (int i = 0; i < len; i++)
00312 {
00313 is >> tmp;
00314 if (is)
00315 a.elem (i) = tmp;
00316 else
00317 break;
00318 }
00319 }
00320 return is;
00321 }
00322
00323
00324
00325 RowVector
00326 linspace (double x1, double x2, int n)
00327 {
00328 RowVector retval;
00329
00330 if (n > 1)
00331 {
00332 retval.resize (n);
00333 double delta = (x2 - x1) / (n - 1);
00334 retval.elem (0) = x1;
00335 for (int i = 1; i < n-1; i++)
00336 retval.elem (i) = x1 + i * delta;
00337 retval.elem (n-1) = x2;
00338 }
00339 else if (n == 1)
00340 {
00341 if (x1 == x2)
00342 {
00343 retval.resize (1);
00344 retval.elem (0) = x1;
00345 }
00346 else
00347 (*current_liboctave_error_handler)
00348 ("linspace: npoints is 1, but x1 != x2");
00349 }
00350 else
00351 (*current_liboctave_error_handler)
00352 ("linspace: npoints must be greater than 0");
00353
00354 return retval;
00355 }
00356
00357
00358
00359 double
00360 operator * (const RowVector& v, const ColumnVector& a)
00361 {
00362 double retval = 0.0;
00363
00364 int len = v.length ();
00365
00366 int a_len = a.length ();
00367
00368 if (len != a_len)
00369 gripe_nonconformant ("operator *", len, a_len);
00370 else if (len != 0)
00371 retval = F77_FUNC (ddot, DDOT) (len, v.data (), 1, a.data (), 1);
00372
00373 return retval;
00374 }
00375
00376 Complex
00377 operator * (const RowVector& v, const ComplexColumnVector& a)
00378 {
00379 ComplexRowVector tmp (v);
00380 return tmp * a;
00381 }
00382
00383
00384
00385
00386
00387