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 <cfloat>
00033
00034 #include <iostream>
00035
00036
00037 #ifdef HAVE_SYS_TYPES_H
00038 #include <sys/types.h>
00039 #endif
00040
00041 #include "Array-util.h"
00042 #include "CMatrix.h"
00043 #include "CmplxAEPBAL.h"
00044 #include "CmplxDET.h"
00045 #include "CmplxSCHUR.h"
00046 #include "CmplxSVD.h"
00047 #include "f77-fcn.h"
00048 #include "lo-error.h"
00049 #include "lo-ieee.h"
00050 #include "lo-mappers.h"
00051 #include "lo-utils.h"
00052 #include "mx-base.h"
00053 #include "mx-cm-dm.h"
00054 #include "mx-dm-cm.h"
00055 #include "mx-cm-s.h"
00056 #include "mx-inlines.cc"
00057 #include "oct-cmplx.h"
00058
00059 #if defined (HAVE_FFTW3)
00060 #include "oct-fftw.h"
00061 #endif
00062
00063
00064
00065 extern "C"
00066 {
00067 F77_RET_T
00068 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
00069 const int&, Complex*, const int&, int&,
00070 int&, double*, int&
00071 F77_CHAR_ARG_LEN_DECL);
00072
00073 F77_RET_T
00074 F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
00075 F77_CONST_CHAR_ARG_DECL,
00076 const int&, const int&, const int&, double*,
00077 const int&, double*, const int&, int&
00078 F77_CHAR_ARG_LEN_DECL
00079 F77_CHAR_ARG_LEN_DECL);
00080
00081 F77_RET_T
00082 F77_FUNC (zgemm, ZGEMM) (F77_CONST_CHAR_ARG_DECL,
00083 F77_CONST_CHAR_ARG_DECL,
00084 const int&, const int&, const int&,
00085 const Complex&, const Complex*, const int&,
00086 const Complex*, const int&, const Complex&,
00087 Complex*, const int&
00088 F77_CHAR_ARG_LEN_DECL
00089 F77_CHAR_ARG_LEN_DECL);
00090
00091 F77_RET_T
00092 F77_FUNC (zgetrf, ZGETRF) (const int&, const int&, Complex*, const int&,
00093 int*, int&);
00094
00095 F77_RET_T
00096 F77_FUNC (zgetrs, ZGETRS) (F77_CONST_CHAR_ARG_DECL,
00097 const int&, const int&, Complex*, const int&,
00098 const int*, Complex*, const int&, int&
00099 F77_CHAR_ARG_LEN_DECL);
00100
00101 F77_RET_T
00102 F77_FUNC (zgetri, ZGETRI) (const int&, Complex*, const int&, const int*,
00103 Complex*, const int&, int&);
00104
00105 F77_RET_T
00106 F77_FUNC (zgecon, ZGECON) (F77_CONST_CHAR_ARG_DECL,
00107 const int&, Complex*,
00108 const int&, const double&, double&,
00109 Complex*, double*, int&
00110 F77_CHAR_ARG_LEN_DECL);
00111
00112 F77_RET_T
00113 F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&,
00114 Complex*, const int&, Complex*,
00115 const int&, double*, double&, int&,
00116 Complex*, const int&, double*, int&);
00117
00118
00119
00120
00121
00122
00123 F77_RET_T
00124 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
00125
00126 F77_RET_T
00127 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
00128
00129 F77_RET_T
00130 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
00131
00132 F77_RET_T
00133 F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&,
00134 double&, Complex&, Complex&);
00135
00136 F77_RET_T
00137 F77_FUNC (ztrsyl, ZTRSYL) (F77_CONST_CHAR_ARG_DECL,
00138 F77_CONST_CHAR_ARG_DECL,
00139 const int&, const int&, const int&,
00140 const Complex*, const int&,
00141 const Complex*, const int&,
00142 const Complex*, const int&, double&, int&
00143 F77_CHAR_ARG_LEN_DECL
00144 F77_CHAR_ARG_LEN_DECL);
00145
00146 F77_RET_T
00147 F77_FUNC (xzlange, XZLANGE) (F77_CONST_CHAR_ARG_DECL,
00148 const int&, const int&, const Complex*,
00149 const int&, double*, double&
00150 F77_CHAR_ARG_LEN_DECL);
00151 }
00152
00153 static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
00154
00155
00156
00157 ComplexMatrix::ComplexMatrix (const Matrix& a)
00158 : MArray2<Complex> (a.rows (), a.cols ())
00159 {
00160 for (int j = 0; j < cols (); j++)
00161 for (int i = 0; i < rows (); i++)
00162 elem (i, j) = a.elem (i, j);
00163 }
00164
00165 ComplexMatrix::ComplexMatrix (const RowVector& rv)
00166 : MArray2<Complex> (1, rv.length (), 0.0)
00167 {
00168 for (int i = 0; i < rv.length (); i++)
00169 elem (0, i) = rv.elem (i);
00170 }
00171
00172 ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
00173 : MArray2<Complex> (cv.length (), 1, 0.0)
00174 {
00175 for (int i = 0; i < cv.length (); i++)
00176 elem (i, 0) = cv.elem (i);
00177 }
00178
00179 ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
00180 : MArray2<Complex> (a.rows (), a.cols (), 0.0)
00181 {
00182 for (int i = 0; i < a.length (); i++)
00183 elem (i, i) = a.elem (i, i);
00184 }
00185
00186 ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
00187 : MArray2<Complex> (1, rv.length (), 0.0)
00188 {
00189 for (int i = 0; i < rv.length (); i++)
00190 elem (0, i) = rv.elem (i);
00191 }
00192
00193 ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
00194 : MArray2<Complex> (cv.length (), 1, 0.0)
00195 {
00196 for (int i = 0; i < cv.length (); i++)
00197 elem (i, 0) = cv.elem (i);
00198 }
00199
00200 ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
00201 : MArray2<Complex> (a.rows (), a.cols (), 0.0)
00202 {
00203 for (int i = 0; i < a.length (); i++)
00204 elem (i, i) = a.elem (i, i);
00205 }
00206
00207
00208
00209
00210 ComplexMatrix::ComplexMatrix (const boolMatrix& a)
00211 : MArray2<Complex> (a.rows (), a.cols (), 0.0)
00212 {
00213 for (int i = 0; i < a.rows (); i++)
00214 for (int j = 0; j < a.cols (); j++)
00215 elem (i, j) = a.elem (i, j);
00216 }
00217
00218 ComplexMatrix::ComplexMatrix (const charMatrix& a)
00219 : MArray2<Complex> (a.rows (), a.cols (), 0.0)
00220 {
00221 for (int i = 0; i < a.rows (); i++)
00222 for (int j = 0; j < a.cols (); j++)
00223 elem (i, j) = a.elem (i, j);
00224 }
00225
00226 bool
00227 ComplexMatrix::operator == (const ComplexMatrix& a) const
00228 {
00229 if (rows () != a.rows () || cols () != a.cols ())
00230 return false;
00231
00232 return mx_inline_equal (data (), a.data (), length ());
00233 }
00234
00235 bool
00236 ComplexMatrix::operator != (const ComplexMatrix& a) const
00237 {
00238 return !(*this == a);
00239 }
00240
00241 bool
00242 ComplexMatrix::is_hermitian (void) const
00243 {
00244 int nr = rows ();
00245 int nc = cols ();
00246
00247 if (is_square () && nr > 0)
00248 {
00249 for (int i = 0; i < nr; i++)
00250 for (int j = i; j < nc; j++)
00251 if (elem (i, j) != conj (elem (j, i)))
00252 return false;
00253
00254 return true;
00255 }
00256
00257 return false;
00258 }
00259
00260
00261
00262 ComplexMatrix&
00263 ComplexMatrix::insert (const Matrix& a, int r, int c)
00264 {
00265 int a_nr = a.rows ();
00266 int a_nc = a.cols ();
00267
00268 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00269 {
00270 (*current_liboctave_error_handler) ("range error for insert");
00271 return *this;
00272 }
00273
00274 if (a_nr >0 && a_nc > 0)
00275 {
00276 make_unique ();
00277
00278 for (int j = 0; j < a_nc; j++)
00279 for (int i = 0; i < a_nr; i++)
00280 xelem (r+i, c+j) = a.elem (i, j);
00281 }
00282
00283 return *this;
00284 }
00285
00286 ComplexMatrix&
00287 ComplexMatrix::insert (const RowVector& a, int r, int c)
00288 {
00289 int a_len = a.length ();
00290
00291 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
00292 {
00293 (*current_liboctave_error_handler) ("range error for insert");
00294 return *this;
00295 }
00296
00297 if (a_len > 0)
00298 {
00299 make_unique ();
00300
00301 for (int i = 0; i < a_len; i++)
00302 xelem (r, c+i) = a.elem (i);
00303 }
00304
00305 return *this;
00306 }
00307
00308 ComplexMatrix&
00309 ComplexMatrix::insert (const ColumnVector& a, int r, int c)
00310 {
00311 int a_len = a.length ();
00312
00313 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
00314 {
00315 (*current_liboctave_error_handler) ("range error for insert");
00316 return *this;
00317 }
00318
00319 if (a_len > 0)
00320 {
00321 make_unique ();
00322
00323 for (int i = 0; i < a_len; i++)
00324 xelem (r+i, c) = a.elem (i);
00325 }
00326
00327 return *this;
00328 }
00329
00330 ComplexMatrix&
00331 ComplexMatrix::insert (const DiagMatrix& a, int r, int c)
00332 {
00333 int a_nr = a.rows ();
00334 int a_nc = a.cols ();
00335
00336 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00337 {
00338 (*current_liboctave_error_handler) ("range error for insert");
00339 return *this;
00340 }
00341
00342 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
00343
00344 int a_len = a.length ();
00345
00346 if (a_len > 0)
00347 {
00348 make_unique ();
00349
00350 for (int i = 0; i < a_len; i++)
00351 xelem (r+i, c+i) = a.elem (i, i);
00352 }
00353
00354 return *this;
00355 }
00356
00357 ComplexMatrix&
00358 ComplexMatrix::insert (const ComplexMatrix& a, int r, int c)
00359 {
00360 Array2<Complex>::insert (a, r, c);
00361 return *this;
00362 }
00363
00364 ComplexMatrix&
00365 ComplexMatrix::insert (const ComplexRowVector& a, int r, int c)
00366 {
00367 int a_len = a.length ();
00368 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
00369 {
00370 (*current_liboctave_error_handler) ("range error for insert");
00371 return *this;
00372 }
00373
00374 for (int i = 0; i < a_len; i++)
00375 elem (r, c+i) = a.elem (i);
00376
00377 return *this;
00378 }
00379
00380 ComplexMatrix&
00381 ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c)
00382 {
00383 int a_len = a.length ();
00384
00385 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
00386 {
00387 (*current_liboctave_error_handler) ("range error for insert");
00388 return *this;
00389 }
00390
00391 if (a_len > 0)
00392 {
00393 make_unique ();
00394
00395 for (int i = 0; i < a_len; i++)
00396 xelem (r+i, c) = a.elem (i);
00397 }
00398
00399 return *this;
00400 }
00401
00402 ComplexMatrix&
00403 ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c)
00404 {
00405 int a_nr = a.rows ();
00406 int a_nc = a.cols ();
00407
00408 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00409 {
00410 (*current_liboctave_error_handler) ("range error for insert");
00411 return *this;
00412 }
00413
00414 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
00415
00416 int a_len = a.length ();
00417
00418 if (a_len > 0)
00419 {
00420 make_unique ();
00421
00422 for (int i = 0; i < a_len; i++)
00423 xelem (r+i, c+i) = a.elem (i, i);
00424 }
00425
00426 return *this;
00427 }
00428
00429 ComplexMatrix&
00430 ComplexMatrix::fill (double val)
00431 {
00432 int nr = rows ();
00433 int nc = cols ();
00434
00435 if (nr > 0 && nc > 0)
00436 {
00437 make_unique ();
00438
00439 for (int j = 0; j < nc; j++)
00440 for (int i = 0; i < nr; i++)
00441 xelem (i, j) = val;
00442 }
00443
00444 return *this;
00445 }
00446
00447 ComplexMatrix&
00448 ComplexMatrix::fill (const Complex& val)
00449 {
00450 int nr = rows ();
00451 int nc = cols ();
00452
00453 if (nr > 0 && nc > 0)
00454 {
00455 make_unique ();
00456
00457 for (int j = 0; j < nc; j++)
00458 for (int i = 0; i < nr; i++)
00459 xelem (i, j) = val;
00460 }
00461
00462 return *this;
00463 }
00464
00465 ComplexMatrix&
00466 ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2)
00467 {
00468 int nr = rows ();
00469 int nc = cols ();
00470
00471 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
00472 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
00473 {
00474 (*current_liboctave_error_handler) ("range error for fill");
00475 return *this;
00476 }
00477
00478 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00479 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00480
00481 if (r2 >= r1 && c2 >= c1)
00482 {
00483 make_unique ();
00484
00485 for (int j = c1; j <= c2; j++)
00486 for (int i = r1; i <= r2; i++)
00487 xelem (i, j) = val;
00488 }
00489
00490 return *this;
00491 }
00492
00493 ComplexMatrix&
00494 ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2)
00495 {
00496 int nr = rows ();
00497 int nc = cols ();
00498
00499 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
00500 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
00501 {
00502 (*current_liboctave_error_handler) ("range error for fill");
00503 return *this;
00504 }
00505
00506 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00507 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00508
00509 if (r2 >= r1 && c2 >=c1)
00510 {
00511 make_unique ();
00512
00513 for (int j = c1; j <= c2; j++)
00514 for (int i = r1; i <= r2; i++)
00515 xelem (i, j) = val;
00516 }
00517
00518 return *this;
00519 }
00520
00521 ComplexMatrix
00522 ComplexMatrix::append (const Matrix& a) const
00523 {
00524 int nr = rows ();
00525 int nc = cols ();
00526 if (nr != a.rows ())
00527 {
00528 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00529 return *this;
00530 }
00531
00532 int nc_insert = nc;
00533 ComplexMatrix retval (nr, nc + a.cols ());
00534 retval.insert (*this, 0, 0);
00535 retval.insert (a, 0, nc_insert);
00536 return retval;
00537 }
00538
00539 ComplexMatrix
00540 ComplexMatrix::append (const RowVector& a) const
00541 {
00542 int nr = rows ();
00543 int nc = cols ();
00544 if (nr != 1)
00545 {
00546 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00547 return *this;
00548 }
00549
00550 int nc_insert = nc;
00551 ComplexMatrix retval (nr, nc + a.length ());
00552 retval.insert (*this, 0, 0);
00553 retval.insert (a, 0, nc_insert);
00554 return retval;
00555 }
00556
00557 ComplexMatrix
00558 ComplexMatrix::append (const ColumnVector& a) const
00559 {
00560 int nr = rows ();
00561 int nc = cols ();
00562 if (nr != a.length ())
00563 {
00564 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00565 return *this;
00566 }
00567
00568 int nc_insert = nc;
00569 ComplexMatrix retval (nr, nc + 1);
00570 retval.insert (*this, 0, 0);
00571 retval.insert (a, 0, nc_insert);
00572 return retval;
00573 }
00574
00575 ComplexMatrix
00576 ComplexMatrix::append (const DiagMatrix& a) const
00577 {
00578 int nr = rows ();
00579 int nc = cols ();
00580 if (nr != a.rows ())
00581 {
00582 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00583 return *this;
00584 }
00585
00586 int nc_insert = nc;
00587 ComplexMatrix retval (nr, nc + a.cols ());
00588 retval.insert (*this, 0, 0);
00589 retval.insert (a, 0, nc_insert);
00590 return retval;
00591 }
00592
00593 ComplexMatrix
00594 ComplexMatrix::append (const ComplexMatrix& a) const
00595 {
00596 int nr = rows ();
00597 int nc = cols ();
00598 if (nr != a.rows ())
00599 {
00600 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00601 return *this;
00602 }
00603
00604 int nc_insert = nc;
00605 ComplexMatrix retval (nr, nc + a.cols ());
00606 retval.insert (*this, 0, 0);
00607 retval.insert (a, 0, nc_insert);
00608 return retval;
00609 }
00610
00611 ComplexMatrix
00612 ComplexMatrix::append (const ComplexRowVector& a) const
00613 {
00614 int nr = rows ();
00615 int nc = cols ();
00616 if (nr != 1)
00617 {
00618 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00619 return *this;
00620 }
00621
00622 int nc_insert = nc;
00623 ComplexMatrix retval (nr, nc + a.length ());
00624 retval.insert (*this, 0, 0);
00625 retval.insert (a, 0, nc_insert);
00626 return retval;
00627 }
00628
00629 ComplexMatrix
00630 ComplexMatrix::append (const ComplexColumnVector& a) const
00631 {
00632 int nr = rows ();
00633 int nc = cols ();
00634 if (nr != a.length ())
00635 {
00636 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00637 return *this;
00638 }
00639
00640 int nc_insert = nc;
00641 ComplexMatrix retval (nr, nc + 1);
00642 retval.insert (*this, 0, 0);
00643 retval.insert (a, 0, nc_insert);
00644 return retval;
00645 }
00646
00647 ComplexMatrix
00648 ComplexMatrix::append (const ComplexDiagMatrix& a) const
00649 {
00650 int nr = rows ();
00651 int nc = cols ();
00652 if (nr != a.rows ())
00653 {
00654 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00655 return *this;
00656 }
00657
00658 int nc_insert = nc;
00659 ComplexMatrix retval (nr, nc + a.cols ());
00660 retval.insert (*this, 0, 0);
00661 retval.insert (a, 0, nc_insert);
00662 return retval;
00663 }
00664
00665 ComplexMatrix
00666 ComplexMatrix::stack (const Matrix& a) const
00667 {
00668 int nr = rows ();
00669 int nc = cols ();
00670 if (nc != a.cols ())
00671 {
00672 (*current_liboctave_error_handler)
00673 ("column dimension mismatch for stack");
00674 return *this;
00675 }
00676
00677 int nr_insert = nr;
00678 ComplexMatrix retval (nr + a.rows (), nc);
00679 retval.insert (*this, 0, 0);
00680 retval.insert (a, nr_insert, 0);
00681 return retval;
00682 }
00683
00684 ComplexMatrix
00685 ComplexMatrix::stack (const RowVector& a) const
00686 {
00687 int nr = rows ();
00688 int nc = cols ();
00689 if (nc != a.length ())
00690 {
00691 (*current_liboctave_error_handler)
00692 ("column dimension mismatch for stack");
00693 return *this;
00694 }
00695
00696 int nr_insert = nr;
00697 ComplexMatrix retval (nr + 1, nc);
00698 retval.insert (*this, 0, 0);
00699 retval.insert (a, nr_insert, 0);
00700 return retval;
00701 }
00702
00703 ComplexMatrix
00704 ComplexMatrix::stack (const ColumnVector& a) const
00705 {
00706 int nr = rows ();
00707 int nc = cols ();
00708 if (nc != 1)
00709 {
00710 (*current_liboctave_error_handler)
00711 ("column dimension mismatch for stack");
00712 return *this;
00713 }
00714
00715 int nr_insert = nr;
00716 ComplexMatrix retval (nr + a.length (), nc);
00717 retval.insert (*this, 0, 0);
00718 retval.insert (a, nr_insert, 0);
00719 return retval;
00720 }
00721
00722 ComplexMatrix
00723 ComplexMatrix::stack (const DiagMatrix& a) const
00724 {
00725 int nr = rows ();
00726 int nc = cols ();
00727 if (nc != a.cols ())
00728 {
00729 (*current_liboctave_error_handler)
00730 ("column dimension mismatch for stack");
00731 return *this;
00732 }
00733
00734 int nr_insert = nr;
00735 ComplexMatrix retval (nr + a.rows (), nc);
00736 retval.insert (*this, 0, 0);
00737 retval.insert (a, nr_insert, 0);
00738 return retval;
00739 }
00740
00741 ComplexMatrix
00742 ComplexMatrix::stack (const ComplexMatrix& a) const
00743 {
00744 int nr = rows ();
00745 int nc = cols ();
00746 if (nc != a.cols ())
00747 {
00748 (*current_liboctave_error_handler)
00749 ("column dimension mismatch for stack");
00750 return *this;
00751 }
00752
00753 int nr_insert = nr;
00754 ComplexMatrix retval (nr + a.rows (), nc);
00755 retval.insert (*this, 0, 0);
00756 retval.insert (a, nr_insert, 0);
00757 return retval;
00758 }
00759
00760 ComplexMatrix
00761 ComplexMatrix::stack (const ComplexRowVector& a) const
00762 {
00763 int nr = rows ();
00764 int nc = cols ();
00765 if (nc != a.length ())
00766 {
00767 (*current_liboctave_error_handler)
00768 ("column dimension mismatch for stack");
00769 return *this;
00770 }
00771
00772 int nr_insert = nr;
00773 ComplexMatrix retval (nr + 1, nc);
00774 retval.insert (*this, 0, 0);
00775 retval.insert (a, nr_insert, 0);
00776 return retval;
00777 }
00778
00779 ComplexMatrix
00780 ComplexMatrix::stack (const ComplexColumnVector& a) const
00781 {
00782 int nr = rows ();
00783 int nc = cols ();
00784 if (nc != 1)
00785 {
00786 (*current_liboctave_error_handler)
00787 ("column dimension mismatch for stack");
00788 return *this;
00789 }
00790
00791 int nr_insert = nr;
00792 ComplexMatrix retval (nr + a.length (), nc);
00793 retval.insert (*this, 0, 0);
00794 retval.insert (a, nr_insert, 0);
00795 return retval;
00796 }
00797
00798 ComplexMatrix
00799 ComplexMatrix::stack (const ComplexDiagMatrix& a) const
00800 {
00801 int nr = rows ();
00802 int nc = cols ();
00803 if (nc != a.cols ())
00804 {
00805 (*current_liboctave_error_handler)
00806 ("column dimension mismatch for stack");
00807 return *this;
00808 }
00809
00810 int nr_insert = nr;
00811 ComplexMatrix retval (nr + a.rows (), nc);
00812 retval.insert (*this, 0, 0);
00813 retval.insert (a, nr_insert, 0);
00814 return retval;
00815 }
00816
00817 ComplexMatrix
00818 ComplexMatrix::hermitian (void) const
00819 {
00820 int nr = rows ();
00821 int nc = cols ();
00822 ComplexMatrix result;
00823 if (length () > 0)
00824 {
00825 result.resize (nc, nr);
00826 for (int j = 0; j < nc; j++)
00827 for (int i = 0; i < nr; i++)
00828 result.elem (j, i) = conj (elem (i, j));
00829 }
00830 return result;
00831 }
00832
00833 ComplexMatrix
00834 conj (const ComplexMatrix& a)
00835 {
00836 int a_len = a.length ();
00837 ComplexMatrix retval;
00838 if (a_len > 0)
00839 retval = ComplexMatrix (mx_inline_conj_dup (a.data (), a_len),
00840 a.rows (), a.cols ());
00841 return retval;
00842 }
00843
00844
00845
00846 ComplexMatrix
00847 ComplexMatrix::extract (int r1, int c1, int r2, int c2) const
00848 {
00849 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00850 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00851
00852 int new_r = r2 - r1 + 1;
00853 int new_c = c2 - c1 + 1;
00854
00855 ComplexMatrix result (new_r, new_c);
00856
00857 for (int j = 0; j < new_c; j++)
00858 for (int i = 0; i < new_r; i++)
00859 result.xelem (i, j) = elem (r1+i, c1+j);
00860
00861 return result;
00862 }
00863
00864 ComplexMatrix
00865 ComplexMatrix::extract_n (int r1, int c1, int nr, int nc) const
00866 {
00867 ComplexMatrix result (nr, nc);
00868
00869 for (int j = 0; j < nc; j++)
00870 for (int i = 0; i < nr; i++)
00871 result.xelem (i, j) = elem (r1+i, c1+j);
00872
00873 return result;
00874 }
00875
00876
00877
00878 ComplexRowVector
00879 ComplexMatrix::row (int i) const
00880 {
00881 int nc = cols ();
00882 if (i < 0 || i >= rows ())
00883 {
00884 (*current_liboctave_error_handler) ("invalid row selection");
00885 return ComplexRowVector ();
00886 }
00887
00888 ComplexRowVector retval (nc);
00889 for (int j = 0; j < cols (); j++)
00890 retval.xelem (j) = elem (i, j);
00891
00892 return retval;
00893 }
00894
00895 ComplexRowVector
00896 ComplexMatrix::row (char *s) const
00897 {
00898 if (! s)
00899 {
00900 (*current_liboctave_error_handler) ("invalid row selection");
00901 return ComplexRowVector ();
00902 }
00903
00904 char c = *s;
00905 if (c == 'f' || c == 'F')
00906 return row (0);
00907 else if (c == 'l' || c == 'L')
00908 return row (rows () - 1);
00909 else
00910 {
00911 (*current_liboctave_error_handler) ("invalid row selection");
00912 return ComplexRowVector ();
00913 }
00914 }
00915
00916 ComplexColumnVector
00917 ComplexMatrix::column (int i) const
00918 {
00919 int nr = rows ();
00920 if (i < 0 || i >= cols ())
00921 {
00922 (*current_liboctave_error_handler) ("invalid column selection");
00923 return ComplexColumnVector ();
00924 }
00925
00926 ComplexColumnVector retval (nr);
00927 for (int j = 0; j < nr; j++)
00928 retval.xelem (j) = elem (j, i);
00929
00930 return retval;
00931 }
00932
00933 ComplexColumnVector
00934 ComplexMatrix::column (char *s) const
00935 {
00936 if (! s)
00937 {
00938 (*current_liboctave_error_handler) ("invalid column selection");
00939 return ComplexColumnVector ();
00940 }
00941
00942 char c = *s;
00943 if (c == 'f' || c == 'F')
00944 return column (0);
00945 else if (c == 'l' || c == 'L')
00946 return column (cols () - 1);
00947 else
00948 {
00949 (*current_liboctave_error_handler) ("invalid column selection");
00950 return ComplexColumnVector ();
00951 }
00952 }
00953
00954 ComplexMatrix
00955 ComplexMatrix::inverse (void) const
00956 {
00957 int info;
00958 double rcond;
00959 return inverse (info, rcond, 0, 0);
00960 }
00961
00962 ComplexMatrix
00963 ComplexMatrix::inverse (int& info) const
00964 {
00965 double rcond;
00966 return inverse (info, rcond, 0, 0);
00967 }
00968
00969 ComplexMatrix
00970 ComplexMatrix::inverse (int& info, double& rcond, int force,
00971 int calc_cond) const
00972 {
00973 ComplexMatrix retval;
00974
00975 int nr = rows ();
00976 int nc = cols ();
00977
00978 if (nr != nc)
00979 (*current_liboctave_error_handler) ("inverse requires square matrix");
00980 else
00981 {
00982 Array<int> ipvt (nr);
00983 int *pipvt = ipvt.fortran_vec ();
00984
00985 retval = *this;
00986 Complex *tmp_data = retval.fortran_vec ();
00987
00988 Array<Complex> z(1);
00989 int lwork = -1;
00990
00991
00992
00993 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
00994 z.fortran_vec (), lwork, info));
00995
00996 if (f77_exception_encountered)
00997 {
00998 (*current_liboctave_error_handler)
00999 ("unrecoverable error in zgetri");
01000 return retval;
01001 }
01002
01003 lwork = static_cast<int> (real(z(0)));
01004 lwork = (lwork < 2 *nc ? 2*nc : lwork);
01005 z.resize (lwork);
01006 Complex *pz = z.fortran_vec ();
01007
01008 info = 0;
01009
01010
01011 double anorm;
01012 if (calc_cond)
01013 anorm = retval.abs().sum().row(0).max();
01014
01015 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info));
01016
01017 if (f77_exception_encountered)
01018 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
01019 else
01020 {
01021
01022 rcond = 0.0;
01023 if (info != 0)
01024 info = -1;
01025 else if (calc_cond)
01026 {
01027
01028 int zgecon_info = 0;
01029 char job = '1';
01030 Array<double> rz (2 * nc);
01031 double *prz = rz.fortran_vec ();
01032 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01033 nc, tmp_data, nr, anorm,
01034 rcond, pz, prz, zgecon_info
01035 F77_CHAR_ARG_LEN (1)));
01036
01037 if (f77_exception_encountered)
01038 (*current_liboctave_error_handler)
01039 ("unrecoverable error in zgecon");
01040
01041 if (zgecon_info != 0)
01042 info = -1;
01043 }
01044
01045 if (info == -1 && ! force)
01046 retval = *this;
01047 else
01048 {
01049 int zgetri_info = 0;
01050
01051 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
01052 pz, lwork, zgetri_info));
01053
01054 if (f77_exception_encountered)
01055 (*current_liboctave_error_handler)
01056 ("unrecoverable error in zgetri");
01057
01058 if (zgetri_info != 0)
01059 info = -1;
01060 }
01061 }
01062 }
01063
01064 return retval;
01065 }
01066
01067 ComplexMatrix
01068 ComplexMatrix::pseudo_inverse (double tol) const
01069 {
01070 ComplexMatrix retval;
01071
01072 ComplexSVD result (*this, SVD::economy);
01073
01074 DiagMatrix S = result.singular_values ();
01075 ComplexMatrix U = result.left_singular_matrix ();
01076 ComplexMatrix V = result.right_singular_matrix ();
01077
01078 ColumnVector sigma = S.diag ();
01079
01080 int r = sigma.length () - 1;
01081 int nr = rows ();
01082 int nc = cols ();
01083
01084 if (tol <= 0.0)
01085 {
01086 if (nr > nc)
01087 tol = nr * sigma.elem (0) * DBL_EPSILON;
01088 else
01089 tol = nc * sigma.elem (0) * DBL_EPSILON;
01090 }
01091
01092 while (r >= 0 && sigma.elem (r) < tol)
01093 r--;
01094
01095 if (r < 0)
01096 retval = ComplexMatrix (nc, nr, 0.0);
01097 else
01098 {
01099 ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
01100 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
01101 ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
01102 retval = Vr * D * Ur.hermitian ();
01103 }
01104
01105 return retval;
01106 }
01107
01108 #if defined (HAVE_FFTW3)
01109
01110 ComplexMatrix
01111 ComplexMatrix::fourier (void) const
01112 {
01113 size_t nr = rows ();
01114 size_t nc = cols ();
01115
01116 ComplexMatrix retval (nr, nc);
01117
01118 size_t npts, nsamples;
01119
01120 if (nr == 1 || nc == 1)
01121 {
01122 npts = nr > nc ? nr : nc;
01123 nsamples = 1;
01124 }
01125 else
01126 {
01127 npts = nr;
01128 nsamples = nc;
01129 }
01130
01131 const Complex *in (data ());
01132 Complex *out (retval.fortran_vec ());
01133
01134 octave_fftw::fft (in, out, npts, nsamples);
01135
01136 return retval;
01137 }
01138
01139 ComplexMatrix
01140 ComplexMatrix::ifourier (void) const
01141 {
01142 size_t nr = rows ();
01143 size_t nc = cols ();
01144
01145 ComplexMatrix retval (nr, nc);
01146
01147 size_t npts, nsamples;
01148
01149 if (nr == 1 || nc == 1)
01150 {
01151 npts = nr > nc ? nr : nc;
01152 nsamples = 1;
01153 }
01154 else
01155 {
01156 npts = nr;
01157 nsamples = nc;
01158 }
01159
01160 const Complex *in (data ());
01161 Complex *out (retval.fortran_vec ());
01162
01163 octave_fftw::ifft (in, out, npts, nsamples);
01164
01165 return retval;
01166 }
01167
01168 ComplexMatrix
01169 ComplexMatrix::fourier2d (void) const
01170 {
01171 dim_vector dv(rows (), cols ());
01172
01173 ComplexMatrix retval (rows (), cols ());
01174 const Complex *in (data ());
01175 Complex *out (retval.fortran_vec ());
01176
01177 octave_fftw::fftNd (in, out, 2, dv);
01178
01179 return retval;
01180 }
01181
01182 ComplexMatrix
01183 ComplexMatrix::ifourier2d (void) const
01184 {
01185 dim_vector dv(rows (), cols ());
01186
01187 ComplexMatrix retval (rows (), cols ());
01188 const Complex *in (data ());
01189 Complex *out (retval.fortran_vec ());
01190
01191 octave_fftw::ifftNd (in, out, 2, dv);
01192
01193 return retval;
01194 }
01195
01196 #else
01197
01198 ComplexMatrix
01199 ComplexMatrix::fourier (void) const
01200 {
01201 ComplexMatrix retval;
01202
01203 int nr = rows ();
01204 int nc = cols ();
01205
01206 int npts, nsamples;
01207
01208 if (nr == 1 || nc == 1)
01209 {
01210 npts = nr > nc ? nr : nc;
01211 nsamples = 1;
01212 }
01213 else
01214 {
01215 npts = nr;
01216 nsamples = nc;
01217 }
01218
01219 int nn = 4*npts+15;
01220
01221 Array<Complex> wsave (nn);
01222 Complex *pwsave = wsave.fortran_vec ();
01223
01224 retval = *this;
01225 Complex *tmp_data = retval.fortran_vec ();
01226
01227 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01228
01229 for (int j = 0; j < nsamples; j++)
01230 {
01231 OCTAVE_QUIT;
01232
01233 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
01234 }
01235
01236 return retval;
01237 }
01238
01239 ComplexMatrix
01240 ComplexMatrix::ifourier (void) const
01241 {
01242 ComplexMatrix retval;
01243
01244 int nr = rows ();
01245 int nc = cols ();
01246
01247 int npts, nsamples;
01248
01249 if (nr == 1 || nc == 1)
01250 {
01251 npts = nr > nc ? nr : nc;
01252 nsamples = 1;
01253 }
01254 else
01255 {
01256 npts = nr;
01257 nsamples = nc;
01258 }
01259
01260 int nn = 4*npts+15;
01261
01262 Array<Complex> wsave (nn);
01263 Complex *pwsave = wsave.fortran_vec ();
01264
01265 retval = *this;
01266 Complex *tmp_data = retval.fortran_vec ();
01267
01268 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01269
01270 for (int j = 0; j < nsamples; j++)
01271 {
01272 OCTAVE_QUIT;
01273
01274 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
01275 }
01276
01277 for (int j = 0; j < npts*nsamples; j++)
01278 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
01279
01280 return retval;
01281 }
01282
01283 ComplexMatrix
01284 ComplexMatrix::fourier2d (void) const
01285 {
01286 ComplexMatrix retval;
01287
01288 int nr = rows ();
01289 int nc = cols ();
01290
01291 int npts, nsamples;
01292
01293 if (nr == 1 || nc == 1)
01294 {
01295 npts = nr > nc ? nr : nc;
01296 nsamples = 1;
01297 }
01298 else
01299 {
01300 npts = nr;
01301 nsamples = nc;
01302 }
01303
01304 int nn = 4*npts+15;
01305
01306 Array<Complex> wsave (nn);
01307 Complex *pwsave = wsave.fortran_vec ();
01308
01309 retval = *this;
01310 Complex *tmp_data = retval.fortran_vec ();
01311
01312 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01313
01314 for (int j = 0; j < nsamples; j++)
01315 {
01316 OCTAVE_QUIT;
01317
01318 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
01319 }
01320
01321 npts = nc;
01322 nsamples = nr;
01323 nn = 4*npts+15;
01324
01325 wsave.resize (nn);
01326 pwsave = wsave.fortran_vec ();
01327
01328 Array<Complex> tmp (npts);
01329 Complex *prow = tmp.fortran_vec ();
01330
01331 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01332
01333 for (int j = 0; j < nsamples; j++)
01334 {
01335 OCTAVE_QUIT;
01336
01337 for (int i = 0; i < npts; i++)
01338 prow[i] = tmp_data[i*nr + j];
01339
01340 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
01341
01342 for (int i = 0; i < npts; i++)
01343 tmp_data[i*nr + j] = prow[i];
01344 }
01345
01346 return retval;
01347 }
01348
01349 ComplexMatrix
01350 ComplexMatrix::ifourier2d (void) const
01351 {
01352 ComplexMatrix retval;
01353
01354 int nr = rows ();
01355 int nc = cols ();
01356
01357 int npts, nsamples;
01358
01359 if (nr == 1 || nc == 1)
01360 {
01361 npts = nr > nc ? nr : nc;
01362 nsamples = 1;
01363 }
01364 else
01365 {
01366 npts = nr;
01367 nsamples = nc;
01368 }
01369
01370 int nn = 4*npts+15;
01371
01372 Array<Complex> wsave (nn);
01373 Complex *pwsave = wsave.fortran_vec ();
01374
01375 retval = *this;
01376 Complex *tmp_data = retval.fortran_vec ();
01377
01378 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01379
01380 for (int j = 0; j < nsamples; j++)
01381 {
01382 OCTAVE_QUIT;
01383
01384 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
01385 }
01386
01387 for (int j = 0; j < npts*nsamples; j++)
01388 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
01389
01390 npts = nc;
01391 nsamples = nr;
01392 nn = 4*npts+15;
01393
01394 wsave.resize (nn);
01395 pwsave = wsave.fortran_vec ();
01396
01397 Array<Complex> tmp (npts);
01398 Complex *prow = tmp.fortran_vec ();
01399
01400 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01401
01402 for (int j = 0; j < nsamples; j++)
01403 {
01404 OCTAVE_QUIT;
01405
01406 for (int i = 0; i < npts; i++)
01407 prow[i] = tmp_data[i*nr + j];
01408
01409 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
01410
01411 for (int i = 0; i < npts; i++)
01412 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts);
01413 }
01414
01415 return retval;
01416 }
01417
01418 #endif
01419
01420 ComplexDET
01421 ComplexMatrix::determinant (void) const
01422 {
01423 int info;
01424 double rcond;
01425 return determinant (info, rcond, 0);
01426 }
01427
01428 ComplexDET
01429 ComplexMatrix::determinant (int& info) const
01430 {
01431 double rcond;
01432 return determinant (info, rcond, 0);
01433 }
01434
01435 ComplexDET
01436 ComplexMatrix::determinant (int& info, double& rcond, int calc_cond) const
01437 {
01438 ComplexDET retval;
01439
01440 int nr = rows ();
01441 int nc = cols ();
01442
01443 if (nr == 0 || nc == 0)
01444 {
01445 Complex d[2];
01446 d[0] = 1.0;
01447 d[1] = 0.0;
01448 retval = ComplexDET (d);
01449 }
01450 else
01451 {
01452 Array<int> ipvt (nr);
01453 int *pipvt = ipvt.fortran_vec ();
01454
01455 ComplexMatrix atmp = *this;
01456 Complex *tmp_data = atmp.fortran_vec ();
01457
01458 info = 0;
01459
01460
01461 double anorm = 0;
01462 if (calc_cond)
01463 anorm = atmp.abs().sum().row(0).max();
01464
01465 F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info));
01466
01467 if (f77_exception_encountered)
01468 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
01469 else
01470 {
01471
01472 rcond = 0.0;
01473 if (info != 0)
01474 {
01475 info = -1;
01476 retval = ComplexDET ();
01477 }
01478 else
01479 {
01480 if (calc_cond)
01481 {
01482
01483 char job = '1';
01484 Array<Complex> z (2*nr);
01485 Complex *pz = z.fortran_vec ();
01486 Array<double> rz (2*nr);
01487 double *prz = rz.fortran_vec ();
01488
01489 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01490 nc, tmp_data, nr, anorm,
01491 rcond, pz, prz, info
01492 F77_CHAR_ARG_LEN (1)));
01493
01494 if (f77_exception_encountered)
01495 (*current_liboctave_error_handler)
01496 ("unrecoverable error in zgecon");
01497 }
01498
01499 if (info != 0)
01500 {
01501 info = -1;
01502 retval = ComplexDET ();
01503 }
01504 else
01505 {
01506 Complex d[2] = { 1., 0.};
01507 for (int i=0; i<nc; i++)
01508 {
01509 if (ipvt(i) != (i+1)) d[0] = -d[0];
01510 d[0] = d[0] * atmp(i,i);
01511 if (d[0] == 0.) break;
01512 while (::abs(d[0]) < 1.)
01513 {
01514 d[0] = 10. * d[0];
01515 d[1] = d[1] - 1.0;
01516 }
01517 while (::abs(d[0]) >= 10.)
01518 {
01519 d[0] = 0.1 * d[0];
01520 d[1] = d[1] + 1.0;
01521 }
01522 }
01523 retval = ComplexDET (d);
01524 }
01525 }
01526 }
01527 }
01528
01529 return retval;
01530 }
01531
01532 ComplexMatrix
01533 ComplexMatrix::solve (const Matrix& b) const
01534 {
01535 int info;
01536 double rcond;
01537 return solve (b, info, rcond, 0);
01538 }
01539
01540 ComplexMatrix
01541 ComplexMatrix::solve (const Matrix& b, int& info) const
01542 {
01543 double rcond;
01544 return solve (b, info, rcond, 0);
01545 }
01546
01547 ComplexMatrix
01548 ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const
01549 {
01550 return solve (b, info, rcond, 0);
01551 }
01552
01553 ComplexMatrix
01554 ComplexMatrix::solve (const Matrix& b, int& info, double& rcond,
01555 solve_singularity_handler sing_handler) const
01556 {
01557 ComplexMatrix tmp (b);
01558 return solve (tmp, info, rcond, sing_handler);
01559 }
01560
01561 ComplexMatrix
01562 ComplexMatrix::solve (const ComplexMatrix& b) const
01563 {
01564 int info;
01565 double rcond;
01566 return solve (b, info, rcond, 0);
01567 }
01568
01569 ComplexMatrix
01570 ComplexMatrix::solve (const ComplexMatrix& b, int& info) const
01571 {
01572 double rcond;
01573 return solve (b, info, rcond, 0);
01574 }
01575
01576 ComplexMatrix
01577 ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
01578 {
01579 return solve (b, info, rcond, 0);
01580 }
01581
01582 ComplexMatrix
01583 ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond,
01584 solve_singularity_handler sing_handler) const
01585 {
01586 ComplexMatrix retval;
01587
01588 int nr = rows ();
01589 int nc = cols ();
01590
01591 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
01592 (*current_liboctave_error_handler)
01593 ("matrix dimension mismatch in solution of linear equations");
01594 else
01595 {
01596 info = 0;
01597
01598 Array<int> ipvt (nr);
01599 int *pipvt = ipvt.fortran_vec ();
01600
01601 ComplexMatrix atmp = *this;
01602 Complex *tmp_data = atmp.fortran_vec ();
01603
01604 Array<Complex> z (2 * nc);
01605 Complex *pz = z.fortran_vec ();
01606 Array<double> rz (2 * nc);
01607 double *prz = rz.fortran_vec ();
01608
01609
01610 double anorm = atmp.abs().sum().row(0).max();
01611
01612 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01613
01614 if (f77_exception_encountered)
01615 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
01616 else
01617 {
01618
01619 rcond = 0.0;
01620 if (info != 0)
01621 {
01622 info = -2;
01623
01624 if (sing_handler)
01625 sing_handler (rcond);
01626 else
01627 (*current_liboctave_error_handler)
01628 ("matrix singular to machine precision");
01629
01630 }
01631 else
01632 {
01633
01634 char job = '1';
01635 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01636 nc, tmp_data, nr, anorm,
01637 rcond, pz, prz, info
01638 F77_CHAR_ARG_LEN (1)));
01639
01640 if (f77_exception_encountered)
01641 (*current_liboctave_error_handler)
01642 ("unrecoverable error in zgecon");
01643
01644 if (info != 0)
01645 info = -2;
01646
01647 volatile double rcond_plus_one = rcond + 1.0;
01648
01649 if (rcond_plus_one == 1.0 || xisnan (rcond))
01650 {
01651 info = -2;
01652
01653 if (sing_handler)
01654 sing_handler (rcond);
01655 else
01656 (*current_liboctave_error_handler)
01657 ("matrix singular to machine precision, rcond = %g",
01658 rcond);
01659 }
01660 else
01661 {
01662 retval = b;
01663 Complex *result = retval.fortran_vec ();
01664
01665 int b_nc = b.cols ();
01666
01667 job = 'N';
01668 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
01669 nr, b_nc, tmp_data, nr,
01670 pipvt, result, b.rows(), info
01671 F77_CHAR_ARG_LEN (1)));
01672
01673 if (f77_exception_encountered)
01674 (*current_liboctave_error_handler)
01675 ("unrecoverable error in zgetrs");
01676 }
01677 }
01678 }
01679 }
01680
01681 return retval;
01682 }
01683
01684 ComplexColumnVector
01685 ComplexMatrix::solve (const ColumnVector& b) const
01686 {
01687 int info;
01688 double rcond;
01689 return solve (ComplexColumnVector (b), info, rcond, 0);
01690 }
01691
01692 ComplexColumnVector
01693 ComplexMatrix::solve (const ColumnVector& b, int& info) const
01694 {
01695 double rcond;
01696 return solve (ComplexColumnVector (b), info, rcond, 0);
01697 }
01698
01699 ComplexColumnVector
01700 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond) const
01701 {
01702 return solve (ComplexColumnVector (b), info, rcond, 0);
01703 }
01704
01705 ComplexColumnVector
01706 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond,
01707 solve_singularity_handler sing_handler) const
01708 {
01709 return solve (ComplexColumnVector (b), info, rcond, sing_handler);
01710 }
01711
01712 ComplexColumnVector
01713 ComplexMatrix::solve (const ComplexColumnVector& b) const
01714 {
01715 int info;
01716 double rcond;
01717 return solve (b, info, rcond, 0);
01718 }
01719
01720 ComplexColumnVector
01721 ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
01722 {
01723 double rcond;
01724 return solve (b, info, rcond, 0);
01725 }
01726
01727 ComplexColumnVector
01728 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
01729 double& rcond) const
01730 {
01731 return solve (b, info, rcond, 0);
01732 }
01733
01734 ComplexColumnVector
01735 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
01736 double& rcond,
01737 solve_singularity_handler sing_handler) const
01738 {
01739 ComplexColumnVector retval;
01740
01741 int nr = rows ();
01742 int nc = cols ();
01743
01744 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
01745 (*current_liboctave_error_handler)
01746 ("matrix dimension mismatch in solution of linear equations");
01747 else
01748 {
01749 info = 0;
01750
01751 Array<int> ipvt (nr);
01752 int *pipvt = ipvt.fortran_vec ();
01753
01754 ComplexMatrix atmp = *this;
01755 Complex *tmp_data = atmp.fortran_vec ();
01756
01757 Array<Complex> z (2 * nc);
01758 Complex *pz = z.fortran_vec ();
01759 Array<double> rz (2 * nc);
01760 double *prz = rz.fortran_vec ();
01761
01762
01763 double anorm = atmp.abs().sum().row(0).max();
01764
01765 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01766
01767 if (f77_exception_encountered)
01768 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
01769 else
01770 {
01771
01772 rcond = 0.0;
01773 if (info != 0)
01774 {
01775 info = -2;
01776
01777 if (sing_handler)
01778 sing_handler (rcond);
01779 else
01780 (*current_liboctave_error_handler)
01781 ("matrix singular to machine precision, rcond = %g",
01782 rcond);
01783 }
01784 else
01785 {
01786
01787 char job = '1';
01788 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01789 nc, tmp_data, nr, anorm,
01790 rcond, pz, prz, info
01791 F77_CHAR_ARG_LEN (1)));
01792
01793 if (f77_exception_encountered)
01794 (*current_liboctave_error_handler)
01795 ("unrecoverable error in zgecon");
01796
01797 if (info != 0)
01798 info = -2;
01799
01800 volatile double rcond_plus_one = rcond + 1.0;
01801
01802 if (rcond_plus_one == 1.0 || xisnan (rcond))
01803 {
01804 info = -2;
01805
01806 if (sing_handler)
01807 sing_handler (rcond);
01808 else
01809 (*current_liboctave_error_handler)
01810 ("matrix singular to machine precision, rcond = %g",
01811 rcond);
01812 }
01813 else
01814 {
01815 retval = b;
01816 Complex *result = retval.fortran_vec ();
01817
01818 job = 'N';
01819 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
01820 nr, 1, tmp_data, nr, pipvt,
01821 result, b.length(), info
01822 F77_CHAR_ARG_LEN (1)));
01823
01824 if (f77_exception_encountered)
01825 (*current_liboctave_error_handler)
01826 ("unrecoverable error in zgetrs");
01827
01828 }
01829 }
01830 }
01831 }
01832 return retval;
01833 }
01834
01835 ComplexMatrix
01836 ComplexMatrix::lssolve (const Matrix& b) const
01837 {
01838 int info;
01839 int rank;
01840 return lssolve (ComplexMatrix (b), info, rank);
01841 }
01842
01843 ComplexMatrix
01844 ComplexMatrix::lssolve (const Matrix& b, int& info) const
01845 {
01846 int rank;
01847 return lssolve (ComplexMatrix (b), info, rank);
01848 }
01849
01850 ComplexMatrix
01851 ComplexMatrix::lssolve (const Matrix& b, int& info, int& rank) const
01852 {
01853 return lssolve (ComplexMatrix (b), info, rank);
01854 }
01855
01856 ComplexMatrix
01857 ComplexMatrix::lssolve (const ComplexMatrix& b) const
01858 {
01859 int info;
01860 int rank;
01861 return lssolve (b, info, rank);
01862 }
01863
01864 ComplexMatrix
01865 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const
01866 {
01867 int rank;
01868 return lssolve (b, info, rank);
01869 }
01870
01871 ComplexMatrix
01872 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
01873 {
01874 ComplexMatrix retval;
01875
01876 int nrhs = b.cols ();
01877
01878 int m = rows ();
01879 int n = cols ();
01880
01881 if (m == 0 || n == 0 || m != b.rows ())
01882 (*current_liboctave_error_handler)
01883 ("matrix dimension mismatch solution of linear equations");
01884 else
01885 {
01886 ComplexMatrix atmp = *this;
01887 Complex *tmp_data = atmp.fortran_vec ();
01888
01889 int nrr = m > n ? m : n;
01890 ComplexMatrix result (nrr, nrhs);
01891
01892 for (int j = 0; j < nrhs; j++)
01893 for (int i = 0; i < m; i++)
01894 result.elem (i, j) = b.elem (i, j);
01895
01896 Complex *presult = result.fortran_vec ();
01897
01898 int len_s = m < n ? m : n;
01899 Array<double> s (len_s);
01900 double *ps = s.fortran_vec ();
01901
01902 double rcond = -1.0;
01903
01904 int lrwork = (5 * (m < n ? m : n)) - 4;
01905 lrwork = lrwork > 1 ? lrwork : 1;
01906 Array<double> rwork (lrwork);
01907 double *prwork = rwork.fortran_vec ();
01908
01909
01910
01911 int lwork = -1;
01912
01913 Array<Complex> work (1);
01914
01915 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
01916 nrr, ps, rcond, rank,
01917 work.fortran_vec (), lwork, prwork,
01918 info));
01919
01920 if (f77_exception_encountered)
01921 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
01922 else
01923 {
01924 lwork = static_cast<int> (real (work(0)));
01925 work.resize (lwork);
01926
01927 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
01928 nrr, ps, rcond, rank,
01929 work.fortran_vec (), lwork,
01930 prwork, info));
01931
01932 if (f77_exception_encountered)
01933 (*current_liboctave_error_handler)
01934 ("unrecoverable error in zgelss");
01935 else
01936 {
01937 retval.resize (n, nrhs);
01938 for (int j = 0; j < nrhs; j++)
01939 for (int i = 0; i < n; i++)
01940 retval.elem (i, j) = result.elem (i, j);
01941 }
01942 }
01943 }
01944
01945 return retval;
01946 }
01947
01948 ComplexColumnVector
01949 ComplexMatrix::lssolve (const ColumnVector& b) const
01950 {
01951 int info;
01952 int rank;
01953 return lssolve (ComplexColumnVector (b), info, rank);
01954 }
01955
01956 ComplexColumnVector
01957 ComplexMatrix::lssolve (const ColumnVector& b, int& info) const
01958 {
01959 int rank;
01960 return lssolve (ComplexColumnVector (b), info, rank);
01961 }
01962
01963 ComplexColumnVector
01964 ComplexMatrix::lssolve (const ColumnVector& b, int& info, int& rank) const
01965 {
01966 return lssolve (ComplexColumnVector (b), info, rank);
01967 }
01968
01969 ComplexColumnVector
01970 ComplexMatrix::lssolve (const ComplexColumnVector& b) const
01971 {
01972 int info;
01973 int rank;
01974 return lssolve (b, info, rank);
01975 }
01976
01977 ComplexColumnVector
01978 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const
01979 {
01980 int rank;
01981 return lssolve (b, info, rank);
01982 }
01983
01984 ComplexColumnVector
01985 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
01986 int& rank) const
01987 {
01988 ComplexColumnVector retval;
01989
01990 int nrhs = 1;
01991
01992 int m = rows ();
01993 int n = cols ();
01994
01995 if (m == 0 || n == 0 || m != b.length ())
01996 (*current_liboctave_error_handler)
01997 ("matrix dimension mismatch solution of least squares problem");
01998 else
01999 {
02000 ComplexMatrix atmp = *this;
02001 Complex *tmp_data = atmp.fortran_vec ();
02002
02003 int nrr = m > n ? m : n;
02004 ComplexColumnVector result (nrr);
02005
02006 for (int i = 0; i < m; i++)
02007 result.elem (i) = b.elem (i);
02008
02009 Complex *presult = result.fortran_vec ();
02010
02011 int len_s = m < n ? m : n;
02012 Array<double> s (len_s);
02013 double *ps = s.fortran_vec ();
02014
02015 double rcond = -1.0;
02016
02017 int lrwork = (5 * (m < n ? m : n)) - 4;
02018 lrwork = lrwork > 1 ? lrwork : 1;
02019 Array<double> rwork (lrwork);
02020 double *prwork = rwork.fortran_vec ();
02021
02022
02023
02024 int lwork = -1;
02025
02026 Array<Complex> work (1);
02027
02028 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
02029 nrr, ps, rcond, rank,
02030 work.fortran_vec (), lwork, prwork,
02031 info));
02032
02033 if (f77_exception_encountered)
02034 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
02035 else
02036 {
02037 lwork = static_cast<int> (real (work(0)));
02038 work.resize (lwork);
02039
02040 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
02041 nrr, ps, rcond, rank,
02042 work.fortran_vec (), lwork,
02043 prwork, info));
02044
02045 if (f77_exception_encountered)
02046 (*current_liboctave_error_handler)
02047 ("unrecoverable error in zgelss");
02048 else
02049 {
02050 retval.resize (n);
02051 for (int i = 0; i < n; i++)
02052 retval.elem (i) = result.elem (i);
02053 }
02054 }
02055 }
02056
02057 return retval;
02058 }
02059
02060
02061
02062 static double padec [] =
02063 {
02064 5.0000000000000000e-1,
02065 1.1666666666666667e-1,
02066 1.6666666666666667e-2,
02067 1.6025641025641026e-3,
02068 1.0683760683760684e-4,
02069 4.8562548562548563e-6,
02070 1.3875013875013875e-7,
02071 1.9270852604185938e-9,
02072 };
02073
02074 ComplexMatrix
02075 ComplexMatrix::expm (void) const
02076 {
02077 ComplexMatrix retval;
02078
02079 ComplexMatrix m = *this;
02080
02081 int nc = columns ();
02082
02083
02084
02085
02086
02087 Complex trshift = 0.0;
02088
02089 for (int i = 0; i < nc; i++)
02090 trshift += m.elem (i, i);
02091
02092 trshift /= nc;
02093
02094 if (trshift.real () < 0.0)
02095 trshift = trshift.imag ();
02096
02097 for (int i = 0; i < nc; i++)
02098 m.elem (i, i) -= trshift;
02099
02100
02101
02102
02103 Complex *mp = m.fortran_vec ();
02104
02105 int info, ilo, ihi,ilos,ihis;
02106 Array<double> dpermute (nc);
02107 Array<double> dscale (nc);
02108
02109
02110
02111
02112 char job = 'P';
02113 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
02114 nc, mp, nc, ilo, ihi,
02115 dpermute.fortran_vec (), info
02116 F77_CHAR_ARG_LEN (1)));
02117
02118 if (f77_exception_encountered)
02119 {
02120 (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
02121 return retval;
02122 }
02123
02124
02125 job = 'S';
02126 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
02127 nc, mp, nc, ilos, ihis,
02128 dscale.fortran_vec (), info
02129 F77_CHAR_ARG_LEN (1)));
02130
02131 if (f77_exception_encountered)
02132 {
02133 (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
02134 return retval;
02135 }
02136
02137
02138
02139 ColumnVector work (nc);
02140 double inf_norm;
02141
02142 F77_XFCN (xzlange, XZLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
02143 nc, nc, m.fortran_vec (), nc,
02144 work.fortran_vec (), inf_norm
02145 F77_CHAR_ARG_LEN (1)));
02146
02147 if (f77_exception_encountered)
02148 {
02149 (*current_liboctave_error_handler) ("unrecoverable error in zlange");
02150 return retval;
02151 }
02152
02153 int sqpow = (inf_norm > 0.0
02154 ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
02155
02156
02157
02158 if (sqpow < 0)
02159 sqpow = 0;
02160
02161 if (sqpow > 0)
02162 {
02163 double scale_factor = 1.0;
02164 for (int i = 0; i < sqpow; i++)
02165 scale_factor *= 2.0;
02166
02167 m = m / scale_factor;
02168 }
02169
02170
02171
02172 ComplexMatrix npp (nc, nc, 0.0);
02173 ComplexMatrix dpp = npp;
02174
02175
02176
02177 int minus_one_j = -1;
02178 for (int j = 7; j >= 0; j--)
02179 {
02180 npp = m * npp + m * padec[j];
02181 dpp = m * dpp + m * (minus_one_j * padec[j]);
02182 minus_one_j *= -1;
02183 }
02184
02185
02186
02187 dpp = -dpp;
02188 for (int j = 0; j < nc; j++)
02189 {
02190 npp.elem (j, j) += 1.0;
02191 dpp.elem (j, j) += 1.0;
02192 }
02193
02194
02195
02196 retval = dpp.solve (npp);
02197
02198
02199
02200 while (sqpow)
02201 {
02202 retval = retval * retval;
02203 sqpow--;
02204 }
02205
02206
02207
02208
02209
02210 for (int i = 0; i < nc; i++)
02211 for (int j = 0; j < nc; j++)
02212 retval(i,j) *= dscale(i) / dscale(j);
02213
02214 OCTAVE_QUIT;
02215
02216
02217 Array<int> iperm (nc);
02218 for (int i = 0; i < nc; i++)
02219 iperm(i) = i;
02220
02221
02222 for (int i = 0; i < (ilo-1); i++)
02223 {
02224 int swapidx = static_cast<int> (dpermute(i)) - 1;
02225 int tmp = iperm(i);
02226 iperm(i) = iperm(swapidx);
02227 iperm(swapidx) = tmp;
02228 }
02229
02230
02231 for (int i = nc - 1; i >= ihi; i--)
02232 {
02233 int swapidx = static_cast<int> (dpermute(i)) - 1;
02234 int tmp = iperm(i);
02235 iperm(i) = iperm(swapidx);
02236 iperm(swapidx) = tmp;
02237 }
02238
02239
02240 Array<int> invpvec (nc);
02241 for (int i = 0; i < nc; i++)
02242 invpvec(iperm(i)) = i;
02243
02244 OCTAVE_QUIT;
02245
02246 ComplexMatrix tmpMat = retval;
02247 for (int i = 0; i < nc; i++)
02248 for (int j = 0; j < nc; j++)
02249 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
02250
02251
02252
02253 return exp (trshift) * retval;
02254 }
02255
02256
02257
02258 ComplexMatrix
02259 operator * (const ColumnVector& v, const ComplexRowVector& a)
02260 {
02261 ComplexColumnVector tmp (v);
02262 return tmp * a;
02263 }
02264
02265 ComplexMatrix
02266 operator * (const ComplexColumnVector& a, const RowVector& b)
02267 {
02268 ComplexRowVector tmp (b);
02269 return a * tmp;
02270 }
02271
02272 ComplexMatrix
02273 operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
02274 {
02275 ComplexMatrix retval;
02276
02277 int len = v.length ();
02278
02279 if (len != 0)
02280 {
02281 int a_len = a.length ();
02282
02283 retval.resize (len, a_len);
02284 Complex *c = retval.fortran_vec ();
02285
02286 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
02287 F77_CONST_CHAR_ARG2 ("N", 1),
02288 len, a_len, 1, 1.0, v.data (), len,
02289 a.data (), 1, 0.0, c, len
02290 F77_CHAR_ARG_LEN (1)
02291 F77_CHAR_ARG_LEN (1)));
02292
02293 if (f77_exception_encountered)
02294 (*current_liboctave_error_handler)
02295 ("unrecoverable error in zgemm");
02296 }
02297
02298 return retval;
02299 }
02300
02301
02302
02303 ComplexMatrix&
02304 ComplexMatrix::operator += (const DiagMatrix& a)
02305 {
02306 int nr = rows ();
02307 int nc = cols ();
02308
02309 int a_nr = rows ();
02310 int a_nc = cols ();
02311
02312 if (nr != a_nr || nc != a_nc)
02313 {
02314 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
02315 return *this;
02316 }
02317
02318 for (int i = 0; i < a.length (); i++)
02319 elem (i, i) += a.elem (i, i);
02320
02321 return *this;
02322 }
02323
02324 ComplexMatrix&
02325 ComplexMatrix::operator -= (const DiagMatrix& a)
02326 {
02327 int nr = rows ();
02328 int nc = cols ();
02329
02330 int a_nr = rows ();
02331 int a_nc = cols ();
02332
02333 if (nr != a_nr || nc != a_nc)
02334 {
02335 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
02336 return *this;
02337 }
02338
02339 for (int i = 0; i < a.length (); i++)
02340 elem (i, i) -= a.elem (i, i);
02341
02342 return *this;
02343 }
02344
02345 ComplexMatrix&
02346 ComplexMatrix::operator += (const ComplexDiagMatrix& a)
02347 {
02348 int nr = rows ();
02349 int nc = cols ();
02350
02351 int a_nr = rows ();
02352 int a_nc = cols ();
02353
02354 if (nr != a_nr || nc != a_nc)
02355 {
02356 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
02357 return *this;
02358 }
02359
02360 for (int i = 0; i < a.length (); i++)
02361 elem (i, i) += a.elem (i, i);
02362
02363 return *this;
02364 }
02365
02366 ComplexMatrix&
02367 ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
02368 {
02369 int nr = rows ();
02370 int nc = cols ();
02371
02372 int a_nr = rows ();
02373 int a_nc = cols ();
02374
02375 if (nr != a_nr || nc != a_nc)
02376 {
02377 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
02378 return *this;
02379 }
02380
02381 for (int i = 0; i < a.length (); i++)
02382 elem (i, i) -= a.elem (i, i);
02383
02384 return *this;
02385 }
02386
02387
02388
02389 ComplexMatrix&
02390 ComplexMatrix::operator += (const Matrix& a)
02391 {
02392 int nr = rows ();
02393 int nc = cols ();
02394
02395 int a_nr = a.rows ();
02396 int a_nc = a.cols ();
02397
02398 if (nr != a_nr || nc != a_nc)
02399 {
02400 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
02401 return *this;
02402 }
02403
02404 if (nr == 0 || nc == 0)
02405 return *this;
02406
02407 Complex *d = fortran_vec ();
02408
02409 mx_inline_add2 (d, a.data (), length ());
02410 return *this;
02411 }
02412
02413 ComplexMatrix&
02414 ComplexMatrix::operator -= (const Matrix& a)
02415 {
02416 int nr = rows ();
02417 int nc = cols ();
02418
02419 int a_nr = a.rows ();
02420 int a_nc = a.cols ();
02421
02422 if (nr != a_nr || nc != a_nc)
02423 {
02424 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
02425 return *this;
02426 }
02427
02428 if (nr == 0 || nc == 0)
02429 return *this;
02430
02431 Complex *d = fortran_vec ();
02432
02433 mx_inline_subtract2 (d, a.data (), length ());
02434 return *this;
02435 }
02436
02437
02438
02439 boolMatrix
02440 ComplexMatrix::operator ! (void) const
02441 {
02442 int nr = rows ();
02443 int nc = cols ();
02444
02445 boolMatrix b (nr, nc);
02446
02447 for (int j = 0; j < nc; j++)
02448 for (int i = 0; i < nr; i++)
02449 b.elem (i, j) = elem (i, j) != 0.0;
02450
02451 return b;
02452 }
02453
02454
02455
02456 ComplexMatrix
02457 ComplexMatrix::map (c_c_Mapper f) const
02458 {
02459 ComplexMatrix b (*this);
02460 return b.apply (f);
02461 }
02462
02463 Matrix
02464 ComplexMatrix::map (d_c_Mapper f) const
02465 {
02466 int nr = rows ();
02467 int nc = cols ();
02468
02469 Matrix retval (nr, nc);
02470
02471 for (int j = 0; j < nc; j++)
02472 for (int i = 0; i < nr; i++)
02473 retval(i,j) = f (elem(i,j));
02474
02475 return retval;
02476 }
02477
02478 boolMatrix
02479 ComplexMatrix::map (b_c_Mapper f) const
02480 {
02481 int nr = rows ();
02482 int nc = cols ();
02483
02484 boolMatrix retval (nr, nc);
02485
02486 for (int j = 0; j < nc; j++)
02487 for (int i = 0; i < nr; i++)
02488 retval(i,j) = f (elem(i,j));
02489
02490 return retval;
02491 }
02492
02493 ComplexMatrix&
02494 ComplexMatrix::apply (c_c_Mapper f)
02495 {
02496 Complex *d = fortran_vec ();
02497
02498 for (int i = 0; i < length (); i++)
02499 d[i] = f (d[i]);
02500
02501 return *this;
02502 }
02503
02504 bool
02505 ComplexMatrix::any_element_is_inf_or_nan (void) const
02506 {
02507 int nr = rows ();
02508 int nc = cols ();
02509
02510 for (int j = 0; j < nc; j++)
02511 for (int i = 0; i < nr; i++)
02512 {
02513 Complex val = elem (i, j);
02514 if (xisinf (val) || xisnan (val))
02515 return true;
02516 }
02517
02518 return false;
02519 }
02520
02521
02522
02523 bool
02524 ComplexMatrix::all_elements_are_real (void) const
02525 {
02526 int nr = rows ();
02527 int nc = cols ();
02528
02529 for (int j = 0; j < nc; j++)
02530 {
02531 for (int i = 0; i < nr; i++)
02532 {
02533 double ip = imag (elem (i, j));
02534
02535 if (ip != 0.0 || lo_ieee_signbit (ip))
02536 return false;
02537 }
02538 }
02539
02540 return true;
02541 }
02542
02543
02544
02545
02546
02547 bool
02548 ComplexMatrix::all_integers (double& max_val, double& min_val) const
02549 {
02550 int nr = rows ();
02551 int nc = cols ();
02552
02553 if (nr > 0 && nc > 0)
02554 {
02555 Complex val = elem (0, 0);
02556
02557 double r_val = real (val);
02558 double i_val = imag (val);
02559
02560 max_val = r_val;
02561 min_val = r_val;
02562
02563 if (i_val > max_val)
02564 max_val = i_val;
02565
02566 if (i_val < max_val)
02567 min_val = i_val;
02568 }
02569 else
02570 return false;
02571
02572 for (int j = 0; j < nc; j++)
02573 for (int i = 0; i < nr; i++)
02574 {
02575 Complex val = elem (i, j);
02576
02577 double r_val = real (val);
02578 double i_val = imag (val);
02579
02580 if (r_val > max_val)
02581 max_val = r_val;
02582
02583 if (i_val > max_val)
02584 max_val = i_val;
02585
02586 if (r_val < min_val)
02587 min_val = r_val;
02588
02589 if (i_val < min_val)
02590 min_val = i_val;
02591
02592 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
02593 return false;
02594 }
02595
02596 return true;
02597 }
02598
02599 bool
02600 ComplexMatrix::too_large_for_float (void) const
02601 {
02602 int nr = rows ();
02603 int nc = cols ();
02604
02605 for (int j = 0; j < nc; j++)
02606 for (int i = 0; i < nr; i++)
02607 {
02608 Complex val = elem (i, j);
02609
02610 double r_val = real (val);
02611 double i_val = imag (val);
02612
02613 if (r_val > FLT_MAX
02614 || i_val > FLT_MAX
02615 || r_val < FLT_MIN
02616 || i_val < FLT_MIN)
02617 return true;
02618 }
02619
02620 return false;
02621 }
02622
02623
02624
02625
02626 boolMatrix
02627 ComplexMatrix::all (int dim) const
02628 {
02629 MX_ALL_OP (dim);
02630 }
02631
02632 boolMatrix
02633 ComplexMatrix::any (int dim) const
02634 {
02635 MX_ANY_OP (dim);
02636 }
02637
02638 ComplexMatrix
02639 ComplexMatrix::cumprod (int dim) const
02640 {
02641 MX_CUMULATIVE_OP (ComplexMatrix, Complex, *=);
02642 }
02643
02644 ComplexMatrix
02645 ComplexMatrix::cumsum (int dim) const
02646 {
02647 MX_CUMULATIVE_OP (ComplexMatrix, Complex, +=);
02648 }
02649
02650 ComplexMatrix
02651 ComplexMatrix::prod (int dim) const
02652 {
02653 MX_REDUCTION_OP (ComplexMatrix, *=, 1.0, 1.0);
02654 }
02655
02656 ComplexMatrix
02657 ComplexMatrix::sum (int dim) const
02658 {
02659 MX_REDUCTION_OP (ComplexMatrix, +=, 0.0, 0.0);
02660 }
02661
02662 ComplexMatrix
02663 ComplexMatrix::sumsq (int dim) const
02664 {
02665 #define ROW_EXPR \
02666 Complex d = elem (i, j); \
02667 retval.elem (i, 0) += d * conj (d)
02668
02669 #define COL_EXPR \
02670 Complex d = elem (i, j); \
02671 retval.elem (0, j) += d * conj (d)
02672
02673 MX_BASE_REDUCTION_OP (ComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
02674
02675 #undef ROW_EXPR
02676 #undef COL_EXPR
02677 }
02678
02679 Matrix ComplexMatrix::abs (void) const
02680 {
02681 int nr = rows ();
02682 int nc = cols ();
02683
02684 Matrix retval (nr, nc);
02685
02686 for (int j = 0; j < nc; j++)
02687 for (int i = 0; i < nr; i++)
02688 retval (i, j) = ::abs (elem (i, j));
02689
02690 return retval;
02691 }
02692
02693 ComplexColumnVector
02694 ComplexMatrix::diag (void) const
02695 {
02696 return diag (0);
02697 }
02698
02699 ComplexColumnVector
02700 ComplexMatrix::diag (int k) const
02701 {
02702 int nnr = rows ();
02703 int nnc = cols ();
02704 if (k > 0)
02705 nnc -= k;
02706 else if (k < 0)
02707 nnr += k;
02708
02709 ComplexColumnVector d;
02710
02711 if (nnr > 0 && nnc > 0)
02712 {
02713 int ndiag = (nnr < nnc) ? nnr : nnc;
02714
02715 d.resize (ndiag);
02716
02717 if (k > 0)
02718 {
02719 for (int i = 0; i < ndiag; i++)
02720 d.elem (i) = elem (i, i+k);
02721 }
02722 else if (k < 0)
02723 {
02724 for (int i = 0; i < ndiag; i++)
02725 d.elem (i) = elem (i-k, i);
02726 }
02727 else
02728 {
02729 for (int i = 0; i < ndiag; i++)
02730 d.elem (i) = elem (i, i);
02731 }
02732 }
02733 else
02734 (*current_liboctave_error_handler)
02735 ("diag: requested diagonal out of range");
02736
02737 return d;
02738 }
02739
02740 bool
02741 ComplexMatrix::row_is_real_only (int i) const
02742 {
02743 bool retval = true;
02744
02745 int nc = columns ();
02746
02747 for (int j = 0; j < nc; j++)
02748 {
02749 if (imag (elem (i, j)) != 0.0)
02750 {
02751 retval = false;
02752 break;
02753 }
02754 }
02755
02756 return retval;
02757 }
02758
02759 bool
02760 ComplexMatrix::column_is_real_only (int j) const
02761 {
02762 bool retval = true;
02763
02764 int nr = rows ();
02765
02766 for (int i = 0; i < nr; i++)
02767 {
02768 if (imag (elem (i, j)) != 0.0)
02769 {
02770 retval = false;
02771 break;
02772 }
02773 }
02774
02775 return retval;
02776 }
02777
02778 ComplexColumnVector
02779 ComplexMatrix::row_min (void) const
02780 {
02781 Array<int> dummy_idx;
02782 return row_min (dummy_idx);
02783 }
02784
02785 ComplexColumnVector
02786 ComplexMatrix::row_min (Array<int>& idx_arg) const
02787 {
02788 ComplexColumnVector result;
02789
02790 int nr = rows ();
02791 int nc = cols ();
02792
02793 if (nr > 0 && nc > 0)
02794 {
02795 result.resize (nr);
02796 idx_arg.resize (nr);
02797
02798 for (int i = 0; i < nr; i++)
02799 {
02800 bool real_only = row_is_real_only (i);
02801
02802 int idx_j;
02803
02804 Complex tmp_min;
02805
02806 double abs_min = octave_NaN;
02807
02808 for (idx_j = 0; idx_j < nc; idx_j++)
02809 {
02810 tmp_min = elem (i, idx_j);
02811
02812 if (! octave_is_NaN_or_NA (tmp_min))
02813 {
02814 abs_min = real_only ? real (tmp_min) : ::abs (tmp_min);
02815 break;
02816 }
02817 }
02818
02819 for (int j = idx_j+1; j < nc; j++)
02820 {
02821 Complex tmp = elem (i, j);
02822
02823 if (octave_is_NaN_or_NA (tmp))
02824 continue;
02825
02826 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
02827
02828 if (abs_tmp < abs_min)
02829 {
02830 idx_j = j;
02831 tmp_min = tmp;
02832 abs_min = abs_tmp;
02833 }
02834 }
02835
02836 if (octave_is_NaN_or_NA (tmp_min))
02837 {
02838 result.elem (i) = Complex_NaN_result;
02839 idx_arg.elem (i) = 0;
02840 }
02841 else
02842 {
02843 result.elem (i) = tmp_min;
02844 idx_arg.elem (i) = idx_j;
02845 }
02846 }
02847 }
02848
02849 return result;
02850 }
02851
02852 ComplexColumnVector
02853 ComplexMatrix::row_max (void) const
02854 {
02855 Array<int> dummy_idx;
02856 return row_max (dummy_idx);
02857 }
02858
02859 ComplexColumnVector
02860 ComplexMatrix::row_max (Array<int>& idx_arg) const
02861 {
02862 ComplexColumnVector result;
02863
02864 int nr = rows ();
02865 int nc = cols ();
02866
02867 if (nr > 0 && nc > 0)
02868 {
02869 result.resize (nr);
02870 idx_arg.resize (nr);
02871
02872 for (int i = 0; i < nr; i++)
02873 {
02874 bool real_only = row_is_real_only (i);
02875
02876 int idx_j;
02877
02878 Complex tmp_max;
02879
02880 double abs_max = octave_NaN;
02881
02882 for (idx_j = 0; idx_j < nc; idx_j++)
02883 {
02884 tmp_max = elem (i, idx_j);
02885
02886 if (! octave_is_NaN_or_NA (tmp_max))
02887 {
02888 abs_max = real_only ? real (tmp_max) : ::abs (tmp_max);
02889 break;
02890 }
02891 }
02892
02893 for (int j = idx_j+1; j < nc; j++)
02894 {
02895 Complex tmp = elem (i, j);
02896
02897 if (octave_is_NaN_or_NA (tmp))
02898 continue;
02899
02900 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
02901
02902 if (abs_tmp > abs_max)
02903 {
02904 idx_j = j;
02905 tmp_max = tmp;
02906 abs_max = abs_tmp;
02907 }
02908 }
02909
02910 if (octave_is_NaN_or_NA (tmp_max))
02911 {
02912 result.elem (i) = Complex_NaN_result;
02913 idx_arg.elem (i) = 0;
02914 }
02915 else
02916 {
02917 result.elem (i) = tmp_max;
02918 idx_arg.elem (i) = idx_j;
02919 }
02920 }
02921 }
02922
02923 return result;
02924 }
02925
02926 ComplexRowVector
02927 ComplexMatrix::column_min (void) const
02928 {
02929 Array<int> dummy_idx;
02930 return column_min (dummy_idx);
02931 }
02932
02933 ComplexRowVector
02934 ComplexMatrix::column_min (Array<int>& idx_arg) const
02935 {
02936 ComplexRowVector result;
02937
02938 int nr = rows ();
02939 int nc = cols ();
02940
02941 if (nr > 0 && nc > 0)
02942 {
02943 result.resize (nc);
02944 idx_arg.resize (nc);
02945
02946 for (int j = 0; j < nc; j++)
02947 {
02948 bool real_only = column_is_real_only (j);
02949
02950 int idx_i;
02951
02952 Complex tmp_min;
02953
02954 double abs_min = octave_NaN;
02955
02956 for (idx_i = 0; idx_i < nr; idx_i++)
02957 {
02958 tmp_min = elem (idx_i, j);
02959
02960 if (! octave_is_NaN_or_NA (tmp_min))
02961 {
02962 abs_min = real_only ? real (tmp_min) : ::abs (tmp_min);
02963 break;
02964 }
02965 }
02966
02967 for (int i = idx_i+1; i < nr; i++)
02968 {
02969 Complex tmp = elem (i, j);
02970
02971 if (octave_is_NaN_or_NA (tmp))
02972 continue;
02973
02974 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
02975
02976 if (abs_tmp < abs_min)
02977 {
02978 idx_i = i;
02979 tmp_min = tmp;
02980 abs_min = abs_tmp;
02981 }
02982 }
02983
02984 if (octave_is_NaN_or_NA (tmp_min))
02985 {
02986 result.elem (j) = Complex_NaN_result;
02987 idx_arg.elem (j) = 0;
02988 }
02989 else
02990 {
02991 result.elem (j) = tmp_min;
02992 idx_arg.elem (j) = idx_i;
02993 }
02994 }
02995 }
02996
02997 return result;
02998 }
02999
03000 ComplexRowVector
03001 ComplexMatrix::column_max (void) const
03002 {
03003 Array<int> dummy_idx;
03004 return column_max (dummy_idx);
03005 }
03006
03007 ComplexRowVector
03008 ComplexMatrix::column_max (Array<int>& idx_arg) const
03009 {
03010 ComplexRowVector result;
03011
03012 int nr = rows ();
03013 int nc = cols ();
03014
03015 if (nr > 0 && nc > 0)
03016 {
03017 result.resize (nc);
03018 idx_arg.resize (nc);
03019
03020 for (int j = 0; j < nc; j++)
03021 {
03022 bool real_only = column_is_real_only (j);
03023
03024 int idx_i;
03025
03026 Complex tmp_max;
03027
03028 double abs_max = octave_NaN;
03029
03030 for (idx_i = 0; idx_i < nr; idx_i++)
03031 {
03032 tmp_max = elem (idx_i, j);
03033
03034 if (! octave_is_NaN_or_NA (tmp_max))
03035 {
03036 abs_max = real_only ? real (tmp_max) : ::abs (tmp_max);
03037 break;
03038 }
03039 }
03040
03041 for (int i = idx_i+1; i < nr; i++)
03042 {
03043 Complex tmp = elem (i, j);
03044
03045 if (octave_is_NaN_or_NA (tmp))
03046 continue;
03047
03048 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
03049
03050 if (abs_tmp > abs_max)
03051 {
03052 idx_i = i;
03053 tmp_max = tmp;
03054 abs_max = abs_tmp;
03055 }
03056 }
03057
03058 if (octave_is_NaN_or_NA (tmp_max))
03059 {
03060 result.elem (j) = Complex_NaN_result;
03061 idx_arg.elem (j) = 0;
03062 }
03063 else
03064 {
03065 result.elem (j) = tmp_max;
03066 idx_arg.elem (j) = idx_i;
03067 }
03068 }
03069 }
03070
03071 return result;
03072 }
03073
03074
03075
03076 std::ostream&
03077 operator << (std::ostream& os, const ComplexMatrix& a)
03078 {
03079 for (int i = 0; i < a.rows (); i++)
03080 {
03081 for (int j = 0; j < a.cols (); j++)
03082 {
03083 os << " ";
03084 octave_write_complex (os, a.elem (i, j));
03085 }
03086 os << "\n";
03087 }
03088 return os;
03089 }
03090
03091 std::istream&
03092 operator >> (std::istream& is, ComplexMatrix& a)
03093 {
03094 int nr = a.rows ();
03095 int nc = a.cols ();
03096
03097 if (nr < 1 || nc < 1)
03098 is.clear (std::ios::badbit);
03099 else
03100 {
03101 Complex tmp;
03102 for (int i = 0; i < nr; i++)
03103 for (int j = 0; j < nc; j++)
03104 {
03105 tmp = octave_read_complex (is);
03106 if (is)
03107 a.elem (i, j) = tmp;
03108 else
03109 goto done;
03110 }
03111 }
03112
03113 done:
03114
03115 return is;
03116 }
03117
03118 ComplexMatrix
03119 Givens (const Complex& x, const Complex& y)
03120 {
03121 double cc;
03122 Complex cs, temp_r;
03123
03124 F77_FUNC (zlartg, ZLARTG) (x, y, cc, cs, temp_r);
03125
03126 ComplexMatrix g (2, 2);
03127
03128 g.elem (0, 0) = cc;
03129 g.elem (1, 1) = cc;
03130 g.elem (0, 1) = cs;
03131 g.elem (1, 0) = -conj (cs);
03132
03133 return g;
03134 }
03135
03136 ComplexMatrix
03137 Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
03138 const ComplexMatrix& c)
03139 {
03140 ComplexMatrix retval;
03141
03142
03143
03144
03145
03146
03147 ComplexSCHUR as (a, "U");
03148 ComplexSCHUR bs (b, "U");
03149
03150
03151
03152 ComplexMatrix ua = as.unitary_matrix ();
03153 ComplexMatrix sch_a = as.schur_matrix ();
03154
03155 ComplexMatrix ub = bs.unitary_matrix ();
03156 ComplexMatrix sch_b = bs.schur_matrix ();
03157
03158 ComplexMatrix cx = ua.hermitian () * c * ub;
03159
03160
03161
03162
03163 int a_nr = a.rows ();
03164 int b_nr = b.rows ();
03165
03166 double scale;
03167 int info;
03168
03169 Complex *pa = sch_a.fortran_vec ();
03170 Complex *pb = sch_b.fortran_vec ();
03171 Complex *px = cx.fortran_vec ();
03172
03173 F77_XFCN (ztrsyl, ZTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
03174 F77_CONST_CHAR_ARG2 ("N", 1),
03175 1, a_nr, b_nr, pa, a_nr, pb,
03176 b_nr, px, a_nr, scale, info
03177 F77_CHAR_ARG_LEN (1)
03178 F77_CHAR_ARG_LEN (1)));
03179
03180 if (f77_exception_encountered)
03181 (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl");
03182 else
03183 {
03184
03185
03186 retval = -ua * cx * ub.hermitian ();
03187 }
03188
03189 return retval;
03190 }
03191
03192 ComplexMatrix
03193 operator * (const ComplexMatrix& m, const Matrix& a)
03194 {
03195 ComplexMatrix tmp (a);
03196 return m * tmp;
03197 }
03198
03199 ComplexMatrix
03200 operator * (const Matrix& m, const ComplexMatrix& a)
03201 {
03202 ComplexMatrix tmp (m);
03203 return tmp * a;
03204 }
03205
03206 ComplexMatrix
03207 operator * (const ComplexMatrix& m, const ComplexMatrix& a)
03208 {
03209 ComplexMatrix retval;
03210
03211 int nr = m.rows ();
03212 int nc = m.cols ();
03213
03214 int a_nr = a.rows ();
03215 int a_nc = a.cols ();
03216
03217 if (nc != a_nr)
03218 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
03219 else
03220 {
03221 if (nr == 0 || nc == 0 || a_nc == 0)
03222 retval.resize (nr, a_nc, 0.0);
03223 else
03224 {
03225 int ld = nr;
03226 int lda = a.rows ();
03227
03228 retval.resize (nr, a_nc);
03229 Complex *c = retval.fortran_vec ();
03230
03231 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
03232 F77_CONST_CHAR_ARG2 ("N", 1),
03233 nr, a_nc, nc, 1.0, m.data (),
03234 ld, a.data (), lda, 0.0, c, nr
03235 F77_CHAR_ARG_LEN (1)
03236 F77_CHAR_ARG_LEN (1)));
03237
03238 if (f77_exception_encountered)
03239 (*current_liboctave_error_handler)
03240 ("unrecoverable error in zgemm");
03241 }
03242 }
03243
03244 return retval;
03245 }
03246
03247
03248
03249
03250 #define EMPTY_RETURN_CHECK(T) \
03251 if (nr == 0 || nc == 0) \
03252 return T (nr, nc);
03253
03254 ComplexMatrix
03255 min (const Complex& c, const ComplexMatrix& m)
03256 {
03257 int nr = m.rows ();
03258 int nc = m.columns ();
03259
03260 EMPTY_RETURN_CHECK (ComplexMatrix);
03261
03262 ComplexMatrix result (nr, nc);
03263
03264 for (int j = 0; j < nc; j++)
03265 for (int i = 0; i < nr; i++)
03266 {
03267 OCTAVE_QUIT;
03268 result (i, j) = xmin (c, m (i, j));
03269 }
03270
03271 return result;
03272 }
03273
03274 ComplexMatrix
03275 min (const ComplexMatrix& m, const Complex& c)
03276 {
03277 int nr = m.rows ();
03278 int nc = m.columns ();
03279
03280 EMPTY_RETURN_CHECK (ComplexMatrix);
03281
03282 ComplexMatrix result (nr, nc);
03283
03284 for (int j = 0; j < nc; j++)
03285 for (int i = 0; i < nr; i++)
03286 {
03287 OCTAVE_QUIT;
03288 result (i, j) = xmin (m (i, j), c);
03289 }
03290
03291 return result;
03292 }
03293
03294 ComplexMatrix
03295 min (const ComplexMatrix& a, const ComplexMatrix& b)
03296 {
03297 int nr = a.rows ();
03298 int nc = a.columns ();
03299
03300 if (nr != b.rows () || nc != b.columns ())
03301 {
03302 (*current_liboctave_error_handler)
03303 ("two-arg min expecting args of same size");
03304 return ComplexMatrix ();
03305 }
03306
03307 EMPTY_RETURN_CHECK (ComplexMatrix);
03308
03309 ComplexMatrix result (nr, nc);
03310
03311 for (int j = 0; j < nc; j++)
03312 {
03313 int columns_are_real_only = 1;
03314 for (int i = 0; i < nr; i++)
03315 {
03316 OCTAVE_QUIT;
03317 if (imag (a (i, j)) != 0.0 || imag (b (i, j)) != 0.0)
03318 {
03319 columns_are_real_only = 0;
03320 break;
03321 }
03322 }
03323
03324 if (columns_are_real_only)
03325 {
03326 for (int i = 0; i < nr; i++)
03327 result (i, j) = xmin (real (a (i, j)), real (b (i, j)));
03328 }
03329 else
03330 {
03331 for (int i = 0; i < nr; i++)
03332 {
03333 OCTAVE_QUIT;
03334 result (i, j) = xmin (a (i, j), b (i, j));
03335 }
03336 }
03337 }
03338
03339 return result;
03340 }
03341
03342 ComplexMatrix
03343 max (const Complex& c, const ComplexMatrix& m)
03344 {
03345 int nr = m.rows ();
03346 int nc = m.columns ();
03347
03348 EMPTY_RETURN_CHECK (ComplexMatrix);
03349
03350 ComplexMatrix result (nr, nc);
03351
03352 for (int j = 0; j < nc; j++)
03353 for (int i = 0; i < nr; i++)
03354 {
03355 OCTAVE_QUIT;
03356 result (i, j) = xmax (c, m (i, j));
03357 }
03358
03359 return result;
03360 }
03361
03362 ComplexMatrix
03363 max (const ComplexMatrix& m, const Complex& c)
03364 {
03365 int nr = m.rows ();
03366 int nc = m.columns ();
03367
03368 EMPTY_RETURN_CHECK (ComplexMatrix);
03369
03370 ComplexMatrix result (nr, nc);
03371
03372 for (int j = 0; j < nc; j++)
03373 for (int i = 0; i < nr; i++)
03374 {
03375 OCTAVE_QUIT;
03376 result (i, j) = xmax (m (i, j), c);
03377 }
03378
03379 return result;
03380 }
03381
03382 ComplexMatrix
03383 max (const ComplexMatrix& a, const ComplexMatrix& b)
03384 {
03385 int nr = a.rows ();
03386 int nc = a.columns ();
03387
03388 if (nr != b.rows () || nc != b.columns ())
03389 {
03390 (*current_liboctave_error_handler)
03391 ("two-arg max expecting args of same size");
03392 return ComplexMatrix ();
03393 }
03394
03395 EMPTY_RETURN_CHECK (ComplexMatrix);
03396
03397 ComplexMatrix result (nr, nc);
03398
03399 for (int j = 0; j < nc; j++)
03400 {
03401 int columns_are_real_only = 1;
03402 for (int i = 0; i < nr; i++)
03403 {
03404 OCTAVE_QUIT;
03405 if (imag (a (i, j)) != 0.0 || imag (b (i, j)) != 0.0)
03406 {
03407 columns_are_real_only = 0;
03408 break;
03409 }
03410 }
03411
03412 if (columns_are_real_only)
03413 {
03414 for (int i = 0; i < nr; i++)
03415 {
03416 OCTAVE_QUIT;
03417 result (i, j) = xmax (real (a (i, j)), real (b (i, j)));
03418 }
03419 }
03420 else
03421 {
03422 for (int i = 0; i < nr; i++)
03423 {
03424 OCTAVE_QUIT;
03425 result (i, j) = xmax (a (i, j), b (i, j));
03426 }
03427 }
03428 }
03429
03430 return result;
03431 }
03432
03433 MS_CMP_OPS(ComplexMatrix, real, Complex, real)
03434 MS_BOOL_OPS(ComplexMatrix, Complex, 0.0)
03435
03436 SM_CMP_OPS(Complex, real, ComplexMatrix, real)
03437 SM_BOOL_OPS(Complex, ComplexMatrix, 0.0)
03438
03439 MM_CMP_OPS(ComplexMatrix, real, ComplexMatrix, real)
03440 MM_BOOL_OPS(ComplexMatrix, ComplexMatrix, 0.0)
03441
03442
03443
03444
03445
03446