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 (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
00047 const int&, const int&, const Complex&,
00048 const Complex*, const int&, const Complex*,
00049 const int&, const Complex&, Complex*, const int&
00050 F77_CHAR_ARG_LEN_DECL);
00051 }
00052
00053
00054
00055 ComplexRowVector::ComplexRowVector (const RowVector& a)
00056 : MArray<Complex> (a.length ())
00057 {
00058 for (int i = 0; i < length (); i++)
00059 elem (i) = a.elem (i);
00060 }
00061
00062 bool
00063 ComplexRowVector::operator == (const ComplexRowVector& a) const
00064 {
00065 int len = length ();
00066 if (len != a.length ())
00067 return 0;
00068 return mx_inline_equal (data (), a.data (), len);
00069 }
00070
00071 bool
00072 ComplexRowVector::operator != (const ComplexRowVector& a) const
00073 {
00074 return !(*this == a);
00075 }
00076
00077
00078
00079 ComplexRowVector&
00080 ComplexRowVector::insert (const RowVector& a, int c)
00081 {
00082 int a_len = a.length ();
00083
00084 if (c < 0 || c + a_len > length ())
00085 {
00086 (*current_liboctave_error_handler) ("range error for insert");
00087 return *this;
00088 }
00089
00090 if (a_len > 0)
00091 {
00092 make_unique ();
00093
00094 for (int i = 0; i < a_len; i++)
00095 xelem (c+i) = a.elem (i);
00096 }
00097
00098 return *this;
00099 }
00100
00101 ComplexRowVector&
00102 ComplexRowVector::insert (const ComplexRowVector& a, int c)
00103 {
00104 int a_len = a.length ();
00105
00106 if (c < 0 || c + a_len > length ())
00107 {
00108 (*current_liboctave_error_handler) ("range error for insert");
00109 return *this;
00110 }
00111
00112 if (a_len > 0)
00113 {
00114 make_unique ();
00115
00116 for (int i = 0; i < a_len; i++)
00117 xelem (c+i) = a.elem (i);
00118 }
00119
00120 return *this;
00121 }
00122
00123 ComplexRowVector&
00124 ComplexRowVector::fill (double val)
00125 {
00126 int len = length ();
00127
00128 if (len > 0)
00129 {
00130 make_unique ();
00131
00132 for (int i = 0; i < len; i++)
00133 xelem (i) = val;
00134 }
00135
00136 return *this;
00137 }
00138
00139 ComplexRowVector&
00140 ComplexRowVector::fill (const Complex& val)
00141 {
00142 int len = length ();
00143
00144 if (len > 0)
00145 {
00146 make_unique ();
00147
00148 for (int i = 0; i < len; i++)
00149 xelem (i) = val;
00150 }
00151
00152 return *this;
00153 }
00154
00155 ComplexRowVector&
00156 ComplexRowVector::fill (double val, int c1, int c2)
00157 {
00158 int len = length ();
00159
00160 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
00161 {
00162 (*current_liboctave_error_handler) ("range error for fill");
00163 return *this;
00164 }
00165
00166 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00167
00168 if (c2 >= c1)
00169 {
00170 make_unique ();
00171
00172 for (int i = c1; i <= c2; i++)
00173 xelem (i) = val;
00174 }
00175
00176 return *this;
00177 }
00178
00179 ComplexRowVector&
00180 ComplexRowVector::fill (const Complex& val, int c1, int c2)
00181 {
00182 int len = length ();
00183
00184 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
00185 {
00186 (*current_liboctave_error_handler) ("range error for fill");
00187 return *this;
00188 }
00189
00190 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00191
00192 if (c2 >= c1)
00193 {
00194 make_unique ();
00195
00196 for (int i = c1; i <= c2; i++)
00197 xelem (i) = val;
00198 }
00199
00200 return *this;
00201 }
00202
00203 ComplexRowVector
00204 ComplexRowVector::append (const RowVector& a) const
00205 {
00206 int len = length ();
00207 int nc_insert = len;
00208 ComplexRowVector retval (len + a.length ());
00209 retval.insert (*this, 0);
00210 retval.insert (a, nc_insert);
00211 return retval;
00212 }
00213
00214 ComplexRowVector
00215 ComplexRowVector::append (const ComplexRowVector& a) const
00216 {
00217 int len = length ();
00218 int nc_insert = len;
00219 ComplexRowVector retval (len + a.length ());
00220 retval.insert (*this, 0);
00221 retval.insert (a, nc_insert);
00222 return retval;
00223 }
00224
00225 ComplexColumnVector
00226 ComplexRowVector::hermitian (void) const
00227 {
00228 int len = length ();
00229 return ComplexColumnVector (mx_inline_conj_dup (data (), len), len);
00230 }
00231
00232 ComplexColumnVector
00233 ComplexRowVector::transpose (void) const
00234 {
00235 return ComplexColumnVector (*this);
00236 }
00237
00238 ComplexRowVector
00239 conj (const ComplexRowVector& a)
00240 {
00241 int a_len = a.length ();
00242 ComplexRowVector retval;
00243 if (a_len > 0)
00244 retval = ComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len);
00245 return retval;
00246 }
00247
00248
00249
00250 ComplexRowVector
00251 ComplexRowVector::extract (int c1, int c2) const
00252 {
00253 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00254
00255 int new_c = c2 - c1 + 1;
00256
00257 ComplexRowVector result (new_c);
00258
00259 for (int i = 0; i < new_c; i++)
00260 result.elem (i) = elem (c1+i);
00261
00262 return result;
00263 }
00264
00265 ComplexRowVector
00266 ComplexRowVector::extract_n (int r1, int n) const
00267 {
00268 ComplexRowVector result (n);
00269
00270 for (int i = 0; i < n; i++)
00271 result.elem (i) = elem (r1+i);
00272
00273 return result;
00274 }
00275
00276
00277
00278 ComplexRowVector&
00279 ComplexRowVector::operator += (const RowVector& a)
00280 {
00281 int len = length ();
00282
00283 int a_len = a.length ();
00284
00285 if (len != a_len)
00286 {
00287 gripe_nonconformant ("operator +=", len, a_len);
00288 return *this;
00289 }
00290
00291 if (len == 0)
00292 return *this;
00293
00294 Complex *d = fortran_vec ();
00295
00296 mx_inline_add2 (d, a.data (), len);
00297 return *this;
00298 }
00299
00300 ComplexRowVector&
00301 ComplexRowVector::operator -= (const RowVector& a)
00302 {
00303 int len = length ();
00304
00305 int a_len = a.length ();
00306
00307 if (len != a_len)
00308 {
00309 gripe_nonconformant ("operator -=", len, a_len);
00310 return *this;
00311 }
00312
00313 if (len == 0)
00314 return *this;
00315
00316 Complex *d = fortran_vec ();
00317
00318 mx_inline_subtract2 (d, a.data (), len);
00319 return *this;
00320 }
00321
00322
00323
00324 ComplexRowVector
00325 operator * (const ComplexRowVector& v, const ComplexMatrix& a)
00326 {
00327 ComplexRowVector retval;
00328
00329 int len = v.length ();
00330
00331 int a_nr = a.rows ();
00332 int a_nc = a.cols ();
00333
00334 if (a_nr != len)
00335 gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
00336 else
00337 {
00338 if (len == 0)
00339 retval.resize (a_nc, 0.0);
00340 else
00341 {
00342
00343
00344 int ld = a_nr;
00345
00346 retval.resize (a_nc);
00347 Complex *y = retval.fortran_vec ();
00348
00349 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
00350 a_nr, a_nc, 1.0, a.data (),
00351 ld, v.data (), 1, 0.0, y, 1
00352 F77_CHAR_ARG_LEN (1)));
00353
00354 if (f77_exception_encountered)
00355 (*current_liboctave_error_handler)
00356 ("unrecoverable error in zgemv");
00357 }
00358 }
00359
00360 return retval;
00361 }
00362
00363 ComplexRowVector
00364 operator * (const RowVector& v, const ComplexMatrix& a)
00365 {
00366 ComplexRowVector tmp (v);
00367 return tmp * a;
00368 }
00369
00370
00371
00372 ComplexRowVector
00373 ComplexRowVector::map (c_c_Mapper f) const
00374 {
00375 ComplexRowVector b (*this);
00376 return b.apply (f);
00377 }
00378
00379 RowVector
00380 ComplexRowVector::map (d_c_Mapper f) const
00381 {
00382 const Complex *d = data ();
00383
00384 int len = length ();
00385
00386 RowVector retval (len);
00387
00388 double *r = retval.fortran_vec ();
00389
00390 for (int i = 0; i < len; i++)
00391 r[i] = f (d[i]);
00392
00393 return retval;
00394 }
00395
00396 ComplexRowVector&
00397 ComplexRowVector::apply (c_c_Mapper f)
00398 {
00399 Complex *d = fortran_vec ();
00400
00401 for (int i = 0; i < length (); i++)
00402 d[i] = f (d[i]);
00403
00404 return *this;
00405 }
00406
00407 Complex
00408 ComplexRowVector::min (void) const
00409 {
00410 int len = length ();
00411 if (len == 0)
00412 return Complex (0.0);
00413
00414 Complex res = elem (0);
00415 double absres = abs (res);
00416
00417 for (int i = 1; i < len; i++)
00418 if (abs (elem (i)) < absres)
00419 {
00420 res = elem (i);
00421 absres = abs (res);
00422 }
00423
00424 return res;
00425 }
00426
00427 Complex
00428 ComplexRowVector::max (void) const
00429 {
00430 int len = length ();
00431 if (len == 0)
00432 return Complex (0.0);
00433
00434 Complex res = elem (0);
00435 double absres = abs (res);
00436
00437 for (int i = 1; i < len; i++)
00438 if (abs (elem (i)) > absres)
00439 {
00440 res = elem (i);
00441 absres = abs (res);
00442 }
00443
00444 return res;
00445 }
00446
00447
00448
00449 std::ostream&
00450 operator << (std::ostream& os, const ComplexRowVector& a)
00451 {
00452
00453 for (int i = 0; i < a.length (); i++)
00454 os << " " << a.elem (i);
00455 return os;
00456 }
00457
00458 std::istream&
00459 operator >> (std::istream& is, ComplexRowVector& a)
00460 {
00461 int len = a.length();
00462
00463 if (len < 1)
00464 is.clear (std::ios::badbit);
00465 else
00466 {
00467 Complex tmp;
00468 for (int i = 0; i < len; i++)
00469 {
00470 is >> tmp;
00471 if (is)
00472 a.elem (i) = tmp;
00473 else
00474 break;
00475 }
00476 }
00477 return is;
00478 }
00479
00480
00481
00482
00483
00484 Complex
00485 operator * (const ComplexRowVector& v, const ColumnVector& a)
00486 {
00487 ComplexColumnVector tmp (a);
00488 return v * tmp;
00489 }
00490
00491 Complex
00492 operator * (const ComplexRowVector& v, const ComplexColumnVector& a)
00493 {
00494 int len = v.length ();
00495
00496 int a_len = a.length ();
00497
00498 if (len != a_len)
00499 {
00500 gripe_nonconformant ("operator *", len, a_len);
00501 return 0.0;
00502 }
00503
00504 Complex retval (0.0, 0.0);
00505
00506 for (int i = 0; i < len; i++)
00507 retval += v.elem (i) * a.elem (i);
00508
00509 return retval;
00510 }
00511
00512
00513
00514 ComplexRowVector
00515 linspace (const Complex& x1, const Complex& x2, int n)
00516 {
00517 ComplexRowVector retval;
00518
00519 if (n > 0)
00520 {
00521 retval.resize (n);
00522 Complex delta = (x2 - x1) / (n - 1.0);
00523 retval.elem (0) = x1;
00524 for (int i = 1; i < n-1; i++)
00525 retval.elem (i) = x1 + 1.0 * i * delta;
00526 retval.elem (n-1) = x2;
00527 }
00528 else if (n == 1)
00529 {
00530 if (x1 == x2)
00531 {
00532 retval.resize (1);
00533 retval.elem (0) = x1;
00534 }
00535 else
00536 (*current_liboctave_error_handler)
00537 ("linspace: npoints is 1, but x1 != x2");
00538 }
00539 else
00540 (*current_liboctave_error_handler)
00541 ("linspace: npoints must be greater than 0");
00542
00543 return retval;
00544 }
00545
00546
00547
00548
00549
00550