Mercurial > hg > octave-lyh
annotate liboctave/dDiagMatrix.cc @ 7924:4976f66d469b
miscellaneous cleanup
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 11 Jul 2008 17:59:28 -0400 |
| parents | 82be108cc558 |
| children | 8b1a2555c4e2 |
| 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 // Diagonal Matrix class. |
| 458 | 38 |
| 2385 | 39 bool |
| 458 | 40 DiagMatrix::operator == (const DiagMatrix& a) const |
| 41 { | |
| 42 if (rows () != a.rows () || cols () != a.cols ()) | |
| 43 return 0; | |
| 44 | |
| 3769 | 45 return mx_inline_equal (data (), a.data (), length ()); |
| 458 | 46 } |
| 47 | |
| 2385 | 48 bool |
| 458 | 49 DiagMatrix::operator != (const DiagMatrix& a) const |
| 50 { | |
| 51 return !(*this == a); | |
| 52 } | |
| 53 | |
| 54 DiagMatrix& | |
| 55 DiagMatrix::fill (double val) | |
| 56 { | |
| 5275 | 57 for (octave_idx_type i = 0; i < length (); i++) |
| 458 | 58 elem (i, i) = val; |
| 59 return *this; | |
| 60 } | |
| 61 | |
| 62 DiagMatrix& | |
| 5275 | 63 DiagMatrix::fill (double val, octave_idx_type beg, octave_idx_type end) |
| 458 | 64 { |
| 65 if (beg < 0 || end >= length () || end < beg) | |
| 66 { | |
| 67 (*current_liboctave_error_handler) ("range error for fill"); | |
| 68 return *this; | |
| 69 } | |
| 70 | |
| 5275 | 71 for (octave_idx_type i = beg; i <= end; i++) |
| 458 | 72 elem (i, i) = val; |
| 73 | |
| 74 return *this; | |
| 75 } | |
| 76 | |
| 77 DiagMatrix& | |
| 78 DiagMatrix::fill (const ColumnVector& a) | |
| 79 { | |
| 5275 | 80 octave_idx_type len = length (); |
| 458 | 81 if (a.length () != len) |
| 82 { | |
| 83 (*current_liboctave_error_handler) ("range error for fill"); | |
| 84 return *this; | |
| 85 } | |
| 86 | |
| 5275 | 87 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 88 elem (i, i) = a.elem (i); |
| 89 | |
| 90 return *this; | |
| 91 } | |
| 92 | |
| 93 DiagMatrix& | |
| 94 DiagMatrix::fill (const RowVector& a) | |
| 95 { | |
| 5275 | 96 octave_idx_type len = length (); |
| 458 | 97 if (a.length () != len) |
| 98 { | |
| 99 (*current_liboctave_error_handler) ("range error for fill"); | |
| 100 return *this; | |
| 101 } | |
| 102 | |
| 5275 | 103 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 104 elem (i, i) = a.elem (i); |
| 105 | |
| 106 return *this; | |
| 107 } | |
| 108 | |
| 109 DiagMatrix& | |
| 5275 | 110 DiagMatrix::fill (const ColumnVector& a, octave_idx_type beg) |
| 458 | 111 { |
| 5275 | 112 octave_idx_type a_len = a.length (); |
| 458 | 113 if (beg < 0 || beg + a_len >= length ()) |
| 114 { | |
| 115 (*current_liboctave_error_handler) ("range error for fill"); | |
| 116 return *this; | |
| 117 } | |
| 118 | |
| 5275 | 119 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 120 elem (i+beg, i+beg) = a.elem (i); |
| 121 | |
| 122 return *this; | |
| 123 } | |
| 124 | |
| 125 DiagMatrix& | |
| 5275 | 126 DiagMatrix::fill (const RowVector& a, octave_idx_type beg) |
| 458 | 127 { |
| 5275 | 128 octave_idx_type a_len = a.length (); |
| 458 | 129 if (beg < 0 || beg + a_len >= length ()) |
| 130 { | |
| 131 (*current_liboctave_error_handler) ("range error for fill"); | |
| 132 return *this; | |
| 133 } | |
| 134 | |
| 5275 | 135 for (octave_idx_type i = 0; i < a_len; i++) |
| 458 | 136 elem (i+beg, i+beg) = a.elem (i); |
| 137 | |
| 138 return *this; | |
| 139 } | |
| 140 | |
| 141 DiagMatrix | |
| 1205 | 142 real (const ComplexDiagMatrix& a) |
| 143 { | |
| 144 DiagMatrix retval; | |
| 5275 | 145 octave_idx_type a_len = a.length (); |
| 1205 | 146 if (a_len > 0) |
| 3769 | 147 retval = DiagMatrix (mx_inline_real_dup (a.data (), a_len), a.rows (), |
| 1205 | 148 a.cols ()); |
| 149 return retval; | |
| 150 } | |
| 151 | |
| 152 DiagMatrix | |
| 153 imag (const ComplexDiagMatrix& a) | |
| 154 { | |
| 155 DiagMatrix retval; | |
| 5275 | 156 octave_idx_type a_len = a.length (); |
| 1205 | 157 if (a_len > 0) |
| 3769 | 158 retval = DiagMatrix (mx_inline_imag_dup (a.data (), a_len), a.rows (), |
| 1205 | 159 a.cols ()); |
| 160 return retval; | |
| 161 } | |
| 162 | |
| 458 | 163 Matrix |
| 5275 | 164 DiagMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const |
| 458 | 165 { |
| 5275 | 166 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
| 167 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
| 458 | 168 |
| 5275 | 169 octave_idx_type new_r = r2 - r1 + 1; |
| 170 octave_idx_type new_c = c2 - c1 + 1; | |
| 458 | 171 |
| 172 Matrix result (new_r, new_c); | |
| 173 | |
| 5275 | 174 for (octave_idx_type j = 0; j < new_c; j++) |
| 175 for (octave_idx_type i = 0; i < new_r; i++) | |
| 458 | 176 result.elem (i, j) = elem (r1+i, c1+j); |
| 177 | |
| 178 return result; | |
| 179 } | |
| 180 | |
| 181 // extract row or column i. | |
| 182 | |
| 183 RowVector | |
| 5275 | 184 DiagMatrix::row (octave_idx_type i) const |
| 458 | 185 { |
| 5275 | 186 octave_idx_type r = rows (); |
| 187 octave_idx_type c = cols (); | |
| 3504 | 188 if (i < 0 || i >= r) |
| 458 | 189 { |
| 190 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 191 return RowVector (); | |
| 192 } | |
| 193 | |
| 3504 | 194 RowVector retval (c, 0.0); |
| 195 if (r <= c || (r > c && i < c)) | |
| 458 | 196 retval.elem (i) = elem (i, i); |
| 197 | |
| 198 return retval; | |
| 199 } | |
| 200 | |
| 201 RowVector | |
| 202 DiagMatrix::row (char *s) const | |
| 203 { | |
| 533 | 204 if (! s) |
| 458 | 205 { |
| 206 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 207 return RowVector (); | |
| 208 } | |
| 209 | |
| 210 char c = *s; | |
| 211 if (c == 'f' || c == 'F') | |
| 5275 | 212 return row (static_cast<octave_idx_type>(0)); |
| 458 | 213 else if (c == 'l' || c == 'L') |
| 214 return row (rows () - 1); | |
| 215 else | |
| 216 { | |
| 217 (*current_liboctave_error_handler) ("invalid row selection"); | |
| 218 return RowVector (); | |
| 219 } | |
| 220 } | |
| 221 | |
| 222 ColumnVector | |
| 5275 | 223 DiagMatrix::column (octave_idx_type i) const |
| 458 | 224 { |
| 5275 | 225 octave_idx_type r = rows (); |
| 226 octave_idx_type c = cols (); | |
| 3504 | 227 if (i < 0 || i >= c) |
| 458 | 228 { |
| 229 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 230 return ColumnVector (); | |
| 231 } | |
| 232 | |
| 3504 | 233 ColumnVector retval (r, 0.0); |
| 234 if (r >= c || (r < c && i < r)) | |
| 458 | 235 retval.elem (i) = elem (i, i); |
| 236 | |
| 237 return retval; | |
| 238 } | |
| 239 | |
| 240 ColumnVector | |
| 241 DiagMatrix::column (char *s) const | |
| 242 { | |
| 533 | 243 if (! s) |
| 458 | 244 { |
| 245 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 246 return ColumnVector (); | |
| 247 } | |
| 248 | |
| 249 char c = *s; | |
| 250 if (c == 'f' || c == 'F') | |
| 5275 | 251 return column (static_cast<octave_idx_type>(0)); |
| 458 | 252 else if (c == 'l' || c == 'L') |
| 253 return column (cols () - 1); | |
| 254 else | |
| 255 { | |
| 256 (*current_liboctave_error_handler) ("invalid column selection"); | |
| 257 return ColumnVector (); | |
| 258 } | |
| 259 } | |
| 260 | |
| 261 DiagMatrix | |
| 262 DiagMatrix::inverse (void) const | |
| 263 { | |
| 264 int info; | |
| 265 return inverse (info); | |
| 266 } | |
| 267 | |
| 268 DiagMatrix | |
| 269 DiagMatrix::inverse (int &info) const | |
| 270 { | |
| 5275 | 271 octave_idx_type r = rows (); |
| 272 octave_idx_type c = cols (); | |
| 273 octave_idx_type len = length (); | |
| 3504 | 274 if (r != c) |
| 458 | 275 { |
| 276 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
| 277 return DiagMatrix (); | |
| 278 } | |
| 279 | |
| 3504 | 280 DiagMatrix retval (r, c); |
| 1627 | 281 |
| 458 | 282 info = 0; |
| 5275 | 283 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 284 { |
| 285 if (elem (i, i) == 0.0) | |
| 286 { | |
| 287 info = -1; | |
| 1627 | 288 return *this; |
| 458 | 289 } |
| 290 else | |
| 1627 | 291 retval.elem (i, i) = 1.0 / elem (i, i); |
| 458 | 292 } |
| 293 | |
| 1627 | 294 return retval; |
| 458 | 295 } |
| 296 | |
| 297 // diagonal matrix by diagonal matrix -> diagonal matrix operations | |
| 298 | |
| 299 // diagonal matrix by diagonal matrix -> diagonal matrix operations | |
| 300 | |
| 301 DiagMatrix | |
| 302 operator * (const DiagMatrix& a, const DiagMatrix& b) | |
| 303 { | |
| 5275 | 304 octave_idx_type a_nr = a.rows (); |
| 305 octave_idx_type a_nc = a.cols (); | |
| 2385 | 306 |
| 5275 | 307 octave_idx_type b_nr = b.rows (); |
| 308 octave_idx_type b_nc = b.cols (); | |
| 2385 | 309 |
| 3504 | 310 if (a_nc != b_nr) |
| 458 | 311 { |
| 3504 | 312 gripe_nonconformant ("operaotr *", a_nr, a_nc, b_nr, b_nc); |
| 458 | 313 return DiagMatrix (); |
| 314 } | |
| 315 | |
| 3504 | 316 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
| 317 return DiagMatrix (a_nr, a_nc, 0.0); | |
| 458 | 318 |
| 3504 | 319 DiagMatrix c (a_nr, b_nc); |
| 458 | 320 |
| 5275 | 321 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
| 458 | 322 |
| 5275 | 323 for (octave_idx_type i = 0; i < len; i++) |
| 458 | 324 { |
| 325 double a_element = a.elem (i, i); | |
| 326 double b_element = b.elem (i, i); | |
| 327 | |
| 328 if (a_element == 0.0 || b_element == 0.0) | |
| 329 c.elem (i, i) = 0.0; | |
| 330 else if (a_element == 1.0) | |
| 331 c.elem (i, i) = b_element; | |
| 332 else if (b_element == 1.0) | |
| 333 c.elem (i, i) = a_element; | |
| 334 else | |
| 335 c.elem (i, i) = a_element * b_element; | |
| 336 } | |
| 337 | |
| 338 return c; | |
| 339 } | |
| 340 | |
| 341 // other operations | |
| 342 | |
| 343 ColumnVector | |
| 5275 | 344 DiagMatrix::diag (octave_idx_type k) const |
| 458 | 345 { |
| 7924 | 346 ColumnVector d; |
| 347 | |
| 5275 | 348 octave_idx_type nnr = rows (); |
| 349 octave_idx_type nnc = cols (); | |
|
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
350 |
| 7924 | 351 if (nnr == 0 || nnc == 0) |
| 352 return d; | |
|
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
353 |
| 458 | 354 if (k > 0) |
| 355 nnc -= k; | |
| 356 else if (k < 0) | |
| 357 nnr += k; | |
| 358 | |
| 359 if (nnr > 0 && nnc > 0) | |
| 360 { | |
| 5275 | 361 octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; |
| 458 | 362 |
| 363 d.resize (ndiag); | |
| 364 | |
| 365 if (k > 0) | |
| 366 { | |
| 5275 | 367 for (octave_idx_type i = 0; i < ndiag; i++) |
| 458 | 368 d.elem (i) = elem (i, i+k); |
| 369 } | |
| 370 else if ( k < 0) | |
| 371 { | |
| 5275 | 372 for (octave_idx_type i = 0; i < ndiag; i++) |
| 458 | 373 d.elem (i) = elem (i-k, i); |
| 374 } | |
| 375 else | |
| 376 { | |
| 5275 | 377 for (octave_idx_type i = 0; i < ndiag; i++) |
| 458 | 378 d.elem (i) = elem (i, i); |
| 379 } | |
| 380 } | |
| 381 else | |
| 4513 | 382 (*current_liboctave_error_handler) |
| 383 ("diag: requested diagonal out of range"); | |
| 458 | 384 |
| 385 return d; | |
| 386 } | |
| 387 | |
| 3504 | 388 std::ostream& |
| 389 operator << (std::ostream& os, const DiagMatrix& a) | |
| 458 | 390 { |
| 391 // int field_width = os.precision () + 7; | |
| 1360 | 392 |
| 5275 | 393 for (octave_idx_type i = 0; i < a.rows (); i++) |
| 458 | 394 { |
| 5275 | 395 for (octave_idx_type j = 0; j < a.cols (); j++) |
| 458 | 396 { |
| 397 if (i == j) | |
| 398 os << " " /* setw (field_width) */ << a.elem (i, i); | |
| 399 else | |
| 400 os << " " /* setw (field_width) */ << 0.0; | |
| 401 } | |
| 402 os << "\n"; | |
| 403 } | |
| 404 return os; | |
| 405 } | |
| 406 | |
| 407 /* | |
| 408 ;;; Local Variables: *** | |
| 409 ;;; mode: C++ *** | |
| 410 ;;; End: *** | |
| 411 */ |