00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <iostream>
00033
00034 #include "Array-util.h"
00035 #include "lo-error.h"
00036 #include "mx-base.h"
00037 #include "mx-inlines.cc"
00038 #include "oct-cmplx.h"
00039
00040
00041
00042 bool
00043 DiagMatrix::operator == (const DiagMatrix& a) const
00044 {
00045 if (rows () != a.rows () || cols () != a.cols ())
00046 return 0;
00047
00048 return mx_inline_equal (data (), a.data (), length ());
00049 }
00050
00051 bool
00052 DiagMatrix::operator != (const DiagMatrix& a) const
00053 {
00054 return !(*this == a);
00055 }
00056
00057 DiagMatrix&
00058 DiagMatrix::fill (double val)
00059 {
00060 for (int i = 0; i < length (); i++)
00061 elem (i, i) = val;
00062 return *this;
00063 }
00064
00065 DiagMatrix&
00066 DiagMatrix::fill (double val, int beg, int end)
00067 {
00068 if (beg < 0 || end >= length () || end < beg)
00069 {
00070 (*current_liboctave_error_handler) ("range error for fill");
00071 return *this;
00072 }
00073
00074 for (int i = beg; i <= end; i++)
00075 elem (i, i) = val;
00076
00077 return *this;
00078 }
00079
00080 DiagMatrix&
00081 DiagMatrix::fill (const ColumnVector& a)
00082 {
00083 int len = length ();
00084 if (a.length () != len)
00085 {
00086 (*current_liboctave_error_handler) ("range error for fill");
00087 return *this;
00088 }
00089
00090 for (int i = 0; i < len; i++)
00091 elem (i, i) = a.elem (i);
00092
00093 return *this;
00094 }
00095
00096 DiagMatrix&
00097 DiagMatrix::fill (const RowVector& a)
00098 {
00099 int len = length ();
00100 if (a.length () != len)
00101 {
00102 (*current_liboctave_error_handler) ("range error for fill");
00103 return *this;
00104 }
00105
00106 for (int i = 0; i < len; i++)
00107 elem (i, i) = a.elem (i);
00108
00109 return *this;
00110 }
00111
00112 DiagMatrix&
00113 DiagMatrix::fill (const ColumnVector& a, int beg)
00114 {
00115 int a_len = a.length ();
00116 if (beg < 0 || beg + a_len >= length ())
00117 {
00118 (*current_liboctave_error_handler) ("range error for fill");
00119 return *this;
00120 }
00121
00122 for (int i = 0; i < a_len; i++)
00123 elem (i+beg, i+beg) = a.elem (i);
00124
00125 return *this;
00126 }
00127
00128 DiagMatrix&
00129 DiagMatrix::fill (const RowVector& a, int beg)
00130 {
00131 int a_len = a.length ();
00132 if (beg < 0 || beg + a_len >= length ())
00133 {
00134 (*current_liboctave_error_handler) ("range error for fill");
00135 return *this;
00136 }
00137
00138 for (int i = 0; i < a_len; i++)
00139 elem (i+beg, i+beg) = a.elem (i);
00140
00141 return *this;
00142 }
00143
00144 DiagMatrix
00145 DiagMatrix::transpose (void) const
00146 {
00147 return DiagMatrix (mx_inline_dup (data (), length ()), cols (), rows ());
00148 }
00149
00150 DiagMatrix
00151 real (const ComplexDiagMatrix& a)
00152 {
00153 DiagMatrix retval;
00154 int a_len = a.length ();
00155 if (a_len > 0)
00156 retval = DiagMatrix (mx_inline_real_dup (a.data (), a_len), a.rows (),
00157 a.cols ());
00158 return retval;
00159 }
00160
00161 DiagMatrix
00162 imag (const ComplexDiagMatrix& a)
00163 {
00164 DiagMatrix retval;
00165 int a_len = a.length ();
00166 if (a_len > 0)
00167 retval = DiagMatrix (mx_inline_imag_dup (a.data (), a_len), a.rows (),
00168 a.cols ());
00169 return retval;
00170 }
00171
00172 Matrix
00173 DiagMatrix::extract (int r1, int c1, int r2, int c2) const
00174 {
00175 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00176 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
00177
00178 int new_r = r2 - r1 + 1;
00179 int new_c = c2 - c1 + 1;
00180
00181 Matrix result (new_r, new_c);
00182
00183 for (int j = 0; j < new_c; j++)
00184 for (int i = 0; i < new_r; i++)
00185 result.elem (i, j) = elem (r1+i, c1+j);
00186
00187 return result;
00188 }
00189
00190
00191
00192 RowVector
00193 DiagMatrix::row (int i) const
00194 {
00195 int r = rows ();
00196 int c = cols ();
00197 if (i < 0 || i >= r)
00198 {
00199 (*current_liboctave_error_handler) ("invalid row selection");
00200 return RowVector ();
00201 }
00202
00203 RowVector retval (c, 0.0);
00204 if (r <= c || (r > c && i < c))
00205 retval.elem (i) = elem (i, i);
00206
00207 return retval;
00208 }
00209
00210 RowVector
00211 DiagMatrix::row (char *s) const
00212 {
00213 if (! s)
00214 {
00215 (*current_liboctave_error_handler) ("invalid row selection");
00216 return RowVector ();
00217 }
00218
00219 char c = *s;
00220 if (c == 'f' || c == 'F')
00221 return row (0);
00222 else if (c == 'l' || c == 'L')
00223 return row (rows () - 1);
00224 else
00225 {
00226 (*current_liboctave_error_handler) ("invalid row selection");
00227 return RowVector ();
00228 }
00229 }
00230
00231 ColumnVector
00232 DiagMatrix::column (int i) const
00233 {
00234 int r = rows ();
00235 int c = cols ();
00236 if (i < 0 || i >= c)
00237 {
00238 (*current_liboctave_error_handler) ("invalid column selection");
00239 return ColumnVector ();
00240 }
00241
00242 ColumnVector retval (r, 0.0);
00243 if (r >= c || (r < c && i < r))
00244 retval.elem (i) = elem (i, i);
00245
00246 return retval;
00247 }
00248
00249 ColumnVector
00250 DiagMatrix::column (char *s) const
00251 {
00252 if (! s)
00253 {
00254 (*current_liboctave_error_handler) ("invalid column selection");
00255 return ColumnVector ();
00256 }
00257
00258 char c = *s;
00259 if (c == 'f' || c == 'F')
00260 return column (0);
00261 else if (c == 'l' || c == 'L')
00262 return column (cols () - 1);
00263 else
00264 {
00265 (*current_liboctave_error_handler) ("invalid column selection");
00266 return ColumnVector ();
00267 }
00268 }
00269
00270 DiagMatrix
00271 DiagMatrix::inverse (void) const
00272 {
00273 int info;
00274 return inverse (info);
00275 }
00276
00277 DiagMatrix
00278 DiagMatrix::inverse (int &info) const
00279 {
00280 int r = rows ();
00281 int c = cols ();
00282 int len = length ();
00283 if (r != c)
00284 {
00285 (*current_liboctave_error_handler) ("inverse requires square matrix");
00286 return DiagMatrix ();
00287 }
00288
00289 DiagMatrix retval (r, c);
00290
00291 info = 0;
00292 for (int i = 0; i < len; i++)
00293 {
00294 if (elem (i, i) == 0.0)
00295 {
00296 info = -1;
00297 return *this;
00298 }
00299 else
00300 retval.elem (i, i) = 1.0 / elem (i, i);
00301 }
00302
00303 return retval;
00304 }
00305
00306
00307
00308
00309
00310 DiagMatrix
00311 operator * (const DiagMatrix& a, const DiagMatrix& b)
00312 {
00313 int a_nr = a.rows ();
00314 int a_nc = a.cols ();
00315
00316 int b_nr = b.rows ();
00317 int b_nc = b.cols ();
00318
00319 if (a_nc != b_nr)
00320 {
00321 gripe_nonconformant ("operaotr *", a_nr, a_nc, b_nr, b_nc);
00322 return DiagMatrix ();
00323 }
00324
00325 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
00326 return DiagMatrix (a_nr, a_nc, 0.0);
00327
00328 DiagMatrix c (a_nr, b_nc);
00329
00330 int len = a_nr < b_nc ? a_nr : b_nc;
00331
00332 for (int i = 0; i < len; i++)
00333 {
00334 double a_element = a.elem (i, i);
00335 double b_element = b.elem (i, i);
00336
00337 if (a_element == 0.0 || b_element == 0.0)
00338 c.elem (i, i) = 0.0;
00339 else if (a_element == 1.0)
00340 c.elem (i, i) = b_element;
00341 else if (b_element == 1.0)
00342 c.elem (i, i) = a_element;
00343 else
00344 c.elem (i, i) = a_element * b_element;
00345 }
00346
00347 return c;
00348 }
00349
00350
00351
00352 ColumnVector
00353 DiagMatrix::diag (void) const
00354 {
00355 return diag (0);
00356 }
00357
00358
00359
00360 ColumnVector
00361 DiagMatrix::diag (int k) const
00362 {
00363 int nnr = rows ();
00364 int nnc = cols ();
00365 if (k > 0)
00366 nnc -= k;
00367 else if (k < 0)
00368 nnr += k;
00369
00370 ColumnVector d;
00371
00372 if (nnr > 0 && nnc > 0)
00373 {
00374 int ndiag = (nnr < nnc) ? nnr : nnc;
00375
00376 d.resize (ndiag);
00377
00378 if (k > 0)
00379 {
00380 for (int i = 0; i < ndiag; i++)
00381 d.elem (i) = elem (i, i+k);
00382 }
00383 else if ( k < 0)
00384 {
00385 for (int i = 0; i < ndiag; i++)
00386 d.elem (i) = elem (i-k, i);
00387 }
00388 else
00389 {
00390 for (int i = 0; i < ndiag; i++)
00391 d.elem (i) = elem (i, i);
00392 }
00393 }
00394 else
00395 (*current_liboctave_error_handler)
00396 ("diag: requested diagonal out of range");
00397
00398 return d;
00399 }
00400
00401 std::ostream&
00402 operator << (std::ostream& os, const DiagMatrix& a)
00403 {
00404
00405
00406 for (int i = 0; i < a.rows (); i++)
00407 {
00408 for (int j = 0; j < a.cols (); j++)
00409 {
00410 if (i == j)
00411 os << " " << a.elem (i, i);
00412 else
00413 os << " " << 0.0;
00414 }
00415 os << "\n";
00416 }
00417 return os;
00418 }
00419
00420
00421
00422
00423
00424