Mercurial > hg > octave-max
annotate liboctave/CDiagMatrix.cc @ 8375:e3c9102431a9
fix design problems of diag & perm matrix classes
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 05 Dec 2008 10:20:18 +0100 |
| parents | c3f7e2549abb |
| children | 20dfb885f877 |
| rev | line source |
|---|---|
| 1993 | 1 // DiagMatrix manipulations. |
| 458 | 2 /* |
| 3 | |
| 7017 | 4 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2004, |
| 5 2005, 2007 John W. Eaton | |
| 458 | 6 |
| 7 This file is part of Octave. | |
| 8 | |
| 9 Octave is free software; you can redistribute it and/or modify it | |
| 10 under the terms of the GNU General Public License as published by the | |
| 7016 | 11 Free Software Foundation; either version 3 of the License, or (at your |
| 12 option) any later version. | |
| 458 | 13 |
| 14 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 17 for more details. | |
| 18 | |
| 19 You should have received a copy of the GNU General Public License | |
| 7016 | 20 along with Octave; see the file COPYING. If not, see |
| 21 <http://www.gnu.org/licenses/>. | |
| 458 | 22 |
| 23 */ | |
| 24 | |
| 25 #ifdef HAVE_CONFIG_H | |
| 1192 | 26 #include <config.h> |
| 458 | 27 #endif |
| 28 | |
| 3503 | 29 #include <iostream> |
| 458 | 30 |
| 4669 | 31 #include "Array-util.h" |
| 1368 | 32 #include "lo-error.h" |
| 458 | 33 #include "mx-base.h" |
| 34 #include "mx-inlines.cc" | |
| 1650 | 35 #include "oct-cmplx.h" |
| 458 | 36 |
| 1360 | 37 // Complex Diagonal Matrix class |
| 458 | 38 |
| 39 ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a) | |
| 1989 | 40 : MDiagArray2<Complex> (a.rows (), a.cols ()) |
| 458 | 41 { |
| 5275 | 42 for (octave_idx_type i = 0; i < length (); i++) |
| 458 | 43 elem (i, i) = a.elem (i, i); |
| 44 } | |
| 45 | |
| 2386 | 46 bool |
| 458 | 47 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const |
| 48 { | |
| 49 if (rows () != a.rows () || cols () != a.cols ()) | |
| 50 return 0; | |
| 51 | |
| 3769 | 52 return mx_inline_equal (data (), a.data (), length ()); |
| 458 | 53 } |
| 54 | |
| 2386 | 55 bool |
| 458 | 56 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const |
| 57 { | |
| 58 return !(*this == a); | |
| 59 } | |
| 60 | |
| 61 ComplexDiagMatrix& | |
| 62 ComplexDiagMatrix::fill (double val) | |
| 63 { | |
| 5275 | 64 for (octave_idx_type i = 0; i < length (); i++) |
| 458 | 65 elem (i, i) = val; |
| 66 return *this; | |
| 67 } | |
| 68 | |
| 69 ComplexDiagMatrix& | |
| 70 ComplexDiagMatrix::fill (const Complex& val) | |
| 71 { | |
| 5275 | 72 for (octave_idx_type i = 0; i < length (); i++) |
| 458 | 73 elem (i, i) = val; |
| 74 return *this; | |
| 75 } | |
| 76 | |
| 77 ComplexDiagMatrix& | |
| 5275 | 78 ComplexDiagMatrix::fill (double val, octave_idx_type beg, octave_idx_type end) |
| 458 | 79 { |
| 80 if (beg < 0 || end >= length () || end < beg) | |
| 81 { | |
| 82 (*current_liboctave_error_handler) ("range error for fill"); | |
| 83 return *this; | |
| 84 } | |
| 85 | |
| 5275 | 86 for (octave_idx_type i = beg; i <= end; i++) |
| 458 | 87 elem (i, i) = val; |
| 88 | |
| 89 return *this; | |
| 90 } | |
| 91 | |
| 92 ComplexDiagMatrix& | |
| 5275 | 93 ComplexDiagMatrix::fill (const Complex& val, octave_idx_type beg, octave_idx_type end) |
| 458 | 94 { |
| 95 if (beg < 0 || end >= length () || end < beg) | |
| 96 { | |
| 97 (*current_liboctave_error_handler) ("range error for fill"); | |
| 98 return *this; | |
| 99 } | |
| 100 | |
| 5275 | 101 for (octave_idx_type i = beg; i <= end; i++) |
| 458 | 102 elem (i, i) = val; |
| 103 | |
| 104 return *this; | |
| 105 } | |
| 106 | |
| 107 ComplexDiagMatrix& | |
| 108 ComplexDiagMatrix::fill (const ColumnVector& a) | |
| 109 { | |
| 5275 | 110 octave_idx_type len = length (); |
| 458 | 111 if (a.length () != len) |
| 112 { | |
| 113 (*current_liboctave_error_handler) ("range error for fill"); | |
| 114 return *this; | |
| 115 } | |
| 116 | |
| 5275 | 117 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 118 elem (i, i) = a.elem (i); |
| 119 | |
| 120 return *this; | |
| 121 } | |
| 122 | |
| 123 ComplexDiagMatrix& | |
| 124 ComplexDiagMatrix::fill (const ComplexColumnVector& a) | |
| 125 { | |
| 5275 | 126 octave_idx_type len = length (); |
| 458 | 127 if (a.length () != len) |
| 128 { | |
| 129 (*current_liboctave_error_handler) ("range error for fill"); | |
| 130 return *this; | |
| 131 } | |
| 132 | |
| 5275 | 133 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 134 elem (i, i) = a.elem (i); |
| 135 | |
| 136 return *this; | |
| 137 } | |
| 138 | |
| 139 ComplexDiagMatrix& | |
| 140 ComplexDiagMatrix::fill (const RowVector& a) | |
| 141 { | |
| 5275 | 142 octave_idx_type len = length (); |
| 458 | 143 if (a.length () != len) |
| 144 { | |
| 145 (*current_liboctave_error_handler) ("range error for fill"); | |
| 146 return *this; | |
| 147 } | |
| 148 | |
| 5275 | 149 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 150 elem (i, i) = a.elem (i); |
| 151 | |
| 152 return *this; | |
| 153 } | |
| 154 | |
| 155 ComplexDiagMatrix& | |
| 156 ComplexDiagMatrix::fill (const ComplexRowVector& a) | |
| 157 { | |
| 5275 | 158 octave_idx_type len = length (); |
| 458 | 159 if (a.length () != len) |
| 160 { | |
| 161 (*current_liboctave_error_handler) ("range error for fill"); | |
| 162 return *this; | |
| 163 } | |
| 164 | |
| 5275 | 165 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 166 elem (i, i) = a.elem (i); |
| 167 | |
| 168 return *this; | |
| 169 } | |
| 170 | |
| 171 ComplexDiagMatrix& | |
| 5275 | 172 ComplexDiagMatrix::fill (const ColumnVector& a, octave_idx_type beg) |
| 458 | 173 { |
| 5275 | 174 octave_idx_type a_len = a.length (); |
| 458 | 175 if (beg < 0 || beg + a_len >= length ()) |
| 176 { | |
| 177 (*current_liboctave_error_handler) ("range error for fill"); | |
| 178 return *this; | |
| 179 } | |
| 180 | |
| 5275 | 181 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 182 elem (i+beg, i+beg) = a.elem (i); |
| 183 | |
| 184 return *this; | |
| 185 } | |
| 186 | |
| 187 ComplexDiagMatrix& | |
| 5275 | 188 ComplexDiagMatrix::fill (const ComplexColumnVector& a, octave_idx_type beg) |
| 458 | 189 { |
| 5275 | 190 octave_idx_type a_len = a.length (); |
| 458 | 191 if (beg < 0 || beg + a_len >= length ()) |
| 192 { | |
| 193 (*current_liboctave_error_handler) ("range error for fill"); | |
| 194 return *this; | |
| 195 } | |
| 196 | |
| 5275 | 197 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 198 elem (i+beg, i+beg) = a.elem (i); |
| 199 | |
| 200 return *this; | |
| 201 } | |
| 202 | |
| 203 ComplexDiagMatrix& | |
| 5275 | 204 ComplexDiagMatrix::fill (const RowVector& a, octave_idx_type beg) |
| 458 | 205 { |
| 5275 | 206 octave_idx_type a_len = a.length (); |
| 458 | 207 if (beg < 0 || beg + a_len >= length ()) |
| 208 { | |
| 209 (*current_liboctave_error_handler) ("range error for fill"); | |
| 210 return *this; | |
| 211 } | |
| 212 | |
| 5275 | 213 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 214 elem (i+beg, i+beg) = a.elem (i); |
| 215 | |
| 216 return *this; | |
| 217 } | |
| 218 | |
| 219 ComplexDiagMatrix& | |
| 5275 | 220 ComplexDiagMatrix::fill (const ComplexRowVector& a, octave_idx_type beg) |
| 458 | 221 { |
| 5275 | 222 octave_idx_type a_len = a.length (); |
| 458 | 223 if (beg < 0 || beg + a_len >= length ()) |
| 224 { | |
| 225 (*current_liboctave_error_handler) ("range error for fill"); | |
| 226 return *this; | |
| 227 } | |
| 228 | |
| 5275 | 229 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 230 elem (i+beg, i+beg) = a.elem (i); |
| 231 | |
| 232 return *this; | |
| 233 } | |
| 234 | |
|
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
235 DiagMatrix |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
236 ComplexDiagMatrix::abs (void) const |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
237 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
238 DiagMatrix retval (rows (), cols ()); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
239 for (octave_idx_type i = 0; i < rows (); i++) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
240 retval(i, i) = std::abs (elem (i, i)); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
241 return retval; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
242 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
243 |
| 458 | 244 ComplexDiagMatrix |
| 245 conj (const ComplexDiagMatrix& a) | |
| 246 { | |
| 247 ComplexDiagMatrix retval; | |
| 5275 | 248 octave_idx_type a_len = a.length (); |
| 458 | 249 if (a_len > 0) |
| 3769 | 250 retval = ComplexDiagMatrix (mx_inline_conj_dup (a.data (), a_len), |
| 458 | 251 a.rows (), a.cols ()); |
| 252 return retval; | |
| 253 } | |
| 254 | |
| 255 // resize is the destructive analog for this one | |
| 256 | |
| 257 ComplexMatrix | |
| 5275 | 258 ComplexDiagMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const |
| 458 | 259 { |
| 5275 | 260 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
| 261 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
| 458 | 262 |
| 5275 | 263 octave_idx_type new_r = r2 - r1 + 1; |
| 264 octave_idx_type new_c = c2 - c1 + 1; | |
| 458 | 265 |
| 266 ComplexMatrix result (new_r, new_c); | |
| 267 | |
| 5275 | 268 for (octave_idx_type j = 0; j < new_c; j++) |
| 269 for (octave_idx_type i = 0; i < new_r; i++) | |
| 458 | 270 result.elem (i, j) = elem (r1+i, c1+j); |
| 271 | |
| 272 return result; | |
| 273 } | |
| 274 | |
| 275 // extract row or column i. | |
| 276 | |
| 277 ComplexRowVector | |
| 5275 | 278 ComplexDiagMatrix::row (octave_idx_type i) const |
| 458 | 279 { |
| 5275 | 280 octave_idx_type r = rows (); |
| 281 octave_idx_type c = cols (); | |
| 3504 | 282 if (i < 0 || i >= r) |
| 458 | 283 { |
| 284 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 3585 | 285 return ComplexRowVector (); |
| 458 | 286 } |
| 287 | |
| 3504 | 288 ComplexRowVector retval (c, 0.0); |
| 289 if (r <= c || (r > c && i < c)) | |
| 458 | 290 retval.elem (i) = elem (i, i); |
| 291 | |
| 292 return retval; | |
| 293 } | |
| 294 | |
| 295 ComplexRowVector | |
| 296 ComplexDiagMatrix::row (char *s) const | |
| 297 { | |
| 533 | 298 if (! s) |
| 458 | 299 { |
| 300 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 301 return ComplexRowVector (); | |
| 302 } | |
| 303 | |
| 304 char c = *s; | |
| 305 if (c == 'f' || c == 'F') | |
| 5275 | 306 return row (static_cast<octave_idx_type>(0)); |
| 458 | 307 else if (c == 'l' || c == 'L') |
| 308 return row (rows () - 1); | |
| 309 else | |
| 310 { | |
| 311 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 312 return ComplexRowVector (); | |
| 313 } | |
| 314 } | |
| 315 | |
| 316 ComplexColumnVector | |
| 5275 | 317 ComplexDiagMatrix::column (octave_idx_type i) const |
| 458 | 318 { |
| 5275 | 319 octave_idx_type r = rows (); |
| 320 octave_idx_type c = cols (); | |
| 3504 | 321 if (i < 0 || i >= c) |
| 458 | 322 { |
| 323 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 3585 | 324 return ComplexColumnVector (); |
| 458 | 325 } |
| 326 | |
| 3504 | 327 ComplexColumnVector retval (r, 0.0); |
| 328 if (r >= c || (r < c && i < r)) | |
| 458 | 329 retval.elem (i) = elem (i, i); |
| 330 | |
| 331 return retval; | |
| 332 } | |
| 333 | |
| 334 ComplexColumnVector | |
| 335 ComplexDiagMatrix::column (char *s) const | |
| 336 { | |
| 533 | 337 if (! s) |
| 458 | 338 { |
| 339 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 3585 | 340 return ComplexColumnVector (); |
| 458 | 341 } |
| 342 | |
| 343 char c = *s; | |
| 344 if (c == 'f' || c == 'F') | |
| 5275 | 345 return column (static_cast<octave_idx_type>(0)); |
| 458 | 346 else if (c == 'l' || c == 'L') |
| 347 return column (cols () - 1); | |
| 348 else | |
| 349 { | |
| 350 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 3585 | 351 return ComplexColumnVector (); |
| 458 | 352 } |
| 353 } | |
| 354 | |
| 355 ComplexDiagMatrix | |
| 356 ComplexDiagMatrix::inverse (void) const | |
| 357 { | |
| 358 int info; | |
| 359 return inverse (info); | |
| 360 } | |
| 361 | |
| 362 ComplexDiagMatrix | |
| 363 ComplexDiagMatrix::inverse (int& info) const | |
| 364 { | |
| 5275 | 365 octave_idx_type r = rows (); |
| 366 octave_idx_type c = cols (); | |
| 3504 | 367 if (r != c) |
| 458 | 368 { |
| 369 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
| 3585 | 370 return ComplexDiagMatrix (); |
| 458 | 371 } |
| 372 | |
| 3504 | 373 ComplexDiagMatrix retval (r, c); |
| 458 | 374 |
| 375 info = 0; | |
| 5275 | 376 for (octave_idx_type i = 0; i < length (); i++) |
| 458 | 377 { |
| 378 if (elem (i, i) == 0.0) | |
| 379 { | |
| 380 info = -1; | |
| 381 return *this; | |
| 382 } | |
| 383 else | |
| 384 retval.elem (i, i) = 1.0 / elem (i, i); | |
| 385 } | |
| 386 | |
| 1627 | 387 return retval; |
| 458 | 388 } |
| 389 | |
|
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
390 bool |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
391 ComplexDiagMatrix::all_elements_are_real (void) const |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
392 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
393 octave_idx_type len = length (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
394 for (octave_idx_type i = 0; i < len; i++) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
395 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
396 double ip = std::imag (elem (i, i)); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
397 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
398 if (ip != 0.0 || lo_ieee_signbit (ip)) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
399 return false; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
400 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
401 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
402 return true; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
403 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
404 |
| 458 | 405 // diagonal matrix by diagonal matrix -> diagonal matrix operations |
| 406 | |
| 407 ComplexDiagMatrix& | |
| 408 ComplexDiagMatrix::operator += (const DiagMatrix& a) | |
| 409 { | |
| 5275 | 410 octave_idx_type r = rows (); |
| 411 octave_idx_type c = cols (); | |
| 2386 | 412 |
| 5275 | 413 octave_idx_type a_nr = a.rows (); |
| 414 octave_idx_type a_nc = a.cols (); | |
| 2386 | 415 |
| 3504 | 416 if (r != a_nr || c != a_nc) |
| 458 | 417 { |
| 3504 | 418 gripe_nonconformant ("operator +=", r, c, a_nr, a_nc); |
| 458 | 419 return *this; |
| 420 } | |
| 421 | |
| 3504 | 422 if (r == 0 || c == 0) |
| 458 | 423 return *this; |
| 424 | |
| 425 Complex *d = fortran_vec (); // Ensures only one reference to my privates! | |
| 426 | |
| 3769 | 427 mx_inline_add2 (d, a.data (), length ()); |
| 458 | 428 return *this; |
| 429 } | |
| 430 | |
| 431 ComplexDiagMatrix | |
| 432 operator * (const ComplexDiagMatrix& a, const DiagMatrix& b) | |
| 433 { | |
| 5275 | 434 octave_idx_type a_nr = a.rows (); |
| 435 octave_idx_type a_nc = a.cols (); | |
| 2386 | 436 |
| 5275 | 437 octave_idx_type b_nr = b.rows (); |
| 438 octave_idx_type b_nc = b.cols (); | |
| 2386 | 439 |
| 3504 | 440 if (a_nc != b_nr) |
| 458 | 441 { |
| 3504 | 442 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
| 458 | 443 return ComplexDiagMatrix (); |
| 444 } | |
| 445 | |
| 3504 | 446 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
| 447 return ComplexDiagMatrix (a_nr, a_nc, 0.0); | |
| 458 | 448 |
| 3504 | 449 ComplexDiagMatrix c (a_nr, b_nc); |
| 458 | 450 |
| 5275 | 451 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
| 458 | 452 |
| 5275 | 453 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 454 { |
| 455 Complex a_element = a.elem (i, i); | |
| 456 double b_element = b.elem (i, i); | |
| 457 | |
|
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
458 c.elem (i, i) = a_element * b_element; |
| 458 | 459 } |
| 460 | |
| 461 return c; | |
| 462 } | |
| 463 | |
| 464 ComplexDiagMatrix | |
| 1205 | 465 operator * (const DiagMatrix& a, const ComplexDiagMatrix& b) |
| 466 { | |
| 5275 | 467 octave_idx_type a_nr = a.rows (); |
| 468 octave_idx_type a_nc = a.cols (); | |
| 2386 | 469 |
| 5275 | 470 octave_idx_type b_nr = b.rows (); |
| 471 octave_idx_type b_nc = b.cols (); | |
| 2386 | 472 |
| 3504 | 473 if (a_nc != b_nr) |
| 1205 | 474 { |
| 3504 | 475 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
| 1205 | 476 return ComplexDiagMatrix (); |
| 477 } | |
| 478 | |
| 3504 | 479 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
| 480 return ComplexDiagMatrix (a_nr, a_nc, 0.0); | |
| 1205 | 481 |
| 3504 | 482 ComplexDiagMatrix c (a_nr, b_nc); |
| 1205 | 483 |
| 5275 | 484 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
| 1205 | 485 |
| 5275 | 486 for (octave_idx_type i = 0; i < len; i++) |
| 1205 | 487 { |
| 488 double a_element = a.elem (i, i); | |
| 489 Complex b_element = b.elem (i, i); | |
| 490 | |
|
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
491 c.elem (i, i) = a_element * b_element; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
492 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
493 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
494 return c; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
495 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
496 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
497 ComplexDiagMatrix |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
498 operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
499 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
500 octave_idx_type a_nr = a.rows (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
501 octave_idx_type a_nc = a.cols (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
502 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
503 octave_idx_type b_nr = b.rows (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
504 octave_idx_type b_nc = b.cols (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
505 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
506 if (a_nc != b_nr) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
507 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
508 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
509 return ComplexDiagMatrix (); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
510 } |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
511 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
512 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
513 return ComplexDiagMatrix (a_nr, a_nc, 0.0); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
514 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
515 ComplexDiagMatrix c (a_nr, b_nc); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
516 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
517 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
518 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
519 for (octave_idx_type i = 0; i < len; i++) |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
520 { |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
521 Complex a_element = a.elem (i, i); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
522 Complex b_element = b.elem (i, i); |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
523 |
|
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
524 c.elem (i, i) = a_element * b_element; |
| 1205 | 525 } |
| 526 | |
| 527 return c; | |
| 528 } | |
| 529 | |
| 458 | 530 // other operations |
| 531 | |
|
8371
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
532 ComplexDET |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
533 ComplexDiagMatrix::determinant (void) const |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
534 { |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
535 ComplexDET det (1.0); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
536 if (rows () != cols ()) |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
537 { |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
538 (*current_liboctave_error_handler) ("determinant requires square matrix"); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
539 det = ComplexDET (0.0); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
540 } |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
541 else |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
542 { |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
543 octave_idx_type len = length (); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
544 for (octave_idx_type i = 0; i < len; i++) |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
545 det *= elem (i, i); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
546 } |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
547 |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
548 return det; |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
549 } |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
550 |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
551 double |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
552 ComplexDiagMatrix::rcond (void) const |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
553 { |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
554 ColumnVector av = diag (0).map (std::abs); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
555 double amx = av.max (), amn = av.min (); |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
556 return amx == 0 ? 0.0 : amn / amx; |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
557 } |
|
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
558 |
| 458 | 559 // i/o |
| 560 | |
| 3504 | 561 std::ostream& |
| 562 operator << (std::ostream& os, const ComplexDiagMatrix& a) | |
| 458 | 563 { |
| 564 Complex ZERO (0.0); | |
| 565 // int field_width = os.precision () + 7; | |
| 5275 | 566 for (octave_idx_type i = 0; i < a.rows (); i++) |
| 458 | 567 { |
| 5275 | 568 for (octave_idx_type j = 0; j < a.cols (); j++) |
| 458 | 569 { |
| 570 if (i == j) | |
| 571 os << " " /* setw (field_width) */ << a.elem (i, i); | |
| 572 else | |
| 573 os << " " /* setw (field_width) */ << ZERO; | |
| 574 } | |
| 575 os << "\n"; | |
| 576 } | |
| 577 return os; | |
| 578 } | |
| 579 | |
| 580 /* | |
| 581 ;;; Local Variables: *** | |
| 582 ;;; mode: C++ *** | |
| 583 ;;; End: *** | |
| 584 */ |