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 ComplexColumnVector::ComplexColumnVector (const ColumnVector& 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 ComplexColumnVector::operator == (const ComplexColumnVector& 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 ComplexColumnVector::operator != (const ComplexColumnVector& a) const
00073 {
00074 return !(*this == a);
00075 }
00076
00077
00078
00079 ComplexColumnVector&
00080 ComplexColumnVector::insert (const ColumnVector& a, int r)
00081 {
00082 int a_len = a.length ();
00083
00084 if (r < 0 || r + 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 (r+i) = a.elem (i);
00096 }
00097
00098 return *this;
00099 }
00100
00101 ComplexColumnVector&
00102 ComplexColumnVector::insert (const ComplexColumnVector& a, int r)
00103 {
00104 int a_len = a.length ();
00105
00106 if (r < 0 || r + 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 (r+i) = a.elem (i);
00118 }
00119
00120 return *this;
00121 }
00122
00123 ComplexColumnVector&
00124 ComplexColumnVector::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 ComplexColumnVector&
00140 ComplexColumnVector::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
00153 return *this;
00154 }
00155
00156 ComplexColumnVector&
00157 ComplexColumnVector::fill (double val, int r1, int r2)
00158 {
00159 int len = length ();
00160
00161 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
00162 {
00163 (*current_liboctave_error_handler) ("range error for fill");
00164 return *this;
00165 }
00166
00167 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00168
00169 if (r2 >= r1)
00170 {
00171 make_unique ();
00172
00173 for (int i = r1; i <= r2; i++)
00174 xelem (i) = val;
00175 }
00176
00177 return *this;
00178 }
00179
00180 ComplexColumnVector&
00181 ComplexColumnVector::fill (const Complex& val, int r1, int r2)
00182 {
00183 int len = length ();
00184
00185 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
00186 {
00187 (*current_liboctave_error_handler) ("range error for fill");
00188 return *this;
00189 }
00190
00191 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00192
00193 if (r2 >= r1)
00194 {
00195 make_unique ();
00196
00197 for (int i = r1; i <= r2; i++)
00198 xelem (i) = val;
00199 }
00200
00201 return *this;
00202 }
00203
00204 ComplexColumnVector
00205 ComplexColumnVector::stack (const ColumnVector& a) const
00206 {
00207 int len = length ();
00208 int nr_insert = len;
00209 ComplexColumnVector retval (len + a.length ());
00210 retval.insert (*this, 0);
00211 retval.insert (a, nr_insert);
00212 return retval;
00213 }
00214
00215 ComplexColumnVector
00216 ComplexColumnVector::stack (const ComplexColumnVector& a) const
00217 {
00218 int len = length ();
00219 int nr_insert = len;
00220 ComplexColumnVector retval (len + a.length ());
00221 retval.insert (*this, 0);
00222 retval.insert (a, nr_insert);
00223 return retval;
00224 }
00225
00226 ComplexRowVector
00227 ComplexColumnVector::hermitian (void) const
00228 {
00229 int len = length ();
00230 return ComplexRowVector (mx_inline_conj_dup (data (), len), len);
00231 }
00232
00233 ComplexRowVector
00234 ComplexColumnVector::transpose (void) const
00235 {
00236 return ComplexRowVector (*this);
00237 }
00238
00239 ComplexColumnVector
00240 conj (const ComplexColumnVector& a)
00241 {
00242 int a_len = a.length ();
00243 ComplexColumnVector retval;
00244 if (a_len > 0)
00245 retval = ComplexColumnVector (mx_inline_conj_dup (a.data (), a_len), a_len);
00246 return retval;
00247 }
00248
00249
00250
00251 ComplexColumnVector
00252 ComplexColumnVector::extract (int r1, int r2) const
00253 {
00254 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00255
00256 int new_r = r2 - r1 + 1;
00257
00258 ComplexColumnVector result (new_r);
00259
00260 for (int i = 0; i < new_r; i++)
00261 result.elem (i) = elem (r1+i);
00262
00263 return result;
00264 }
00265
00266 ComplexColumnVector
00267 ComplexColumnVector::extract_n (int r1, int n) const
00268 {
00269 ComplexColumnVector result (n);
00270
00271 for (int i = 0; i < n; i++)
00272 result.elem (i) = elem (r1+i);
00273
00274 return result;
00275 }
00276
00277
00278
00279 ComplexColumnVector&
00280 ComplexColumnVector::operator += (const ColumnVector& a)
00281 {
00282 int len = length ();
00283
00284 int a_len = a.length ();
00285
00286 if (len != a_len)
00287 {
00288 gripe_nonconformant ("operator +=", len, a_len);
00289 return *this;
00290 }
00291
00292 if (len == 0)
00293 return *this;
00294
00295 Complex *d = fortran_vec ();
00296
00297 mx_inline_add2 (d, a.data (), len);
00298 return *this;
00299 }
00300
00301 ComplexColumnVector&
00302 ComplexColumnVector::operator -= (const ColumnVector& a)
00303 {
00304 int len = length ();
00305
00306 int a_len = a.length ();
00307
00308 if (len != a_len)
00309 {
00310 gripe_nonconformant ("operator -=", len, a_len);
00311 return *this;
00312 }
00313
00314 if (len == 0)
00315 return *this;
00316
00317 Complex *d = fortran_vec ();
00318
00319 mx_inline_subtract2 (d, a.data (), len);
00320 return *this;
00321 }
00322
00323
00324
00325 ComplexColumnVector
00326 operator * (const ComplexMatrix& m, const ColumnVector& a)
00327 {
00328 ComplexColumnVector tmp (a);
00329 return m * tmp;
00330 }
00331
00332 ComplexColumnVector
00333 operator * (const ComplexMatrix& m, const ComplexColumnVector& a)
00334 {
00335 ComplexColumnVector retval;
00336
00337 int nr = m.rows ();
00338 int nc = m.cols ();
00339
00340 int a_len = a.length ();
00341
00342 if (nc != a_len)
00343 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00344 else
00345 {
00346 if (nc == 0 || nr == 0)
00347 retval.resize (nr, 0.0);
00348 else
00349 {
00350 int ld = nr;
00351
00352 retval.resize (nr);
00353 Complex *y = retval.fortran_vec ();
00354
00355 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
00356 nr, nc, 1.0, m.data (), ld,
00357 a.data (), 1, 0.0, y, 1
00358 F77_CHAR_ARG_LEN (1)));
00359
00360 if (f77_exception_encountered)
00361 (*current_liboctave_error_handler)
00362 ("unrecoverable error in zgemv");
00363 }
00364 }
00365
00366 return retval;
00367 }
00368
00369
00370
00371 ComplexColumnVector
00372 operator * (const Matrix& m, const ComplexColumnVector& a)
00373 {
00374 ComplexMatrix tmp (m);
00375 return tmp * a;
00376 }
00377
00378
00379
00380 ComplexColumnVector
00381 operator * (const DiagMatrix& m, const ComplexColumnVector& a)
00382 {
00383 int nr = m.rows ();
00384 int nc = m.cols ();
00385
00386 int a_len = a.length ();
00387
00388 if (nc != a_len)
00389 {
00390 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00391 return ComplexColumnVector ();
00392 }
00393
00394 if (nc == 0 || nr == 0)
00395 return ComplexColumnVector (0);
00396
00397 ComplexColumnVector result (nr);
00398
00399 for (int i = 0; i < a_len; i++)
00400 result.elem (i) = a.elem (i) * m.elem (i, i);
00401
00402 for (int i = a_len; i < nr; i++)
00403 result.elem (i) = 0.0;
00404
00405 return result;
00406 }
00407
00408 ComplexColumnVector
00409 operator * (const ComplexDiagMatrix& m, const ColumnVector& a)
00410 {
00411 int nr = m.rows ();
00412 int nc = m.cols ();
00413
00414 int a_len = a.length ();
00415
00416 if (nc != a_len)
00417 {
00418 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00419 return ComplexColumnVector ();
00420 }
00421
00422 if (nc == 0 || nr == 0)
00423 return ComplexColumnVector (0);
00424
00425 ComplexColumnVector result (nr);
00426
00427 for (int i = 0; i < a_len; i++)
00428 result.elem (i) = a.elem (i) * m.elem (i, i);
00429
00430 for (int i = a_len; i < nr; i++)
00431 result.elem (i) = 0.0;
00432
00433 return result;
00434 }
00435
00436 ComplexColumnVector
00437 operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a)
00438 {
00439 int nr = m.rows ();
00440 int nc = m.cols ();
00441
00442 int a_len = a.length ();
00443
00444 if (nc != a_len)
00445 {
00446 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00447 return ComplexColumnVector ();
00448 }
00449
00450 if (nc == 0 || nr == 0)
00451 return ComplexColumnVector (0);
00452
00453 ComplexColumnVector result (nr);
00454
00455 for (int i = 0; i < a_len; i++)
00456 result.elem (i) = a.elem (i) * m.elem (i, i);
00457
00458 for (int i = a_len; i < nr; i++)
00459 result.elem (i) = 0.0;
00460
00461 return result;
00462 }
00463
00464
00465
00466 ComplexColumnVector
00467 ComplexColumnVector::map (c_c_Mapper f) const
00468 {
00469 ComplexColumnVector b (*this);
00470 return b.apply (f);
00471 }
00472
00473 ColumnVector
00474 ComplexColumnVector::map (d_c_Mapper f) const
00475 {
00476 const Complex *d = data ();
00477
00478 int len = length ();
00479
00480 ColumnVector retval (len);
00481
00482 double *r = retval.fortran_vec ();
00483
00484 for (int i = 0; i < len; i++)
00485 r[i] = f (d[i]);
00486
00487 return retval;
00488 }
00489
00490 ComplexColumnVector&
00491 ComplexColumnVector::apply (c_c_Mapper f)
00492 {
00493 Complex *d = fortran_vec ();
00494
00495 for (int i = 0; i < length (); i++)
00496 d[i] = f (d[i]);
00497
00498 return *this;
00499 }
00500
00501 Complex
00502 ComplexColumnVector::min (void) const
00503 {
00504 int len = length ();
00505 if (len == 0)
00506 return 0.0;
00507
00508 Complex res = elem (0);
00509 double absres = abs (res);
00510
00511 for (int i = 1; i < len; i++)
00512 if (abs (elem (i)) < absres)
00513 {
00514 res = elem (i);
00515 absres = abs (res);
00516 }
00517
00518 return res;
00519 }
00520
00521 Complex
00522 ComplexColumnVector::max (void) const
00523 {
00524 int len = length ();
00525 if (len == 0)
00526 return 0.0;
00527
00528 Complex res = elem (0);
00529 double absres = abs (res);
00530
00531 for (int i = 1; i < len; i++)
00532 if (abs (elem (i)) > absres)
00533 {
00534 res = elem (i);
00535 absres = abs (res);
00536 }
00537
00538 return res;
00539 }
00540
00541
00542
00543 std::ostream&
00544 operator << (std::ostream& os, const ComplexColumnVector& a)
00545 {
00546
00547 for (int i = 0; i < a.length (); i++)
00548 os << a.elem (i) << "\n";
00549 return os;
00550 }
00551
00552 std::istream&
00553 operator >> (std::istream& is, ComplexColumnVector& a)
00554 {
00555 int len = a.length();
00556
00557 if (len < 1)
00558 is.clear (std::ios::badbit);
00559 else
00560 {
00561 double tmp;
00562 for (int i = 0; i < len; i++)
00563 {
00564 is >> tmp;
00565 if (is)
00566 a.elem (i) = tmp;
00567 else
00568 break;
00569 }
00570 }
00571 return is;
00572 }
00573
00574
00575
00576
00577
00578