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 #include "Array-util.h"
00037 #include "byte-swap.h"
00038 #include "dMatrix.h"
00039 #include "dbleAEPBAL.h"
00040 #include "dbleDET.h"
00041 #include "dbleSCHUR.h"
00042 #include "dbleSVD.h"
00043 #include "f77-fcn.h"
00044 #include "lo-error.h"
00045 #include "lo-ieee.h"
00046 #include "lo-mappers.h"
00047 #include "lo-utils.h"
00048 #include "mx-base.h"
00049 #include "mx-m-dm.h"
00050 #include "mx-dm-m.h"
00051 #include "mx-inlines.cc"
00052 #include "oct-cmplx.h"
00053 #include "quit.h"
00054
00055 #if defined (HAVE_FFTW3)
00056 #include "oct-fftw.h"
00057 #endif
00058
00059
00060
00061 extern "C"
00062 {
00063 F77_RET_T
00064 F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
00065 const int&, double*, const int&, int&,
00066 int&, double*, int&
00067 F77_CHAR_ARG_LEN_DECL);
00068
00069 F77_RET_T
00070 F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
00071 F77_CONST_CHAR_ARG_DECL,
00072 const int&, const int&, const int&, double*,
00073 const int&, double*, const int&, int&
00074 F77_CHAR_ARG_LEN_DECL
00075 F77_CHAR_ARG_LEN_DECL);
00076
00077
00078 F77_RET_T
00079 F77_FUNC (dgemm, DGEMM) (F77_CONST_CHAR_ARG_DECL,
00080 F77_CONST_CHAR_ARG_DECL,
00081 const int&, const int&, const int&,
00082 const double&, const double*, const int&,
00083 const double*, const int&, const double&,
00084 double*, const int&
00085 F77_CHAR_ARG_LEN_DECL
00086 F77_CHAR_ARG_LEN_DECL);
00087
00088 F77_RET_T
00089 F77_FUNC (dgetrf, DGETRF) (const int&, const int&, double*, const int&,
00090 int*, int&);
00091
00092 F77_RET_T
00093 F77_FUNC (dgetrs, DGETRS) (F77_CONST_CHAR_ARG_DECL, const int&, const int&,
00094 const double*, const int&,
00095 const int*, double*, const int&, int&
00096 F77_CHAR_ARG_LEN_DECL);
00097
00098 F77_RET_T
00099 F77_FUNC (dgetri, DGETRI) (const int&, double*, const int&, const int*,
00100 double*, const int&, int&);
00101
00102 F77_RET_T
00103 F77_FUNC (dgecon, DGECON) (F77_CONST_CHAR_ARG_DECL, const int&, double*,
00104 const int&, const double&, double&,
00105 double*, int*, int&
00106 F77_CHAR_ARG_LEN_DECL);
00107
00108 F77_RET_T
00109 F77_FUNC (dgelss, DGELSS) (const int&, const int&, const int&,
00110 double*, const int&, double*,
00111 const int&, double*, double&, int&,
00112 double*, const int&, int&);
00113
00114
00115
00116
00117
00118
00119 F77_RET_T
00120 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
00121
00122 F77_RET_T
00123 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
00124
00125 F77_RET_T
00126 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
00127
00128 F77_RET_T
00129 F77_FUNC (dlartg, DLARTG) (const double&, const double&, double&,
00130 double&, double&);
00131
00132 F77_RET_T
00133 F77_FUNC (dtrsyl, DTRSYL) (F77_CONST_CHAR_ARG_DECL,
00134 F77_CONST_CHAR_ARG_DECL,
00135 const int&, const int&, const int&,
00136 const double*, const int&, const double*,
00137 const int&, const double*, const int&,
00138 double&, int&
00139 F77_CHAR_ARG_LEN_DECL
00140 F77_CHAR_ARG_LEN_DECL);
00141
00142 F77_RET_T
00143 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL, const int&,
00144 const int&, const double*,
00145 const int&, double*, double&
00146 F77_CHAR_ARG_LEN_DECL);
00147 }
00148
00149
00150
00151 Matrix::Matrix (const RowVector& rv)
00152 : MArray2<double> (1, rv.length (), 0.0)
00153 {
00154 for (int i = 0; i < rv.length (); i++)
00155 elem (0, i) = rv.elem (i);
00156 }
00157
00158 Matrix::Matrix (const ColumnVector& cv)
00159 : MArray2<double> (cv.length (), 1, 0.0)
00160 {
00161 for (int i = 0; i < cv.length (); i++)
00162 elem (i, 0) = cv.elem (i);
00163 }
00164
00165 Matrix::Matrix (const DiagMatrix& a)
00166 : MArray2<double> (a.rows (), a.cols (), 0.0)
00167 {
00168 for (int i = 0; i < a.length (); i++)
00169 elem (i, i) = a.elem (i, i);
00170 }
00171
00172
00173
00174
00175 Matrix::Matrix (const boolMatrix& a)
00176 : MArray2<double> (a.rows (), a.cols ())
00177 {
00178 for (int i = 0; i < a.rows (); i++)
00179 for (int j = 0; j < a.cols (); j++)
00180 elem (i, j) = a.elem (i, j);
00181 }
00182
00183 Matrix::Matrix (const charMatrix& a)
00184 : MArray2<double> (a.rows (), a.cols ())
00185 {
00186 for (int i = 0; i < a.rows (); i++)
00187 for (int j = 0; j < a.cols (); j++)
00188 elem (i, j) = a.elem (i, j);
00189 }
00190
00191 bool
00192 Matrix::operator == (const Matrix& a) const
00193 {
00194 if (rows () != a.rows () || cols () != a.cols ())
00195 return false;
00196
00197 return mx_inline_equal (data (), a.data (), length ());
00198 }
00199
00200 bool
00201 Matrix::operator != (const Matrix& a) const
00202 {
00203 return !(*this == a);
00204 }
00205
00206 bool
00207 Matrix::is_symmetric (void) const
00208 {
00209 if (is_square () && rows () > 0)
00210 {
00211 for (int i = 0; i < rows (); i++)
00212 for (int j = i+1; j < cols (); j++)
00213 if (elem (i, j) != elem (j, i))
00214 return false;
00215
00216 return true;
00217 }
00218
00219 return false;
00220 }
00221
00222 Matrix&
00223 Matrix::insert (const Matrix& a, int r, int c)
00224 {
00225 Array2<double>::insert (a, r, c);
00226 return *this;
00227 }
00228
00229 Matrix&
00230 Matrix::insert (const RowVector& a, int r, int c)
00231 {
00232 int a_len = a.length ();
00233
00234 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
00235 {
00236 (*current_liboctave_error_handler) ("range error for insert");
00237 return *this;
00238 }
00239
00240 if (a_len > 0)
00241 {
00242 make_unique ();
00243
00244 for (int i = 0; i < a_len; i++)
00245 xelem (r, c+i) = a.elem (i);
00246 }
00247
00248 return *this;
00249 }
00250
00251 Matrix&
00252 Matrix::insert (const ColumnVector& a, int r, int c)
00253 {
00254 int a_len = a.length ();
00255
00256 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
00257 {
00258 (*current_liboctave_error_handler) ("range error for insert");
00259 return *this;
00260 }
00261
00262 if (a_len > 0)
00263 {
00264 make_unique ();
00265
00266 for (int i = 0; i < a_len; i++)
00267 xelem (r+i, c) = a.elem (i);
00268 }
00269
00270 return *this;
00271 }
00272
00273 Matrix&
00274 Matrix::insert (const DiagMatrix& a, int r, int c)
00275 {
00276 int a_nr = a.rows ();
00277 int a_nc = a.cols ();
00278
00279 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
00280 {
00281 (*current_liboctave_error_handler) ("range error for insert");
00282 return *this;
00283 }
00284
00285 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
00286
00287 int a_len = a.length ();
00288
00289 if (a_len > 0)
00290 {
00291 make_unique ();
00292
00293 for (int i = 0; i < a_len; i++)
00294 xelem (r+i, c+i) = a.elem (i, i);
00295 }
00296
00297 return *this;
00298 }
00299
00300 Matrix&
00301 Matrix::fill (double val)
00302 {
00303 int nr = rows ();
00304 int nc = cols ();
00305
00306 if (nr > 0 && nc > 0)
00307 {
00308 make_unique ();
00309
00310 for (int j = 0; j < nc; j++)
00311 for (int i = 0; i < nr; i++)
00312 xelem (i, j) = val;
00313 }
00314
00315 return *this;
00316 }
00317
00318 Matrix&
00319 Matrix::fill (double val, int r1, int c1, int r2, int c2)
00320 {
00321 int nr = rows ();
00322 int nc = cols ();
00323
00324 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
00325 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
00326 {
00327 (*current_liboctave_error_handler) ("range error for fill");
00328 return *this;
00329 }
00330
00331 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00332 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00333
00334 if (r2 >= r1 && c2 >= c1)
00335 {
00336 make_unique ();
00337
00338 for (int j = c1; j <= c2; j++)
00339 for (int i = r1; i <= r2; i++)
00340 xelem (i, j) = val;
00341 }
00342
00343 return *this;
00344 }
00345
00346 Matrix
00347 Matrix::append (const Matrix& a) const
00348 {
00349 int nr = rows ();
00350 int nc = cols ();
00351 if (nr != a.rows ())
00352 {
00353 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00354 return Matrix ();
00355 }
00356
00357 int nc_insert = nc;
00358 Matrix retval (nr, nc + a.cols ());
00359 retval.insert (*this, 0, 0);
00360 retval.insert (a, 0, nc_insert);
00361 return retval;
00362 }
00363
00364 Matrix
00365 Matrix::append (const RowVector& a) const
00366 {
00367 int nr = rows ();
00368 int nc = cols ();
00369 if (nr != 1)
00370 {
00371 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00372 return Matrix ();
00373 }
00374
00375 int nc_insert = nc;
00376 Matrix retval (nr, nc + a.length ());
00377 retval.insert (*this, 0, 0);
00378 retval.insert (a, 0, nc_insert);
00379 return retval;
00380 }
00381
00382 Matrix
00383 Matrix::append (const ColumnVector& a) const
00384 {
00385 int nr = rows ();
00386 int nc = cols ();
00387 if (nr != a.length ())
00388 {
00389 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00390 return Matrix ();
00391 }
00392
00393 int nc_insert = nc;
00394 Matrix retval (nr, nc + 1);
00395 retval.insert (*this, 0, 0);
00396 retval.insert (a, 0, nc_insert);
00397 return retval;
00398 }
00399
00400 Matrix
00401 Matrix::append (const DiagMatrix& a) const
00402 {
00403 int nr = rows ();
00404 int nc = cols ();
00405 if (nr != a.rows ())
00406 {
00407 (*current_liboctave_error_handler) ("row dimension mismatch for append");
00408 return *this;
00409 }
00410
00411 int nc_insert = nc;
00412 Matrix retval (nr, nc + a.cols ());
00413 retval.insert (*this, 0, 0);
00414 retval.insert (a, 0, nc_insert);
00415 return retval;
00416 }
00417
00418 Matrix
00419 Matrix::stack (const Matrix& a) const
00420 {
00421 int nr = rows ();
00422 int nc = cols ();
00423 if (nc != a.cols ())
00424 {
00425 (*current_liboctave_error_handler)
00426 ("column dimension mismatch for stack");
00427 return Matrix ();
00428 }
00429
00430 int nr_insert = nr;
00431 Matrix retval (nr + a.rows (), nc);
00432 retval.insert (*this, 0, 0);
00433 retval.insert (a, nr_insert, 0);
00434 return retval;
00435 }
00436
00437 Matrix
00438 Matrix::stack (const RowVector& a) const
00439 {
00440 int nr = rows ();
00441 int nc = cols ();
00442 if (nc != a.length ())
00443 {
00444 (*current_liboctave_error_handler)
00445 ("column dimension mismatch for stack");
00446 return Matrix ();
00447 }
00448
00449 int nr_insert = nr;
00450 Matrix retval (nr + 1, nc);
00451 retval.insert (*this, 0, 0);
00452 retval.insert (a, nr_insert, 0);
00453 return retval;
00454 }
00455
00456 Matrix
00457 Matrix::stack (const ColumnVector& a) const
00458 {
00459 int nr = rows ();
00460 int nc = cols ();
00461 if (nc != 1)
00462 {
00463 (*current_liboctave_error_handler)
00464 ("column dimension mismatch for stack");
00465 return Matrix ();
00466 }
00467
00468 int nr_insert = nr;
00469 Matrix retval (nr + a.length (), nc);
00470 retval.insert (*this, 0, 0);
00471 retval.insert (a, nr_insert, 0);
00472 return retval;
00473 }
00474
00475 Matrix
00476 Matrix::stack (const DiagMatrix& a) const
00477 {
00478 int nr = rows ();
00479 int nc = cols ();
00480 if (nc != a.cols ())
00481 {
00482 (*current_liboctave_error_handler)
00483 ("column dimension mismatch for stack");
00484 return Matrix ();
00485 }
00486
00487 int nr_insert = nr;
00488 Matrix retval (nr + a.rows (), nc);
00489 retval.insert (*this, 0, 0);
00490 retval.insert (a, nr_insert, 0);
00491 return retval;
00492 }
00493
00494 Matrix
00495 real (const ComplexMatrix& a)
00496 {
00497 int a_len = a.length ();
00498 Matrix retval;
00499 if (a_len > 0)
00500 retval = Matrix (mx_inline_real_dup (a.data (), a_len),
00501 a.rows (), a.cols ());
00502 return retval;
00503 }
00504
00505 Matrix
00506 imag (const ComplexMatrix& a)
00507 {
00508 int a_len = a.length ();
00509 Matrix retval;
00510 if (a_len > 0)
00511 retval = Matrix (mx_inline_imag_dup (a.data (), a_len),
00512 a.rows (), a.cols ());
00513 return retval;
00514 }
00515
00516 Matrix
00517 Matrix::extract (int r1, int c1, int r2, int c2) const
00518 {
00519 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00520 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00521
00522 int new_r = r2 - r1 + 1;
00523 int new_c = c2 - c1 + 1;
00524
00525 Matrix result (new_r, new_c);
00526
00527 for (int j = 0; j < new_c; j++)
00528 for (int i = 0; i < new_r; i++)
00529 result.xelem (i, j) = elem (r1+i, c1+j);
00530
00531 return result;
00532 }
00533
00534 Matrix
00535 Matrix::extract_n (int r1, int c1, int nr, int nc) const
00536 {
00537 Matrix result (nr, nc);
00538
00539 for (int j = 0; j < nc; j++)
00540 for (int i = 0; i < nr; i++)
00541 result.xelem (i, j) = elem (r1+i, c1+j);
00542
00543 return result;
00544 }
00545
00546
00547
00548 RowVector
00549 Matrix::row (int i) const
00550 {
00551 int nc = cols ();
00552 if (i < 0 || i >= rows ())
00553 {
00554 (*current_liboctave_error_handler) ("invalid row selection");
00555 return RowVector ();
00556 }
00557
00558 RowVector retval (nc);
00559 for (int j = 0; j < nc; j++)
00560 retval.xelem (j) = elem (i, j);
00561
00562 return retval;
00563 }
00564
00565 RowVector
00566 Matrix::row (char *s) const
00567 {
00568 if (! s)
00569 {
00570 (*current_liboctave_error_handler) ("invalid row selection");
00571 return RowVector ();
00572 }
00573
00574 char c = *s;
00575 if (c == 'f' || c == 'F')
00576 return row (0);
00577 else if (c == 'l' || c == 'L')
00578 return row (rows () - 1);
00579 else
00580 {
00581 (*current_liboctave_error_handler) ("invalid row selection");
00582 return RowVector ();
00583 }
00584 }
00585
00586 ColumnVector
00587 Matrix::column (int i) const
00588 {
00589 int nr = rows ();
00590 if (i < 0 || i >= cols ())
00591 {
00592 (*current_liboctave_error_handler) ("invalid column selection");
00593 return ColumnVector ();
00594 }
00595
00596 ColumnVector retval (nr);
00597 for (int j = 0; j < nr; j++)
00598 retval.xelem (j) = elem (j, i);
00599
00600 return retval;
00601 }
00602
00603 ColumnVector
00604 Matrix::column (char *s) const
00605 {
00606 if (! s)
00607 {
00608 (*current_liboctave_error_handler) ("invalid column selection");
00609 return ColumnVector ();
00610 }
00611
00612 char c = *s;
00613 if (c == 'f' || c == 'F')
00614 return column (0);
00615 else if (c == 'l' || c == 'L')
00616 return column (cols () - 1);
00617 else
00618 {
00619 (*current_liboctave_error_handler) ("invalid column selection");
00620 return ColumnVector ();
00621 }
00622 }
00623
00624 Matrix
00625 Matrix::inverse (void) const
00626 {
00627 int info;
00628 double rcond;
00629 return inverse (info, rcond, 0, 0);
00630 }
00631
00632 Matrix
00633 Matrix::inverse (int& info) const
00634 {
00635 double rcond;
00636 return inverse (info, rcond, 0, 0);
00637 }
00638
00639 Matrix
00640 Matrix::inverse (int& info, double& rcond, int force, int calc_cond) const
00641 {
00642 Matrix retval;
00643
00644 int nr = rows ();
00645 int nc = cols ();
00646
00647 if (nr != nc || nr == 0 || nc == 0)
00648 (*current_liboctave_error_handler) ("inverse requires square matrix");
00649 else
00650 {
00651 Array<int> ipvt (nr);
00652 int *pipvt = ipvt.fortran_vec ();
00653
00654 retval = *this;
00655 double *tmp_data = retval.fortran_vec ();
00656
00657 Array<double> z(1);
00658 int lwork = -1;
00659
00660
00661 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
00662 z.fortran_vec (), lwork, info));
00663
00664 if (f77_exception_encountered)
00665 {
00666 (*current_liboctave_error_handler)
00667 ("unrecoverable error in dgetri");
00668 return retval;
00669 }
00670
00671 lwork = static_cast<int> (z(0));
00672 lwork = (lwork < 2 *nc ? 2*nc : lwork);
00673 z.resize (lwork);
00674 double *pz = z.fortran_vec ();
00675
00676 info = 0;
00677
00678
00679 double anorm = 0;
00680 if (calc_cond)
00681 anorm = retval.abs().sum().row(0).max();
00682
00683 F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, info));
00684
00685 if (f77_exception_encountered)
00686 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
00687 else
00688 {
00689
00690 rcond = 0.0;
00691 if (info != 0)
00692 info = -1;
00693 else if (calc_cond)
00694 {
00695 int dgecon_info = 0;
00696
00697
00698 char job = '1';
00699 Array<int> iz (nc);
00700 int *piz = iz.fortran_vec ();
00701 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
00702 nc, tmp_data, nr, anorm,
00703 rcond, pz, piz, dgecon_info
00704 F77_CHAR_ARG_LEN (1)));
00705
00706 if (f77_exception_encountered)
00707 (*current_liboctave_error_handler)
00708 ("unrecoverable error in dgecon");
00709
00710 if (dgecon_info != 0)
00711 info = -1;
00712 }
00713
00714 if (info == -1 && ! force)
00715 retval = *this;
00716 else
00717 {
00718 int dgetri_info = 0;
00719
00720 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
00721 pz, lwork, dgetri_info));
00722
00723 if (f77_exception_encountered)
00724 (*current_liboctave_error_handler)
00725 ("unrecoverable error in dgetri");
00726
00727 if (dgetri_info != 0)
00728 info = -1;
00729 }
00730 }
00731 }
00732
00733 return retval;
00734 }
00735
00736 Matrix
00737 Matrix::pseudo_inverse (double tol) const
00738 {
00739 SVD result (*this, SVD::economy);
00740
00741 DiagMatrix S = result.singular_values ();
00742 Matrix U = result.left_singular_matrix ();
00743 Matrix V = result.right_singular_matrix ();
00744
00745 ColumnVector sigma = S.diag ();
00746
00747 int r = sigma.length () - 1;
00748 int nr = rows ();
00749 int nc = cols ();
00750
00751 if (tol <= 0.0)
00752 {
00753 if (nr > nc)
00754 tol = nr * sigma.elem (0) * DBL_EPSILON;
00755 else
00756 tol = nc * sigma.elem (0) * DBL_EPSILON;
00757 }
00758
00759 while (r >= 0 && sigma.elem (r) < tol)
00760 r--;
00761
00762 if (r < 0)
00763 return Matrix (nc, nr, 0.0);
00764 else
00765 {
00766 Matrix Ur = U.extract (0, 0, nr-1, r);
00767 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
00768 Matrix Vr = V.extract (0, 0, nc-1, r);
00769 return Vr * D * Ur.transpose ();
00770 }
00771 }
00772
00773 #if defined (HAVE_FFTW3)
00774
00775 ComplexMatrix
00776 Matrix::fourier (void) const
00777 {
00778 size_t nr = rows ();
00779 size_t nc = cols ();
00780
00781 ComplexMatrix retval (nr, nc);
00782
00783 size_t npts, nsamples;
00784
00785 if (nr == 1 || nc == 1)
00786 {
00787 npts = nr > nc ? nr : nc;
00788 nsamples = 1;
00789 }
00790 else
00791 {
00792 npts = nr;
00793 nsamples = nc;
00794 }
00795
00796 const double *in (fortran_vec ());
00797 Complex *out (retval.fortran_vec ());
00798
00799 octave_fftw::fft (in, out, npts, nsamples);
00800
00801 return retval;
00802 }
00803
00804 ComplexMatrix
00805 Matrix::ifourier (void) const
00806 {
00807 size_t nr = rows ();
00808 size_t nc = cols ();
00809
00810 ComplexMatrix retval (nr, nc);
00811
00812 size_t npts, nsamples;
00813
00814 if (nr == 1 || nc == 1)
00815 {
00816 npts = nr > nc ? nr : nc;
00817 nsamples = 1;
00818 }
00819 else
00820 {
00821 npts = nr;
00822 nsamples = nc;
00823 }
00824
00825 ComplexMatrix tmp (*this);
00826 Complex *in (tmp.fortran_vec ());
00827 Complex *out (retval.fortran_vec ());
00828
00829 octave_fftw::ifft (in, out, npts, nsamples);
00830
00831 return retval;
00832 }
00833
00834 ComplexMatrix
00835 Matrix::fourier2d (void) const
00836 {
00837 dim_vector dv(rows (), cols ());
00838
00839 const double *in = fortran_vec ();
00840 ComplexMatrix retval (rows (), cols ());
00841 octave_fftw::fftNd (in, retval.fortran_vec (), 2, dv);
00842
00843 return retval;
00844 }
00845
00846 ComplexMatrix
00847 Matrix::ifourier2d (void) const
00848 {
00849 dim_vector dv(rows (), cols ());
00850
00851 ComplexMatrix retval (*this);
00852 Complex *out (retval.fortran_vec ());
00853
00854 octave_fftw::ifftNd (out, out, 2, dv);
00855
00856 return retval;
00857 }
00858
00859 #else
00860
00861 ComplexMatrix
00862 Matrix::fourier (void) const
00863 {
00864 ComplexMatrix retval;
00865
00866 int nr = rows ();
00867 int nc = cols ();
00868
00869 int npts, nsamples;
00870
00871 if (nr == 1 || nc == 1)
00872 {
00873 npts = nr > nc ? nr : nc;
00874 nsamples = 1;
00875 }
00876 else
00877 {
00878 npts = nr;
00879 nsamples = nc;
00880 }
00881
00882 int nn = 4*npts+15;
00883
00884 Array<Complex> wsave (nn);
00885 Complex *pwsave = wsave.fortran_vec ();
00886
00887 retval = ComplexMatrix (*this);
00888 Complex *tmp_data = retval.fortran_vec ();
00889
00890 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00891
00892 for (int j = 0; j < nsamples; j++)
00893 {
00894 OCTAVE_QUIT;
00895
00896 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
00897 }
00898
00899 return retval;
00900 }
00901
00902 ComplexMatrix
00903 Matrix::ifourier (void) const
00904 {
00905 ComplexMatrix retval;
00906
00907 int nr = rows ();
00908 int nc = cols ();
00909
00910 int npts, nsamples;
00911
00912 if (nr == 1 || nc == 1)
00913 {
00914 npts = nr > nc ? nr : nc;
00915 nsamples = 1;
00916 }
00917 else
00918 {
00919 npts = nr;
00920 nsamples = nc;
00921 }
00922
00923 int nn = 4*npts+15;
00924
00925 Array<Complex> wsave (nn);
00926 Complex *pwsave = wsave.fortran_vec ();
00927
00928 retval = ComplexMatrix (*this);
00929 Complex *tmp_data = retval.fortran_vec ();
00930
00931 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00932
00933 for (int j = 0; j < nsamples; j++)
00934 {
00935 OCTAVE_QUIT;
00936
00937 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
00938 }
00939
00940 for (int j = 0; j < npts*nsamples; j++)
00941 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
00942
00943 return retval;
00944 }
00945
00946 ComplexMatrix
00947 Matrix::fourier2d (void) const
00948 {
00949 ComplexMatrix retval;
00950
00951 int nr = rows ();
00952 int nc = cols ();
00953
00954 int npts, nsamples;
00955
00956 if (nr == 1 || nc == 1)
00957 {
00958 npts = nr > nc ? nr : nc;
00959 nsamples = 1;
00960 }
00961 else
00962 {
00963 npts = nr;
00964 nsamples = nc;
00965 }
00966
00967 int nn = 4*npts+15;
00968
00969 Array<Complex> wsave (nn);
00970 Complex *pwsave = wsave.fortran_vec ();
00971
00972 retval = ComplexMatrix (*this);
00973 Complex *tmp_data = retval.fortran_vec ();
00974
00975 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00976
00977 for (int j = 0; j < nsamples; j++)
00978 {
00979 OCTAVE_QUIT;
00980
00981 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
00982 }
00983
00984 npts = nc;
00985 nsamples = nr;
00986 nn = 4*npts+15;
00987
00988 wsave.resize (nn);
00989 pwsave = wsave.fortran_vec ();
00990
00991 Array<Complex> tmp (npts);
00992 Complex *prow = tmp.fortran_vec ();
00993
00994 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00995
00996 for (int j = 0; j < nsamples; j++)
00997 {
00998 OCTAVE_QUIT;
00999
01000 for (int i = 0; i < npts; i++)
01001 prow[i] = tmp_data[i*nr + j];
01002
01003 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
01004
01005 for (int i = 0; i < npts; i++)
01006 tmp_data[i*nr + j] = prow[i];
01007 }
01008
01009 return retval;
01010 }
01011
01012 ComplexMatrix
01013 Matrix::ifourier2d (void) const
01014 {
01015 ComplexMatrix retval;
01016
01017 int nr = rows ();
01018 int nc = cols ();
01019
01020 int npts, nsamples;
01021
01022 if (nr == 1 || nc == 1)
01023 {
01024 npts = nr > nc ? nr : nc;
01025 nsamples = 1;
01026 }
01027 else
01028 {
01029 npts = nr;
01030 nsamples = nc;
01031 }
01032
01033 int nn = 4*npts+15;
01034
01035 Array<Complex> wsave (nn);
01036 Complex *pwsave = wsave.fortran_vec ();
01037
01038 retval = ComplexMatrix (*this);
01039 Complex *tmp_data = retval.fortran_vec ();
01040
01041 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01042
01043 for (int j = 0; j < nsamples; j++)
01044 {
01045 OCTAVE_QUIT;
01046
01047 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
01048 }
01049
01050 for (int j = 0; j < npts*nsamples; j++)
01051 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
01052
01053 npts = nc;
01054 nsamples = nr;
01055 nn = 4*npts+15;
01056
01057 wsave.resize (nn);
01058 pwsave = wsave.fortran_vec ();
01059
01060 Array<Complex> tmp (npts);
01061 Complex *prow = tmp.fortran_vec ();
01062
01063 F77_FUNC (cffti, CFFTI) (npts, pwsave);
01064
01065 for (int j = 0; j < nsamples; j++)
01066 {
01067 OCTAVE_QUIT;
01068
01069 for (int i = 0; i < npts; i++)
01070 prow[i] = tmp_data[i*nr + j];
01071
01072 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
01073
01074 for (int i = 0; i < npts; i++)
01075 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts);
01076 }
01077
01078 return retval;
01079 }
01080
01081 #endif
01082
01083 DET
01084 Matrix::determinant (void) const
01085 {
01086 int info;
01087 double rcond;
01088 return determinant (info, rcond, 0);
01089 }
01090
01091 DET
01092 Matrix::determinant (int& info) const
01093 {
01094 double rcond;
01095 return determinant (info, rcond, 0);
01096 }
01097
01098 DET
01099 Matrix::determinant (int& info, double& rcond, int calc_cond) const
01100 {
01101 DET retval;
01102
01103 int nr = rows ();
01104 int nc = cols ();
01105
01106 if (nr == 0 || nc == 0)
01107 {
01108 double d[2];
01109 d[0] = 1.0;
01110 d[1] = 0.0;
01111 retval = DET (d);
01112 }
01113 else
01114 {
01115 Array<int> ipvt (nr);
01116 int *pipvt = ipvt.fortran_vec ();
01117
01118 Matrix atmp = *this;
01119 double *tmp_data = atmp.fortran_vec ();
01120
01121 info = 0;
01122
01123
01124 double anorm = 0;
01125 if (calc_cond)
01126 anorm = atmp.abs().sum().row(0).max();
01127
01128 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01129
01130 if (f77_exception_encountered)
01131 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
01132 else
01133 {
01134
01135 rcond = 0.0;
01136 if (info != 0)
01137 {
01138 info = -1;
01139 retval = DET ();
01140 }
01141 else
01142 {
01143 if (calc_cond)
01144 {
01145
01146 char job = '1';
01147 Array<double> z (4 * nc);
01148 double *pz = z.fortran_vec ();
01149 Array<int> iz (nc);
01150 int *piz = iz.fortran_vec ();
01151
01152 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01153 nc, tmp_data, nr, anorm,
01154 rcond, pz, piz, info
01155 F77_CHAR_ARG_LEN (1)));
01156
01157 if (f77_exception_encountered)
01158 (*current_liboctave_error_handler)
01159 ("unrecoverable error in dgecon");
01160 }
01161
01162 if (info != 0)
01163 {
01164 info = -1;
01165 retval = DET ();
01166 }
01167 else
01168 {
01169 double d[2] = { 1., 0.};
01170 for (int i=0; i<nc; i++)
01171 {
01172 if (ipvt(i) != (i+1)) d[0] = -d[0];
01173 d[0] *= atmp(i,i);
01174 if (d[0] == 0.) break;
01175 while (fabs(d[0]) < 1.)
01176 {
01177 d[0] = 10. * d[0];
01178 d[1] = d[1] - 1.0;
01179 }
01180 while (fabs(d[0]) >= 10.)
01181 {
01182 d[0] = 0.1 * d[0];
01183 d[1] = d[1] + 1.0;
01184 }
01185 }
01186 retval = DET (d);
01187 }
01188 }
01189 }
01190 }
01191
01192 return retval;
01193 }
01194
01195 Matrix
01196 Matrix::solve (const Matrix& b) const
01197 {
01198 int info;
01199 double rcond;
01200 return solve (b, info, rcond, 0);
01201 }
01202
01203 Matrix
01204 Matrix::solve (const Matrix& b, int& info) const
01205 {
01206 double rcond;
01207 return solve (b, info, rcond, 0);
01208 }
01209
01210 Matrix
01211 Matrix::solve (const Matrix& b, int& info, double& rcond) const
01212 {
01213 return solve (b, info, rcond, 0);
01214 }
01215
01216 Matrix
01217 Matrix::solve (const Matrix& b, int& info, double& rcond,
01218 solve_singularity_handler sing_handler) const
01219 {
01220 Matrix retval;
01221
01222 int nr = rows ();
01223 int nc = cols ();
01224
01225 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
01226 (*current_liboctave_error_handler)
01227 ("matrix dimension mismatch solution of linear equations");
01228 else
01229 {
01230 info = 0;
01231
01232 Array<int> ipvt (nr);
01233 int *pipvt = ipvt.fortran_vec ();
01234
01235 Matrix atmp = *this;
01236 double *tmp_data = atmp.fortran_vec ();
01237
01238 Array<double> z (4 * nc);
01239 double *pz = z.fortran_vec ();
01240 Array<int> iz (nc);
01241 int *piz = iz.fortran_vec ();
01242
01243
01244 double anorm = atmp.abs().sum().row(0).max();
01245
01246 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01247
01248 if (f77_exception_encountered)
01249 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
01250 else
01251 {
01252
01253 rcond = 0.0;
01254 if (info != 0)
01255 {
01256 info = -2;
01257
01258 if (sing_handler)
01259 sing_handler (rcond);
01260 else
01261 (*current_liboctave_error_handler)
01262 ("matrix singular to machine precision");
01263
01264 }
01265 else
01266 {
01267
01268 char job = '1';
01269 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01270 nc, tmp_data, nr, anorm,
01271 rcond, pz, piz, info
01272 F77_CHAR_ARG_LEN (1)));
01273
01274 if (f77_exception_encountered)
01275 (*current_liboctave_error_handler)
01276 ("unrecoverable error in dgecon");
01277
01278 if (info != 0)
01279 info = -2;
01280
01281 volatile double rcond_plus_one = rcond + 1.0;
01282
01283 if (rcond_plus_one == 1.0 || xisnan (rcond))
01284 {
01285 info = -2;
01286
01287 if (sing_handler)
01288 sing_handler (rcond);
01289 else
01290 (*current_liboctave_error_handler)
01291 ("matrix singular to machine precision, rcond = %g",
01292 rcond);
01293 }
01294 else
01295 {
01296 retval = b;
01297 double *result = retval.fortran_vec ();
01298
01299 int b_nc = b.cols ();
01300
01301 job = 'N';
01302 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
01303 nr, b_nc, tmp_data, nr,
01304 pipvt, result, b.rows(), info
01305 F77_CHAR_ARG_LEN (1)));
01306
01307 if (f77_exception_encountered)
01308 (*current_liboctave_error_handler)
01309 ("unrecoverable error in dgetrs");
01310 }
01311 }
01312 }
01313 }
01314
01315 return retval;
01316 }
01317
01318 ComplexMatrix
01319 Matrix::solve (const ComplexMatrix& b) const
01320 {
01321 ComplexMatrix tmp (*this);
01322 return tmp.solve (b);
01323 }
01324
01325 ComplexMatrix
01326 Matrix::solve (const ComplexMatrix& b, int& info) const
01327 {
01328 ComplexMatrix tmp (*this);
01329 return tmp.solve (b, info);
01330 }
01331
01332 ComplexMatrix
01333 Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
01334 {
01335 ComplexMatrix tmp (*this);
01336 return tmp.solve (b, info, rcond);
01337 }
01338
01339 ComplexMatrix
01340 Matrix::solve (const ComplexMatrix& b, int& info, double& rcond,
01341 solve_singularity_handler sing_handler) const
01342 {
01343 ComplexMatrix tmp (*this);
01344 return tmp.solve (b, info, rcond, sing_handler);
01345 }
01346
01347 ColumnVector
01348 Matrix::solve (const ColumnVector& b) const
01349 {
01350 int info; double rcond;
01351 return solve (b, info, rcond);
01352 }
01353
01354 ColumnVector
01355 Matrix::solve (const ColumnVector& b, int& info) const
01356 {
01357 double rcond;
01358 return solve (b, info, rcond);
01359 }
01360
01361 ColumnVector
01362 Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
01363 {
01364 return solve (b, info, rcond, 0);
01365 }
01366
01367 ColumnVector
01368 Matrix::solve (const ColumnVector& b, int& info, double& rcond,
01369 solve_singularity_handler sing_handler) const
01370 {
01371 ColumnVector retval;
01372
01373 int nr = rows ();
01374 int nc = cols ();
01375
01376 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
01377 (*current_liboctave_error_handler)
01378 ("matrix dimension mismatch solution of linear equations");
01379 else
01380 {
01381 info = 0;
01382
01383 Array<int> ipvt (nr);
01384 int *pipvt = ipvt.fortran_vec ();
01385
01386 Matrix atmp = *this;
01387 double *tmp_data = atmp.fortran_vec ();
01388
01389 Array<double> z (4 * nc);
01390 double *pz = z.fortran_vec ();
01391 Array<int> iz (nc);
01392 int *piz = iz.fortran_vec ();
01393
01394
01395 double anorm = atmp.abs().sum().row(0).max();
01396
01397 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
01398
01399 if (f77_exception_encountered)
01400 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
01401 else
01402 {
01403
01404 rcond = 0.0;
01405 if (info > 0)
01406 {
01407 info = -2;
01408
01409 if (sing_handler)
01410 sing_handler (rcond);
01411 else
01412 (*current_liboctave_error_handler)
01413 ("matrix singular to machine precision");
01414
01415 }
01416 else
01417 {
01418
01419 char job = '1';
01420 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
01421 nc, tmp_data, nr, anorm,
01422 rcond, pz, piz, info
01423 F77_CHAR_ARG_LEN (1)));
01424
01425 if (f77_exception_encountered)
01426 (*current_liboctave_error_handler)
01427 ("unrecoverable error in dgecon");
01428
01429 if (info != 0)
01430 info = -2;
01431
01432 volatile double rcond_plus_one = rcond + 1.0;
01433
01434 if (rcond_plus_one == 1.0 || xisnan (rcond))
01435 {
01436 info = -2;
01437
01438 if (sing_handler)
01439 sing_handler (rcond);
01440 else
01441 (*current_liboctave_error_handler)
01442 ("matrix singular to machine precision, rcond = %g",
01443 rcond);
01444 }
01445 else
01446 {
01447 retval = b;
01448 double *result = retval.fortran_vec ();
01449
01450 job = 'N';
01451 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
01452 nr, 1, tmp_data, nr, pipvt,
01453 result, b.length(), info
01454 F77_CHAR_ARG_LEN (1)));
01455
01456 if (f77_exception_encountered)
01457 (*current_liboctave_error_handler)
01458 ("unrecoverable error in dgetrs");
01459 }
01460 }
01461 }
01462 }
01463
01464 return retval;
01465 }
01466
01467 ComplexColumnVector
01468 Matrix::solve (const ComplexColumnVector& b) const
01469 {
01470 ComplexMatrix tmp (*this);
01471 return tmp.solve (b);
01472 }
01473
01474 ComplexColumnVector
01475 Matrix::solve (const ComplexColumnVector& b, int& info) const
01476 {
01477 ComplexMatrix tmp (*this);
01478 return tmp.solve (b, info);
01479 }
01480
01481 ComplexColumnVector
01482 Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
01483 {
01484 ComplexMatrix tmp (*this);
01485 return tmp.solve (b, info, rcond);
01486 }
01487
01488 ComplexColumnVector
01489 Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond,
01490 solve_singularity_handler sing_handler) const
01491 {
01492 ComplexMatrix tmp (*this);
01493 return tmp.solve (b, info, rcond, sing_handler);
01494 }
01495
01496 Matrix
01497 Matrix::lssolve (const Matrix& b) const
01498 {
01499 int info;
01500 int rank;
01501 return lssolve (b, info, rank);
01502 }
01503
01504 Matrix
01505 Matrix::lssolve (const Matrix& b, int& info) const
01506 {
01507 int rank;
01508 return lssolve (b, info, rank);
01509 }
01510
01511 Matrix
01512 Matrix::lssolve (const Matrix& b, int& info, int& rank) const
01513 {
01514 Matrix retval;
01515
01516 int nrhs = b.cols ();
01517
01518 int m = rows ();
01519 int n = cols ();
01520
01521 if (m == 0 || n == 0 || m != b.rows ())
01522 (*current_liboctave_error_handler)
01523 ("matrix dimension mismatch in solution of least squares problem");
01524 else
01525 {
01526 Matrix atmp = *this;
01527 double *tmp_data = atmp.fortran_vec ();
01528
01529 int nrr = m > n ? m : n;
01530 Matrix result (nrr, nrhs, 0.0);
01531
01532 for (int j = 0; j < nrhs; j++)
01533 for (int i = 0; i < m; i++)
01534 result.elem (i, j) = b.elem (i, j);
01535
01536 double *presult = result.fortran_vec ();
01537
01538 int len_s = m < n ? m : n;
01539 Array<double> s (len_s);
01540 double *ps = s.fortran_vec ();
01541
01542 double rcond = -1.0;
01543
01544
01545
01546 int lwork = -1;
01547
01548 Array<double> work (1);
01549
01550 F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps,
01551 rcond, rank, work.fortran_vec (),
01552 lwork, info));
01553
01554 if (f77_exception_encountered)
01555 (*current_liboctave_error_handler) ("unrecoverable error in dgelss");
01556 else
01557 {
01558 lwork = static_cast<int> (work(0));
01559 work.resize (lwork);
01560
01561 F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult,
01562 nrr, ps, rcond, rank,
01563 work.fortran_vec (), lwork, info));
01564
01565 if (f77_exception_encountered)
01566 (*current_liboctave_error_handler)
01567 ("unrecoverable error in dgelss");
01568 else
01569 {
01570 retval.resize (n, nrhs);
01571 for (int j = 0; j < nrhs; j++)
01572 for (int i = 0; i < n; i++)
01573 retval.elem (i, j) = result.elem (i, j);
01574 }
01575 }
01576 }
01577
01578 return retval;
01579 }
01580
01581 ComplexMatrix
01582 Matrix::lssolve (const ComplexMatrix& b) const
01583 {
01584 ComplexMatrix tmp (*this);
01585 int info;
01586 int rank;
01587 return tmp.lssolve (b, info, rank);
01588 }
01589
01590 ComplexMatrix
01591 Matrix::lssolve (const ComplexMatrix& b, int& info) const
01592 {
01593 ComplexMatrix tmp (*this);
01594 int rank;
01595 return tmp.lssolve (b, info, rank);
01596 }
01597
01598 ComplexMatrix
01599 Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
01600 {
01601 ComplexMatrix tmp (*this);
01602 return tmp.lssolve (b, info, rank);
01603 }
01604
01605 ColumnVector
01606 Matrix::lssolve (const ColumnVector& b) const
01607 {
01608 int info;
01609 int rank;
01610 return lssolve (b, info, rank);
01611 }
01612
01613 ColumnVector
01614 Matrix::lssolve (const ColumnVector& b, int& info) const
01615 {
01616 int rank;
01617 return lssolve (b, info, rank);
01618 }
01619
01620 ColumnVector
01621 Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const
01622 {
01623 ColumnVector retval;
01624
01625 int nrhs = 1;
01626
01627 int m = rows ();
01628 int n = cols ();
01629
01630 if (m == 0 || n == 0 || m != b.length ())
01631 (*current_liboctave_error_handler)
01632 ("matrix dimension mismatch in solution of least squares problem");
01633 else
01634 {
01635 Matrix atmp = *this;
01636 double *tmp_data = atmp.fortran_vec ();
01637
01638 int nrr = m > n ? m : n;
01639 ColumnVector result (nrr);
01640
01641 for (int i = 0; i < m; i++)
01642 result.elem (i) = b.elem (i);
01643
01644 double *presult = result.fortran_vec ();
01645
01646 int len_s = m < n ? m : n;
01647 Array<double> s (len_s);
01648 double *ps = s.fortran_vec ();
01649
01650 double rcond = -1.0;
01651
01652
01653
01654 int lwork = -1;
01655
01656 Array<double> work (1);
01657
01658 F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps,
01659 rcond, rank, work.fortran_vec (),
01660 lwork, info));
01661
01662 if (f77_exception_encountered)
01663 (*current_liboctave_error_handler) ("unrecoverable error in dgelss");
01664 else
01665 {
01666 lwork = static_cast<int> (work(0));
01667 work.resize (lwork);
01668
01669 F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult,
01670 nrr, ps, rcond, rank,
01671 work.fortran_vec (), lwork, info));
01672
01673 if (f77_exception_encountered)
01674 (*current_liboctave_error_handler)
01675 ("unrecoverable error in dgelss");
01676 else
01677 {
01678 retval.resize (n);
01679 for (int i = 0; i < n; i++)
01680 retval.elem (i) = result.elem (i);
01681 }
01682 }
01683 }
01684
01685 return retval;
01686 }
01687
01688 ComplexColumnVector
01689 Matrix::lssolve (const ComplexColumnVector& b) const
01690 {
01691 ComplexMatrix tmp (*this);
01692 return tmp.lssolve (b);
01693 }
01694
01695 ComplexColumnVector
01696 Matrix::lssolve (const ComplexColumnVector& b, int& info) const
01697 {
01698 ComplexMatrix tmp (*this);
01699 return tmp.lssolve (b, info);
01700 }
01701
01702 ComplexColumnVector
01703 Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
01704 {
01705 ComplexMatrix tmp (*this);
01706 return tmp.lssolve (b, info, rank);
01707 }
01708
01709
01710
01711 static double padec [] =
01712 {
01713 5.0000000000000000e-1,
01714 1.1666666666666667e-1,
01715 1.6666666666666667e-2,
01716 1.6025641025641026e-3,
01717 1.0683760683760684e-4,
01718 4.8562548562548563e-6,
01719 1.3875013875013875e-7,
01720 1.9270852604185938e-9,
01721 };
01722
01723 Matrix
01724 Matrix::expm (void) const
01725 {
01726 Matrix retval;
01727
01728 Matrix m = *this;
01729
01730 int nc = columns ();
01731
01732
01733
01734
01735
01736 volatile double trshift = 0.0;
01737
01738 for (int i = 0; i < nc; i++)
01739 trshift += m.elem (i, i);
01740
01741 trshift /= nc;
01742
01743 if (trshift > 0.0)
01744 {
01745 for (int i = 0; i < nc; i++)
01746 m.elem (i, i) -= trshift;
01747 }
01748
01749
01750
01751
01752 double *p_m = m.fortran_vec ();
01753
01754 int info, ilo, ihi, ilos, ihis;
01755 Array<double> dpermute (nc);
01756 Array<double> dscale (nc);
01757
01758
01759 char job = 'P';
01760 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
01761 nc, p_m, nc, ilo, ihi,
01762 dpermute.fortran_vec (), info
01763 F77_CHAR_ARG_LEN (1)));
01764
01765
01766 job = 'S';
01767 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
01768 nc, p_m, nc, ilos, ihis,
01769 dscale.fortran_vec (), info
01770 F77_CHAR_ARG_LEN (1)));
01771
01772 if (f77_exception_encountered)
01773 {
01774 (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
01775 return retval;
01776 }
01777
01778
01779
01780 ColumnVector work(nc);
01781 double inf_norm;
01782
01783 F77_XFCN (xdlange, XDLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
01784 nc, nc, m.fortran_vec (), nc,
01785 work.fortran_vec (), inf_norm
01786 F77_CHAR_ARG_LEN (1)));
01787
01788 if (f77_exception_encountered)
01789 {
01790 (*current_liboctave_error_handler) ("unrecoverable error in dlange");
01791 return retval;
01792 }
01793
01794 int sqpow = (int) (inf_norm > 0.0
01795 ? (1.0 + log (inf_norm) / log (2.0))
01796 : 0.0);
01797
01798
01799
01800 if (sqpow < 0)
01801 sqpow = 0;
01802
01803 if (sqpow > 0)
01804 {
01805 double scale_factor = 1.0;
01806 for (int i = 0; i < sqpow; i++)
01807 scale_factor *= 2.0;
01808
01809 m = m / scale_factor;
01810 }
01811
01812
01813
01814 Matrix npp (nc, nc, 0.0);
01815 Matrix dpp = npp;
01816
01817
01818
01819 int minus_one_j = -1;
01820 for (int j = 7; j >= 0; j--)
01821 {
01822 npp = m * npp + padec[j] * m;
01823 dpp = m * dpp + (minus_one_j * padec[j]) * m;
01824 minus_one_j *= -1;
01825 }
01826
01827
01828
01829 dpp = -dpp;
01830 for (int j = 0; j < nc; j++)
01831 {
01832 npp.elem (j, j) += 1.0;
01833 dpp.elem (j, j) += 1.0;
01834 }
01835
01836
01837
01838 retval = dpp.solve (npp, info);
01839
01840
01841
01842 while (sqpow)
01843 {
01844 retval = retval * retval;
01845 sqpow--;
01846 }
01847
01848
01849
01850 for (int i = 0; i < nc; i++)
01851 for (int j = 0; j < nc; j++)
01852 retval(i,j) *= dscale(i) / dscale(j);
01853
01854 OCTAVE_QUIT;
01855
01856
01857 Array<int> iperm (nc);
01858 for (int i = 0; i < nc; i++)
01859 iperm(i) = i;
01860
01861
01862 for (int i = 0; i < (ilo-1); i++)
01863 {
01864 int swapidx = static_cast<int> (dpermute(i)) - 1;
01865 int tmp = iperm(i);
01866 iperm(i) = iperm (swapidx);
01867 iperm(swapidx) = tmp;
01868 }
01869
01870
01871 for (int i = nc - 1; i >= ihi; i--)
01872 {
01873 int swapidx = static_cast<int> (dpermute(i)) - 1;
01874 int tmp = iperm(i);
01875 iperm(i) = iperm(swapidx);
01876 iperm(swapidx) = tmp;
01877 }
01878
01879
01880 Array<int> invpvec (nc);
01881 for (int i = 0; i < nc; i++)
01882 invpvec(iperm(i)) = i;
01883
01884 OCTAVE_QUIT;
01885
01886 Matrix tmpMat = retval;
01887 for (int i = 0; i < nc; i++)
01888 for (int j = 0; j < nc; j++)
01889 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
01890
01891
01892
01893 if (trshift > 0.0)
01894 retval = exp (trshift) * retval;
01895
01896 return retval;
01897 }
01898
01899 Matrix&
01900 Matrix::operator += (const DiagMatrix& a)
01901 {
01902 int nr = rows ();
01903 int nc = cols ();
01904
01905 int a_nr = a.rows ();
01906 int a_nc = a.cols ();
01907
01908 if (nr != a_nr || nc != a_nc)
01909 {
01910 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
01911 return *this;
01912 }
01913
01914 for (int i = 0; i < a.length (); i++)
01915 elem (i, i) += a.elem (i, i);
01916
01917 return *this;
01918 }
01919
01920 Matrix&
01921 Matrix::operator -= (const DiagMatrix& a)
01922 {
01923 int nr = rows ();
01924 int nc = cols ();
01925
01926 int a_nr = a.rows ();
01927 int a_nc = a.cols ();
01928
01929 if (nr != a_nr || nc != a_nc)
01930 {
01931 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
01932 return *this;
01933 }
01934
01935 for (int i = 0; i < a.length (); i++)
01936 elem (i, i) -= a.elem (i, i);
01937
01938 return *this;
01939 }
01940
01941
01942
01943 boolMatrix
01944 Matrix::operator ! (void) const
01945 {
01946 int nr = rows ();
01947 int nc = cols ();
01948
01949 boolMatrix b (nr, nc);
01950
01951 for (int j = 0; j < nc; j++)
01952 for (int i = 0; i < nr; i++)
01953 b.elem (i, j) = ! elem (i, j);
01954
01955 return b;
01956 }
01957
01958
01959
01960 Matrix
01961 operator * (const ColumnVector& v, const RowVector& a)
01962 {
01963 Matrix retval;
01964
01965 int len = v.length ();
01966
01967 if (len != 0)
01968 {
01969 int a_len = a.length ();
01970
01971 retval.resize (len, a_len);
01972 double *c = retval.fortran_vec ();
01973
01974 F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
01975 F77_CONST_CHAR_ARG2 ("N", 1),
01976 len, a_len, 1, 1.0, v.data (), len,
01977 a.data (), 1, 0.0, c, len
01978 F77_CHAR_ARG_LEN (1)
01979 F77_CHAR_ARG_LEN (1)));
01980
01981 if (f77_exception_encountered)
01982 (*current_liboctave_error_handler)
01983 ("unrecoverable error in dgemm");
01984 }
01985
01986 return retval;
01987 }
01988
01989
01990
01991 Matrix
01992 Matrix::map (d_d_Mapper f) const
01993 {
01994 Matrix b (*this);
01995 return b.apply (f);
01996 }
01997
01998 boolMatrix
01999 Matrix::map (b_d_Mapper f) const
02000 {
02001 int nr = rows ();
02002 int nc = cols ();
02003
02004 boolMatrix retval (nr, nc);
02005
02006 for (int j = 0; j < nc; j++)
02007 for (int i = 0; i < nr; i++)
02008 retval(i,j) = f (elem(i,j));
02009
02010 return retval;
02011 }
02012
02013 Matrix&
02014 Matrix::apply (d_d_Mapper f)
02015 {
02016 double *d = fortran_vec ();
02017
02018 for (int i = 0; i < length (); i++)
02019 d[i] = f (d[i]);
02020
02021 return *this;
02022 }
02023
02024 bool
02025 Matrix::any_element_is_negative (bool neg_zero) const
02026 {
02027 int nel = nelem ();
02028
02029 if (neg_zero)
02030 {
02031 for (int i = 0; i < nel; i++)
02032 if (lo_ieee_signbit (elem (i)))
02033 return true;
02034 }
02035 else
02036 {
02037 for (int i = 0; i < nel; i++)
02038 if (elem (i) < 0)
02039 return true;
02040 }
02041
02042 return false;
02043 }
02044
02045
02046 bool
02047 Matrix::any_element_is_inf_or_nan (void) const
02048 {
02049 int nel = nelem ();
02050
02051 for (int i = 0; i < nel; i++)
02052 {
02053 double val = elem (i);
02054 if (xisinf (val) || xisnan (val))
02055 return true;
02056 }
02057
02058 return false;
02059 }
02060
02061 bool
02062 Matrix::all_elements_are_int_or_inf_or_nan (void) const
02063 {
02064 int nel = nelem ();
02065
02066 for (int i = 0; i < nel; i++)
02067 {
02068 double val = elem (i);
02069 if (xisnan (val) || D_NINT (val) == val)
02070 continue;
02071 else
02072 return false;
02073 }
02074
02075 return true;
02076 }
02077
02078
02079
02080
02081 bool
02082 Matrix::all_integers (double& max_val, double& min_val) const
02083 {
02084 int nel = nelem ();
02085
02086 if (nel > 0)
02087 {
02088 max_val = elem (0);
02089 min_val = elem (0);
02090 }
02091 else
02092 return false;
02093
02094 for (int i = 0; i < nel; i++)
02095 {
02096 double val = elem (i);
02097
02098 if (val > max_val)
02099 max_val = val;
02100
02101 if (val < min_val)
02102 min_val = val;
02103
02104 if (D_NINT (val) != val)
02105 return false;
02106 }
02107
02108 return true;
02109 }
02110
02111 bool
02112 Matrix::too_large_for_float (void) const
02113 {
02114 int nel = nelem ();
02115
02116 for (int i = 0; i < nel; i++)
02117 {
02118 double val = elem (i);
02119
02120 if (val > FLT_MAX || val < FLT_MIN)
02121 return true;
02122 }
02123
02124 return false;
02125 }
02126
02127
02128
02129
02130 boolMatrix
02131 Matrix::all (int dim) const
02132 {
02133 MX_ALL_OP (dim);
02134 }
02135
02136 boolMatrix
02137 Matrix::any (int dim) const
02138 {
02139 MX_ANY_OP (dim);
02140 }
02141
02142 Matrix
02143 Matrix::cumprod (int dim) const
02144 {
02145 MX_CUMULATIVE_OP (Matrix, double, *=);
02146 }
02147
02148 Matrix
02149 Matrix::cumsum (int dim) const
02150 {
02151 MX_CUMULATIVE_OP (Matrix, double, +=);
02152 }
02153
02154 Matrix
02155 Matrix::prod (int dim) const
02156 {
02157 MX_REDUCTION_OP (Matrix, *=, 1.0, 1.0);
02158 }
02159
02160 Matrix
02161 Matrix::sum (int dim) const
02162 {
02163 MX_REDUCTION_OP (Matrix, +=, 0.0, 0.0);
02164 }
02165
02166 Matrix
02167 Matrix::sumsq (int dim) const
02168 {
02169 #define ROW_EXPR \
02170 double d = elem (i, j); \
02171 retval.elem (i, 0) += d * d
02172
02173 #define COL_EXPR \
02174 double d = elem (i, j); \
02175 retval.elem (0, j) += d * d
02176
02177 MX_BASE_REDUCTION_OP (Matrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
02178
02179 #undef ROW_EXPR
02180 #undef COL_EXPR
02181 }
02182
02183 Matrix
02184 Matrix::abs (void) const
02185 {
02186 int nr = rows ();
02187 int nc = cols ();
02188
02189 Matrix retval (nr, nc);
02190
02191 for (int j = 0; j < nc; j++)
02192 for (int i = 0; i < nr; i++)
02193 retval (i, j) = fabs (elem (i, j));
02194
02195 return retval;
02196 }
02197
02198 ColumnVector
02199 Matrix::diag (void) const
02200 {
02201 return diag (0);
02202 }
02203
02204 ColumnVector
02205 Matrix::diag (int k) const
02206 {
02207 int nnr = rows ();
02208 int nnc = cols ();
02209 if (k > 0)
02210 nnc -= k;
02211 else if (k < 0)
02212 nnr += k;
02213
02214 ColumnVector d;
02215
02216 if (nnr > 0 && nnc > 0)
02217 {
02218 int ndiag = (nnr < nnc) ? nnr : nnc;
02219
02220 d.resize (ndiag);
02221
02222 if (k > 0)
02223 {
02224 for (int i = 0; i < ndiag; i++)
02225 d.elem (i) = elem (i, i+k);
02226 }
02227 else if (k < 0)
02228 {
02229 for (int i = 0; i < ndiag; i++)
02230 d.elem (i) = elem (i-k, i);
02231 }
02232 else
02233 {
02234 for (int i = 0; i < ndiag; i++)
02235 d.elem (i) = elem (i, i);
02236 }
02237 }
02238 else
02239 (*current_liboctave_error_handler)
02240 ("diag: requested diagonal out of range");
02241
02242 return d;
02243 }
02244
02245 ColumnVector
02246 Matrix::row_min (void) const
02247 {
02248 Array<int> dummy_idx;
02249 return row_min (dummy_idx);
02250 }
02251
02252 ColumnVector
02253 Matrix::row_min (Array<int>& idx_arg) const
02254 {
02255 ColumnVector result;
02256
02257 int nr = rows ();
02258 int nc = cols ();
02259
02260 if (nr > 0 && nc > 0)
02261 {
02262 result.resize (nr);
02263 idx_arg.resize (nr);
02264
02265 for (int i = 0; i < nr; i++)
02266 {
02267 int idx_j;
02268
02269 double tmp_min = octave_NaN;
02270
02271 for (idx_j = 0; idx_j < nc; idx_j++)
02272 {
02273 tmp_min = elem (i, idx_j);
02274
02275 if (! octave_is_NaN_or_NA (tmp_min))
02276 break;
02277 }
02278
02279 for (int j = idx_j+1; j < nc; j++)
02280 {
02281 double tmp = elem (i, j);
02282
02283 if (octave_is_NaN_or_NA (tmp))
02284 continue;
02285 else if (tmp < tmp_min)
02286 {
02287 idx_j = j;
02288 tmp_min = tmp;
02289 }
02290 }
02291
02292 result.elem (i) = tmp_min;
02293 idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j;
02294 }
02295 }
02296
02297 return result;
02298 }
02299
02300 ColumnVector
02301 Matrix::row_max (void) const
02302 {
02303 Array<int> dummy_idx;
02304 return row_max (dummy_idx);
02305 }
02306
02307 ColumnVector
02308 Matrix::row_max (Array<int>& idx_arg) const
02309 {
02310 ColumnVector result;
02311
02312 int nr = rows ();
02313 int nc = cols ();
02314
02315 if (nr > 0 && nc > 0)
02316 {
02317 result.resize (nr);
02318 idx_arg.resize (nr);
02319
02320 for (int i = 0; i < nr; i++)
02321 {
02322 int idx_j;
02323
02324 double tmp_max = octave_NaN;
02325
02326 for (idx_j = 0; idx_j < nc; idx_j++)
02327 {
02328 tmp_max = elem (i, idx_j);
02329
02330 if (! octave_is_NaN_or_NA (tmp_max))
02331 break;
02332 }
02333
02334 for (int j = idx_j+1; j < nc; j++)
02335 {
02336 double tmp = elem (i, j);
02337
02338 if (octave_is_NaN_or_NA (tmp))
02339 continue;
02340 else if (tmp > tmp_max)
02341 {
02342 idx_j = j;
02343 tmp_max = tmp;
02344 }
02345 }
02346
02347 result.elem (i) = tmp_max;
02348 idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j;
02349 }
02350 }
02351
02352 return result;
02353 }
02354
02355 RowVector
02356 Matrix::column_min (void) const
02357 {
02358 Array<int> dummy_idx;
02359 return column_min (dummy_idx);
02360 }
02361
02362 RowVector
02363 Matrix::column_min (Array<int>& idx_arg) const
02364 {
02365 RowVector result;
02366
02367 int nr = rows ();
02368 int nc = cols ();
02369
02370 if (nr > 0 && nc > 0)
02371 {
02372 result.resize (nc);
02373 idx_arg.resize (nc);
02374
02375 for (int j = 0; j < nc; j++)
02376 {
02377 int idx_i;
02378
02379 double tmp_min = octave_NaN;
02380
02381 for (idx_i = 0; idx_i < nr; idx_i++)
02382 {
02383 tmp_min = elem (idx_i, j);
02384
02385 if (! octave_is_NaN_or_NA (tmp_min))
02386 break;
02387 }
02388
02389 for (int i = idx_i+1; i < nr; i++)
02390 {
02391 double tmp = elem (i, j);
02392
02393 if (octave_is_NaN_or_NA (tmp))
02394 continue;
02395 else if (tmp < tmp_min)
02396 {
02397 idx_i = i;
02398 tmp_min = tmp;
02399 }
02400 }
02401
02402 result.elem (j) = tmp_min;
02403 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_i;
02404 }
02405 }
02406
02407 return result;
02408 }
02409
02410 RowVector
02411 Matrix::column_max (void) const
02412 {
02413 Array<int> dummy_idx;
02414 return column_max (dummy_idx);
02415 }
02416
02417 RowVector
02418 Matrix::column_max (Array<int>& idx_arg) const
02419 {
02420 RowVector result;
02421
02422 int nr = rows ();
02423 int nc = cols ();
02424
02425 if (nr > 0 && nc > 0)
02426 {
02427 result.resize (nc);
02428 idx_arg.resize (nc);
02429
02430 for (int j = 0; j < nc; j++)
02431 {
02432 int idx_i;
02433
02434 double tmp_max = octave_NaN;
02435
02436 for (idx_i = 0; idx_i < nr; idx_i++)
02437 {
02438 tmp_max = elem (idx_i, j);
02439
02440 if (! octave_is_NaN_or_NA (tmp_max))
02441 break;
02442 }
02443
02444 for (int i = idx_i+1; i < nr; i++)
02445 {
02446 double tmp = elem (i, j);
02447
02448 if (octave_is_NaN_or_NA (tmp))
02449 continue;
02450 else if (tmp > tmp_max)
02451 {
02452 idx_i = i;
02453 tmp_max = tmp;
02454 }
02455 }
02456
02457 result.elem (j) = tmp_max;
02458 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_i;
02459 }
02460 }
02461
02462 return result;
02463 }
02464
02465 std::ostream&
02466 operator << (std::ostream& os, const Matrix& a)
02467 {
02468 for (int i = 0; i < a.rows (); i++)
02469 {
02470 for (int j = 0; j < a.cols (); j++)
02471 {
02472 os << " ";
02473 octave_write_double (os, a.elem (i, j));
02474 }
02475 os << "\n";
02476 }
02477 return os;
02478 }
02479
02480 std::istream&
02481 operator >> (std::istream& is, Matrix& a)
02482 {
02483 int nr = a.rows ();
02484 int nc = a.cols ();
02485
02486 if (nr < 1 || nc < 1)
02487 is.clear (std::ios::badbit);
02488 else
02489 {
02490 double tmp;
02491 for (int i = 0; i < nr; i++)
02492 for (int j = 0; j < nc; j++)
02493 {
02494 tmp = octave_read_double (is);
02495 if (is)
02496 a.elem (i, j) = tmp;
02497 else
02498 goto done;
02499 }
02500 }
02501
02502 done:
02503
02504 return is;
02505 }
02506
02507 Matrix
02508 Givens (double x, double y)
02509 {
02510 double cc, s, temp_r;
02511
02512 F77_FUNC (dlartg, DLARTG) (x, y, cc, s, temp_r);
02513
02514 Matrix g (2, 2);
02515
02516 g.elem (0, 0) = cc;
02517 g.elem (1, 1) = cc;
02518 g.elem (0, 1) = s;
02519 g.elem (1, 0) = -s;
02520
02521 return g;
02522 }
02523
02524 Matrix
02525 Sylvester (const Matrix& a, const Matrix& b, const Matrix& c)
02526 {
02527 Matrix retval;
02528
02529
02530
02531
02532
02533
02534 SCHUR as (a, "U");
02535 SCHUR bs (b, "U");
02536
02537
02538
02539 Matrix ua = as.unitary_matrix ();
02540 Matrix sch_a = as.schur_matrix ();
02541
02542 Matrix ub = bs.unitary_matrix ();
02543 Matrix sch_b = bs.schur_matrix ();
02544
02545 Matrix cx = ua.transpose () * c * ub;
02546
02547
02548
02549
02550 int a_nr = a.rows ();
02551 int b_nr = b.rows ();
02552
02553 double scale;
02554 int info;
02555
02556 double *pa = sch_a.fortran_vec ();
02557 double *pb = sch_b.fortran_vec ();
02558 double *px = cx.fortran_vec ();
02559
02560 F77_XFCN (dtrsyl, DTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
02561 F77_CONST_CHAR_ARG2 ("N", 1),
02562 1, a_nr, b_nr, pa, a_nr, pb,
02563 b_nr, px, a_nr, scale, info
02564 F77_CHAR_ARG_LEN (1)
02565 F77_CHAR_ARG_LEN (1)));
02566
02567
02568 if (f77_exception_encountered)
02569 (*current_liboctave_error_handler) ("unrecoverable error in dtrsyl");
02570 else
02571 {
02572
02573
02574 retval = -ua*cx*ub.transpose ();
02575 }
02576
02577 return retval;
02578 }
02579
02580
02581
02582 Matrix
02583 operator * (const Matrix& m, const Matrix& a)
02584 {
02585 Matrix retval;
02586
02587 int nr = m.rows ();
02588 int nc = m.cols ();
02589
02590 int a_nr = a.rows ();
02591 int a_nc = a.cols ();
02592
02593 if (nc != a_nr)
02594 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
02595 else
02596 {
02597 if (nr == 0 || nc == 0 || a_nc == 0)
02598 retval.resize (nr, a_nc, 0.0);
02599 else
02600 {
02601 int ld = nr;
02602 int lda = a_nr;
02603
02604 retval.resize (nr, a_nc);
02605 double *c = retval.fortran_vec ();
02606
02607 F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
02608 F77_CONST_CHAR_ARG2 ("N", 1),
02609 nr, a_nc, nc, 1.0, m.data (),
02610 ld, a.data (), lda, 0.0, c, nr
02611 F77_CHAR_ARG_LEN (1)
02612 F77_CHAR_ARG_LEN (1)));
02613
02614 if (f77_exception_encountered)
02615 (*current_liboctave_error_handler)
02616 ("unrecoverable error in dgemm");
02617 }
02618 }
02619
02620 return retval;
02621 }
02622
02623
02624
02625
02626 #define EMPTY_RETURN_CHECK(T) \
02627 if (nr == 0 || nc == 0) \
02628 return T (nr, nc);
02629
02630 Matrix
02631 min (double d, const Matrix& m)
02632 {
02633 int nr = m.rows ();
02634 int nc = m.columns ();
02635
02636 EMPTY_RETURN_CHECK (Matrix);
02637
02638 Matrix result (nr, nc);
02639
02640 for (int j = 0; j < nc; j++)
02641 for (int i = 0; i < nr; i++)
02642 {
02643 OCTAVE_QUIT;
02644 result (i, j) = xmin (d, m (i, j));
02645 }
02646
02647 return result;
02648 }
02649
02650 Matrix
02651 min (const Matrix& m, double d)
02652 {
02653 int nr = m.rows ();
02654 int nc = m.columns ();
02655
02656 EMPTY_RETURN_CHECK (Matrix);
02657
02658 Matrix result (nr, nc);
02659
02660 for (int j = 0; j < nc; j++)
02661 for (int i = 0; i < nr; i++)
02662 {
02663 OCTAVE_QUIT;
02664 result (i, j) = xmin (m (i, j), d);
02665 }
02666
02667 return result;
02668 }
02669
02670 Matrix
02671 min (const Matrix& a, const Matrix& b)
02672 {
02673 int nr = a.rows ();
02674 int nc = a.columns ();
02675
02676 if (nr != b.rows () || nc != b.columns ())
02677 {
02678 (*current_liboctave_error_handler)
02679 ("two-arg min expecting args of same size");
02680 return Matrix ();
02681 }
02682
02683 EMPTY_RETURN_CHECK (Matrix);
02684
02685 Matrix result (nr, nc);
02686
02687 for (int j = 0; j < nc; j++)
02688 for (int i = 0; i < nr; i++)
02689 {
02690 OCTAVE_QUIT;
02691 result (i, j) = xmin (a (i, j), b (i, j));
02692 }
02693
02694 return result;
02695 }
02696
02697 Matrix
02698 max (double d, const Matrix& m)
02699 {
02700 int nr = m.rows ();
02701 int nc = m.columns ();
02702
02703 EMPTY_RETURN_CHECK (Matrix);
02704
02705 Matrix result (nr, nc);
02706
02707 for (int j = 0; j < nc; j++)
02708 for (int i = 0; i < nr; i++)
02709 {
02710 OCTAVE_QUIT;
02711 result (i, j) = xmax (d, m (i, j));
02712 }
02713
02714 return result;
02715 }
02716
02717 Matrix
02718 max (const Matrix& m, double d)
02719 {
02720 int nr = m.rows ();
02721 int nc = m.columns ();
02722
02723 EMPTY_RETURN_CHECK (Matrix);
02724
02725 Matrix result (nr, nc);
02726
02727 for (int j = 0; j < nc; j++)
02728 for (int i = 0; i < nr; i++)
02729 {
02730 OCTAVE_QUIT;
02731 result (i, j) = xmax (m (i, j), d);
02732 }
02733
02734 return result;
02735 }
02736
02737 Matrix
02738 max (const Matrix& a, const Matrix& b)
02739 {
02740 int nr = a.rows ();
02741 int nc = a.columns ();
02742
02743 if (nr != b.rows () || nc != b.columns ())
02744 {
02745 (*current_liboctave_error_handler)
02746 ("two-arg max expecting args of same size");
02747 return Matrix ();
02748 }
02749
02750 EMPTY_RETURN_CHECK (Matrix);
02751
02752 Matrix result (nr, nc);
02753
02754 for (int j = 0; j < nc; j++)
02755 for (int i = 0; i < nr; i++)
02756 {
02757 OCTAVE_QUIT;
02758 result (i, j) = xmax (a (i, j), b (i, j));
02759 }
02760
02761 return result;
02762 }
02763
02764 MS_CMP_OPS(Matrix, , double, )
02765 MS_BOOL_OPS(Matrix, double, 0.0)
02766
02767 SM_CMP_OPS(double, , Matrix, )
02768 SM_BOOL_OPS(double, Matrix, 0.0)
02769
02770 MM_CMP_OPS(Matrix, , Matrix, )
02771 MM_BOOL_OPS(Matrix, Matrix, 0.0)
02772
02773
02774
02775
02776
02777