|
3
|
1 // Matrix manipulations. -*- C++ -*- |
|
|
2 /* |
|
|
3 |
|
|
4 Copyright (C) 1992, 1993 John W. Eaton |
|
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 // I\'m not sure how this is supposed to work if the .h file declares |
|
|
25 // several classes, each of which is defined in a separate file... |
|
|
26 // |
|
|
27 // #ifdef __GNUG__ |
|
|
28 // #pragma implementation |
|
|
29 // #endif |
|
|
30 |
|
|
31 #include "Matrix.h" |
|
|
32 #include "mx-inlines.cc" |
|
227
|
33 #include "lo-error.h" |
|
232
|
34 #include "f77-uscore.h" |
|
|
35 |
|
|
36 // Fortran functions we call. |
|
|
37 |
|
|
38 extern "C" |
|
|
39 { |
|
|
40 int F77_FCN (dgemm) (const char*, const char*, const int*, |
|
|
41 const int*, const int*, const double*, |
|
|
42 const double*, const int*, const double*, |
|
|
43 const int*, const double*, double*, const int*, |
|
|
44 long, long); |
|
|
45 |
|
|
46 int F77_FCN (dgemv) (const char*, const int*, const int*, |
|
|
47 const double*, const double*, const int*, |
|
|
48 const double*, const int*, const double*, |
|
|
49 double*, const int*, long); |
|
|
50 |
|
|
51 int F77_FCN (dgeco) (double*, const int*, const int*, int*, double*, |
|
|
52 double*); |
|
|
53 |
|
|
54 int F77_FCN (dgesl) (const double*, const int*, const int*, |
|
|
55 const int*, double*, const int*); |
|
|
56 |
|
|
57 int F77_FCN (dgedi) (double*, const int*, const int*, const int*, |
|
|
58 double*, double*, const int*); |
|
|
59 |
|
|
60 int F77_FCN (dgelss) (const int*, const int*, const int*, double*, |
|
|
61 const int*, double*, const int*, double*, |
|
|
62 const double*, int*, double*, const int*, |
|
|
63 int*); |
|
|
64 |
|
|
65 /* |
|
|
66 * f2c translates complex*16 as |
|
|
67 * |
|
|
68 * typedef struct { doublereal re, im; } doublecomplex; |
|
|
69 * |
|
|
70 * and Complex.h from libg++ uses |
|
|
71 * |
|
|
72 * protected: |
|
|
73 * double re; |
|
|
74 * double im; |
|
|
75 * |
|
|
76 * as the only data members, so this should work (fingers crossed that |
|
|
77 * things don't change). |
|
|
78 */ |
|
|
79 |
|
|
80 int F77_FCN (zgemm) (const char*, const char*, const int*, |
|
|
81 const int*, const int*, const Complex*, |
|
|
82 const Complex*, const int*, const Complex*, |
|
|
83 const int*, const Complex*, Complex*, const int*, |
|
|
84 long, long); |
|
|
85 |
|
|
86 int F77_FCN (zgemv) (const char*, const int*, const int*, |
|
|
87 const Complex*, const Complex*, const int*, |
|
|
88 const Complex*, const int*, const Complex*, |
|
|
89 Complex*, const int*, long); |
|
|
90 |
|
|
91 int F77_FCN (zgeco) (Complex*, const int*, const int*, int*, |
|
|
92 double*, Complex*); |
|
|
93 |
|
|
94 int F77_FCN (zgedi) (Complex*, const int*, const int*, int*, |
|
|
95 Complex*, Complex*, const int*); |
|
|
96 |
|
|
97 int F77_FCN (zgesl) (Complex*, const int*, const int*, int*, |
|
|
98 Complex*, const int*); |
|
|
99 |
|
|
100 int F77_FCN (zgelss) (const int*, const int*, const int*, Complex*, |
|
|
101 const int*, Complex*, const int*, double*, |
|
|
102 const double*, int*, Complex*, const int*, |
|
|
103 double*, int*); |
|
|
104 |
|
|
105 // Note that the original complex fft routines were not written for |
|
|
106 // double complex arguments. They have been modified by adding an |
|
|
107 // implicit double precision (a-h,o-z) statement at the beginning of |
|
|
108 // each subroutine. |
|
|
109 |
|
|
110 int F77_FCN (cffti) (const int*, Complex*); |
|
|
111 |
|
|
112 int F77_FCN (cfftf) (const int*, Complex*, Complex*); |
|
|
113 |
|
|
114 int F77_FCN (cfftb) (const int*, Complex*, Complex*); |
|
|
115 } |
|
3
|
116 |
|
|
117 /* |
|
|
118 * Matrix class. |
|
|
119 */ |
|
|
120 |
|
|
121 Matrix::Matrix (int r, int c) |
|
|
122 { |
|
|
123 if (r < 0 || c < 0) |
|
227
|
124 { |
|
|
125 (*current_liboctave_error_handler) |
|
|
126 ("can't construct matrix with negative dimensions"); |
|
|
127 nr = 0; |
|
|
128 nc = 0; |
|
|
129 len = 0; |
|
|
130 data = (double *) NULL; |
|
|
131 return; |
|
|
132 } |
|
3
|
133 |
|
|
134 nr = r; |
|
|
135 nc = c; |
|
|
136 len = nr * nc; |
|
|
137 if (len > 0) |
|
|
138 data = new double [len]; |
|
|
139 else |
|
|
140 data = (double *) NULL; |
|
|
141 } |
|
|
142 |
|
|
143 Matrix::Matrix (int r, int c, double val) |
|
|
144 { |
|
|
145 if (r < 0 || c < 0) |
|
227
|
146 { |
|
|
147 (*current_liboctave_error_handler) |
|
|
148 ("can't construct matrix with negative dimensions"); |
|
|
149 nr = 0; |
|
|
150 nc = 0; |
|
|
151 len = 0; |
|
|
152 data = (double *) NULL; |
|
|
153 return; |
|
|
154 } |
|
3
|
155 |
|
|
156 nr = r; |
|
|
157 nc = c; |
|
|
158 len = nr * nc; |
|
|
159 if (len > 0) |
|
|
160 { |
|
|
161 data = new double [len]; |
|
|
162 copy (data, len, val); |
|
|
163 } |
|
|
164 else |
|
|
165 data = (double *) NULL; |
|
|
166 } |
|
|
167 |
|
|
168 Matrix::Matrix (const Matrix& a) |
|
|
169 { |
|
|
170 nr = a.nr; |
|
|
171 nc = a.nc; |
|
|
172 len = a.len; |
|
|
173 if (len > 0) |
|
|
174 { |
|
|
175 data = new double [len]; |
|
|
176 copy (data, a.data, len); |
|
|
177 } |
|
|
178 else |
|
|
179 data = (double *) NULL; |
|
|
180 } |
|
|
181 |
|
|
182 Matrix::Matrix (const DiagMatrix& a) |
|
|
183 { |
|
|
184 nr = a.nr; |
|
|
185 nc = a.nc; |
|
|
186 len = nr * nc; |
|
|
187 if (len > 0) |
|
|
188 { |
|
|
189 data = new double [len]; |
|
|
190 copy (data, len, 0.0); |
|
|
191 for (int i = 0; i < a.len; i++) |
|
|
192 data[nr*i+i] = a.data[i]; |
|
|
193 } |
|
|
194 else |
|
|
195 data = (double *) NULL; |
|
|
196 } |
|
|
197 |
|
|
198 Matrix::Matrix (double a) |
|
|
199 { |
|
|
200 nr = 1; |
|
|
201 nc = 1; |
|
|
202 len = 1; |
|
|
203 data = new double [1]; |
|
|
204 data[0] = a; |
|
|
205 } |
|
|
206 |
|
|
207 Matrix& |
|
|
208 Matrix::operator = (const Matrix& a) |
|
|
209 { |
|
|
210 if (this != &a) |
|
|
211 { |
|
|
212 delete [] data; |
|
|
213 nr = a.nr; |
|
|
214 nc = a.nc; |
|
|
215 len = a.len; |
|
|
216 if (len > 0) |
|
|
217 { |
|
|
218 data = new double [len]; |
|
|
219 copy (data, a.data, len); |
|
|
220 } |
|
|
221 else |
|
|
222 data = (double *) NULL; |
|
|
223 } |
|
|
224 return *this; |
|
|
225 } |
|
|
226 |
|
227
|
227 double& |
|
|
228 Matrix::checkelem (int r, int c) |
|
|
229 { |
|
|
230 #ifndef NO_RANGE_CHECK |
|
|
231 if (r < 0 || r >= nr || c < 0 || c >= nc) |
|
|
232 { |
|
|
233 (*current_liboctave_error_handler) ("range error"); |
|
|
234 static double foo = 0.0; |
|
|
235 return foo; |
|
|
236 } |
|
|
237 #endif |
|
|
238 |
|
|
239 return elem (r, c); |
|
|
240 } |
|
|
241 |
|
|
242 double |
|
|
243 Matrix::checkelem (int r, int c) const |
|
|
244 { |
|
|
245 #ifndef NO_RANGE_CHECK |
|
|
246 if (r < 0 || r >= nr || c < 0 || c >= nc) |
|
|
247 { |
|
|
248 (*current_liboctave_error_handler) ("range error"); |
|
|
249 return 0.0; |
|
|
250 } |
|
|
251 #endif |
|
|
252 |
|
|
253 return elem (r, c); |
|
|
254 } |
|
|
255 |
|
3
|
256 Matrix& |
|
|
257 Matrix::resize (int r, int c) |
|
|
258 { |
|
|
259 if (r < 0 || c < 0) |
|
227
|
260 { |
|
|
261 (*current_liboctave_error_handler) |
|
|
262 ("can't resize to negative dimensions"); |
|
|
263 return *this; |
|
|
264 } |
|
3
|
265 |
|
|
266 int new_len = r * c; |
|
|
267 double* new_data = (double *) NULL; |
|
|
268 if (new_len > 0) |
|
|
269 { |
|
|
270 new_data = new double [new_len]; |
|
|
271 |
|
|
272 int min_r = nr < r ? nr : r; |
|
|
273 int min_c = nc < c ? nc : c; |
|
|
274 |
|
|
275 for (int j = 0; j < min_c; j++) |
|
|
276 for (int i = 0; i < min_r; i++) |
|
|
277 new_data[r*j+i] = elem (i, j); |
|
|
278 } |
|
|
279 |
|
|
280 delete [] data; |
|
|
281 nr = r; |
|
|
282 nc = c; |
|
|
283 len = new_len; |
|
|
284 data = new_data; |
|
|
285 |
|
|
286 return *this; |
|
|
287 } |
|
|
288 |
|
|
289 Matrix& |
|
|
290 Matrix::resize (int r, int c, double val) |
|
|
291 { |
|
|
292 if (r < 0 || c < 0) |
|
227
|
293 { |
|
|
294 (*current_liboctave_error_handler) |
|
|
295 ("can't resize to negative dimensions"); |
|
|
296 return *this; |
|
|
297 } |
|
3
|
298 |
|
|
299 int new_len = r * c; |
|
|
300 double *new_data = (double *) NULL; |
|
|
301 if (new_len > 0) |
|
|
302 { |
|
|
303 new_data = new double [new_len]; |
|
|
304 |
|
|
305 // There may be faster or cleaner ways to do this. |
|
|
306 |
|
|
307 if (r > nr || c > nc) |
|
|
308 copy (new_data, new_len, val); |
|
|
309 |
|
|
310 int min_r = nr < r ? nr : r; |
|
|
311 int min_c = nc < c ? nc : c; |
|
|
312 |
|
|
313 for (int j = 0; j < min_c; j++) |
|
|
314 for (int i = 0; i < min_r; i++) |
|
|
315 new_data[r*j+i] = elem (i, j); |
|
|
316 } |
|
|
317 |
|
|
318 delete [] data; |
|
|
319 nr = r; |
|
|
320 nc = c; |
|
|
321 len = new_len; |
|
|
322 data = new_data; |
|
|
323 |
|
|
324 return *this; |
|
|
325 } |
|
|
326 |
|
|
327 int |
|
|
328 Matrix::operator == (const Matrix& a) const |
|
|
329 { |
|
|
330 if (nr != a.nr || nc != a.nc) |
|
|
331 return 0; |
|
|
332 |
|
|
333 return equal (data, a.data, len); |
|
|
334 } |
|
|
335 |
|
|
336 int |
|
|
337 Matrix::operator != (const Matrix& a) const |
|
|
338 { |
|
|
339 return !(*this == a); |
|
|
340 } |
|
|
341 |
|
|
342 Matrix& |
|
|
343 Matrix::insert (const Matrix& a, int r, int c) |
|
|
344 { |
|
|
345 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
346 { |
|
|
347 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
348 return *this; |
|
|
349 } |
|
3
|
350 |
|
|
351 for (int j = 0; j < a.nc; j++) |
|
|
352 for (int i = 0; i < a.nr; i++) |
|
|
353 elem (r+i, c+j) = a.elem (i, j); |
|
|
354 |
|
|
355 return *this; |
|
|
356 } |
|
|
357 |
|
|
358 Matrix& |
|
|
359 Matrix::insert (const RowVector& a, int r, int c) |
|
|
360 { |
|
|
361 if (r < 0 || r >= nr || c < 0 || c + a.len - 1 > nc) |
|
227
|
362 { |
|
|
363 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
364 return *this; |
|
|
365 } |
|
3
|
366 |
|
|
367 for (int i = 0; i < a.len; i++) |
|
|
368 elem (r, c+i) = a.data[i]; |
|
|
369 |
|
|
370 return *this; |
|
|
371 } |
|
|
372 |
|
|
373 Matrix& |
|
|
374 Matrix::insert (const ColumnVector& a, int r, int c) |
|
|
375 { |
|
|
376 if (r < 0 || r + a.len - 1 > nr || c < 0 || c >= nc) |
|
227
|
377 { |
|
|
378 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
379 return *this; |
|
|
380 } |
|
3
|
381 |
|
|
382 for (int i = 0; i < a.len; i++) |
|
|
383 elem (r+i, c) = a.data[i]; |
|
|
384 |
|
|
385 return *this; |
|
|
386 } |
|
|
387 |
|
|
388 Matrix& |
|
|
389 Matrix::insert (const DiagMatrix& a, int r, int c) |
|
|
390 { |
|
|
391 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
392 { |
|
|
393 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
394 return *this; |
|
|
395 } |
|
3
|
396 |
|
|
397 for (int i = 0; i < a.len; i++) |
|
|
398 elem (r+i, c+i) = a.data[i]; |
|
|
399 |
|
|
400 return *this; |
|
|
401 } |
|
|
402 |
|
|
403 Matrix& |
|
|
404 Matrix::fill (double val) |
|
|
405 { |
|
|
406 if (nr > 0 && nc > 0) |
|
|
407 copy (data, len, val); |
|
|
408 return *this; |
|
|
409 } |
|
|
410 |
|
|
411 Matrix& |
|
|
412 Matrix::fill (double val, int r1, int c1, int r2, int c2) |
|
|
413 { |
|
|
414 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 |
|
|
415 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) |
|
227
|
416 { |
|
|
417 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
418 return *this; |
|
|
419 } |
|
3
|
420 |
|
|
421 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
422 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
423 |
|
|
424 for (int j = c1; j <= c2; j++) |
|
|
425 for (int i = r1; i <= r2; i++) |
|
|
426 elem (i, j) = val; |
|
|
427 |
|
|
428 return *this; |
|
|
429 } |
|
|
430 |
|
|
431 Matrix |
|
|
432 Matrix::append (const Matrix& a) const |
|
|
433 { |
|
|
434 if (nr != a.nr) |
|
227
|
435 { |
|
|
436 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
437 return Matrix (); |
|
|
438 } |
|
3
|
439 |
|
|
440 int nc_insert = nc; |
|
|
441 Matrix retval (nr, nc + a.nc); |
|
|
442 retval.insert (*this, 0, 0); |
|
|
443 retval.insert (a, 0, nc_insert); |
|
227
|
444 return retval; |
|
3
|
445 } |
|
|
446 |
|
|
447 Matrix |
|
|
448 Matrix::append (const RowVector& a) const |
|
|
449 { |
|
|
450 if (nr != 1) |
|
227
|
451 { |
|
|
452 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
453 return Matrix (); |
|
|
454 } |
|
3
|
455 |
|
|
456 int nc_insert = nc; |
|
|
457 Matrix retval (nr, nc + a.len); |
|
|
458 retval.insert (*this, 0, 0); |
|
|
459 retval.insert (a, 0, nc_insert); |
|
|
460 return retval; |
|
|
461 } |
|
|
462 |
|
|
463 Matrix |
|
|
464 Matrix::append (const ColumnVector& a) const |
|
|
465 { |
|
|
466 if (nr != a.len) |
|
227
|
467 { |
|
|
468 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
469 return Matrix (); |
|
|
470 } |
|
3
|
471 |
|
|
472 int nc_insert = nc; |
|
|
473 Matrix retval (nr, nc + 1); |
|
|
474 retval.insert (*this, 0, 0); |
|
|
475 retval.insert (a, 0, nc_insert); |
|
|
476 return retval; |
|
|
477 } |
|
|
478 |
|
|
479 Matrix |
|
|
480 Matrix::append (const DiagMatrix& a) const |
|
|
481 { |
|
|
482 if (nr != a.nr) |
|
227
|
483 { |
|
|
484 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
485 return *this; |
|
|
486 } |
|
3
|
487 |
|
|
488 int nc_insert = nc; |
|
|
489 Matrix retval (nr, nc + a.nc); |
|
|
490 retval.insert (*this, 0, 0); |
|
|
491 retval.insert (a, 0, nc_insert); |
|
|
492 return retval; |
|
|
493 } |
|
|
494 |
|
|
495 Matrix |
|
|
496 Matrix::stack (const Matrix& a) const |
|
|
497 { |
|
|
498 if (nc != a.nc) |
|
227
|
499 { |
|
|
500 (*current_liboctave_error_handler) |
|
|
501 ("column dimension mismatch for stack"); |
|
|
502 return Matrix (); |
|
|
503 } |
|
3
|
504 |
|
|
505 int nr_insert = nr; |
|
|
506 Matrix retval (nr + a.nr, nc); |
|
|
507 retval.insert (*this, 0, 0); |
|
|
508 retval.insert (a, nr_insert, 0); |
|
|
509 return retval; |
|
|
510 } |
|
|
511 |
|
|
512 Matrix |
|
|
513 Matrix::stack (const RowVector& a) const |
|
|
514 { |
|
|
515 if (nc != a.len) |
|
227
|
516 { |
|
|
517 (*current_liboctave_error_handler) |
|
|
518 ("column dimension mismatch for stack"); |
|
|
519 return Matrix (); |
|
|
520 } |
|
3
|
521 |
|
|
522 int nr_insert = nr; |
|
|
523 Matrix retval (nr + 1, nc); |
|
|
524 retval.insert (*this, 0, 0); |
|
|
525 retval.insert (a, nr_insert, 0); |
|
|
526 return retval; |
|
|
527 } |
|
|
528 |
|
|
529 Matrix |
|
|
530 Matrix::stack (const ColumnVector& a) const |
|
|
531 { |
|
|
532 if (nc != 1) |
|
227
|
533 { |
|
|
534 (*current_liboctave_error_handler) |
|
|
535 ("column dimension mismatch for stack"); |
|
|
536 return Matrix (); |
|
|
537 } |
|
3
|
538 |
|
|
539 int nr_insert = nr; |
|
|
540 Matrix retval (nr + a.len, nc); |
|
|
541 retval.insert (*this, 0, 0); |
|
|
542 retval.insert (a, nr_insert, 0); |
|
|
543 return retval; |
|
|
544 } |
|
|
545 |
|
|
546 Matrix |
|
|
547 Matrix::stack (const DiagMatrix& a) const |
|
|
548 { |
|
|
549 if (nc != a.nc) |
|
227
|
550 { |
|
|
551 (*current_liboctave_error_handler) |
|
|
552 ("column dimension mismatch for stack"); |
|
|
553 return Matrix (); |
|
|
554 } |
|
3
|
555 |
|
|
556 int nr_insert = nr; |
|
|
557 Matrix retval (nr + a.nr, nc); |
|
|
558 retval.insert (*this, 0, 0); |
|
|
559 retval.insert (a, nr_insert, 0); |
|
|
560 return retval; |
|
|
561 } |
|
|
562 |
|
|
563 Matrix |
|
|
564 Matrix::transpose (void) const |
|
|
565 { |
|
88
|
566 Matrix result (nc, nr); |
|
3
|
567 if (len > 0) |
|
|
568 { |
|
|
569 for (int j = 0; j < nc; j++) |
|
|
570 for (int i = 0; i < nr; i++) |
|
|
571 result.data[nc*i+j] = data[nr*j+i]; |
|
|
572 } |
|
|
573 return result; |
|
|
574 } |
|
|
575 |
|
|
576 Matrix |
|
|
577 Matrix::extract (int r1, int c1, int r2, int c2) const |
|
|
578 { |
|
|
579 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
580 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
581 |
|
|
582 int new_r = r2 - r1 + 1; |
|
|
583 int new_c = c2 - c1 + 1; |
|
|
584 |
|
|
585 Matrix result (new_r, new_c); |
|
|
586 |
|
|
587 for (int j = 0; j < new_c; j++) |
|
|
588 for (int i = 0; i < new_r; i++) |
|
|
589 result.data[new_r*j+i] = elem (r1+i, c1+j); |
|
|
590 |
|
|
591 return result; |
|
|
592 } |
|
|
593 |
|
|
594 // extract row or column i. |
|
|
595 |
|
|
596 RowVector |
|
|
597 Matrix::row (int i) const |
|
|
598 { |
|
|
599 if (i < 0 || i >= nr) |
|
227
|
600 { |
|
|
601 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
602 return RowVector (); |
|
|
603 } |
|
3
|
604 |
|
|
605 RowVector retval (nc); |
|
|
606 for (int j = 0; j < nc; j++) |
|
|
607 retval.elem (j) = elem (i, j); |
|
|
608 |
|
|
609 return retval; |
|
|
610 } |
|
|
611 |
|
|
612 RowVector |
|
|
613 Matrix::row (char *s) const |
|
|
614 { |
|
|
615 if (s == (char *) NULL) |
|
227
|
616 { |
|
|
617 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
618 return RowVector (); |
|
|
619 } |
|
3
|
620 |
|
|
621 char c = *s; |
|
|
622 if (c == 'f' || c == 'F') |
|
|
623 return row (0); |
|
|
624 else if (c == 'l' || c == 'L') |
|
|
625 return row (nr - 1); |
|
|
626 else |
|
227
|
627 { |
|
|
628 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
629 return RowVector (); |
|
|
630 } |
|
3
|
631 } |
|
|
632 |
|
|
633 ColumnVector |
|
|
634 Matrix::column (int i) const |
|
|
635 { |
|
|
636 if (i < 0 || i >= nc) |
|
227
|
637 { |
|
|
638 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
639 return ColumnVector (); |
|
|
640 } |
|
3
|
641 |
|
|
642 ColumnVector retval (nr); |
|
|
643 for (int j = 0; j < nr; j++) |
|
|
644 retval.elem (j) = elem (j, i); |
|
|
645 |
|
|
646 return retval; |
|
|
647 } |
|
|
648 |
|
|
649 ColumnVector |
|
|
650 Matrix::column (char *s) const |
|
|
651 { |
|
|
652 if (s == (char *) NULL) |
|
227
|
653 { |
|
|
654 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
655 return ColumnVector (); |
|
|
656 } |
|
3
|
657 |
|
|
658 char c = *s; |
|
|
659 if (c == 'f' || c == 'F') |
|
|
660 return column (0); |
|
|
661 else if (c == 'l' || c == 'L') |
|
|
662 return column (nc - 1); |
|
|
663 else |
|
227
|
664 { |
|
|
665 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
666 return ColumnVector (); |
|
|
667 } |
|
3
|
668 } |
|
|
669 |
|
|
670 Matrix |
|
|
671 Matrix::inverse (int& info, double& rcond) const |
|
|
672 { |
|
227
|
673 if (nr != nc || nr == 0 || nc == 0) |
|
|
674 { |
|
|
675 (*current_liboctave_error_handler) ("inverse requires square matrix"); |
|
|
676 return Matrix (); |
|
|
677 } |
|
3
|
678 |
|
|
679 info = 0; |
|
|
680 |
|
|
681 int *ipvt = new int [nr]; |
|
|
682 double *z = new double [nr]; |
|
|
683 double *tmp_data = dup (data, len); |
|
|
684 |
|
|
685 F77_FCN (dgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z); |
|
|
686 |
|
|
687 if (rcond + 1.0 == 1.0) |
|
|
688 { |
|
|
689 info = -1; |
|
|
690 copy (tmp_data, data, len); // Restore matrix contents. |
|
|
691 } |
|
|
692 else |
|
|
693 { |
|
|
694 int job = 1; |
|
|
695 double dummy; |
|
|
696 |
|
|
697 F77_FCN (dgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job); |
|
|
698 } |
|
|
699 |
|
|
700 delete [] ipvt; |
|
|
701 delete [] z; |
|
|
702 |
|
|
703 return Matrix (tmp_data, nr, nc); |
|
|
704 } |
|
|
705 |
|
|
706 Matrix |
|
|
707 Matrix::inverse (int& info) const |
|
|
708 { |
|
|
709 double rcond; |
|
|
710 return inverse (info, rcond); |
|
|
711 } |
|
|
712 |
|
|
713 Matrix |
|
|
714 Matrix::inverse (void) const |
|
|
715 { |
|
|
716 int info; |
|
|
717 double rcond; |
|
|
718 return inverse (info, rcond); |
|
|
719 } |
|
|
720 |
|
|
721 ComplexMatrix |
|
|
722 Matrix::fourier (void) const |
|
|
723 { |
|
|
724 int npts, nsamples; |
|
|
725 if (nr == 1 || nc == 1) |
|
|
726 { |
|
|
727 npts = nr > nc ? nr : nc; |
|
|
728 nsamples = 1; |
|
|
729 } |
|
|
730 else |
|
|
731 { |
|
|
732 npts = nr; |
|
|
733 nsamples = nc; |
|
|
734 } |
|
|
735 |
|
|
736 int nn = 4*npts+15; |
|
|
737 Complex *wsave = new Complex [nn]; |
|
|
738 Complex *tmp_data = make_complex (data, len); |
|
|
739 |
|
|
740 F77_FCN (cffti) (&npts, wsave); |
|
|
741 |
|
|
742 for (int j = 0; j < nsamples; j++) |
|
|
743 F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave); |
|
|
744 |
|
|
745 delete [] wsave; |
|
|
746 |
|
|
747 return ComplexMatrix (tmp_data, nr, nc); |
|
|
748 } |
|
|
749 |
|
|
750 ComplexMatrix |
|
|
751 Matrix::ifourier (void) const |
|
|
752 { |
|
|
753 int npts, nsamples; |
|
|
754 if (nr == 1 || nc == 1) |
|
|
755 { |
|
|
756 npts = nr > nc ? nr : nc; |
|
|
757 nsamples = 1; |
|
|
758 } |
|
|
759 else |
|
|
760 { |
|
|
761 npts = nr; |
|
|
762 nsamples = nc; |
|
|
763 } |
|
|
764 |
|
|
765 int nn = 4*npts+15; |
|
|
766 Complex *wsave = new Complex [nn]; |
|
|
767 Complex *tmp_data = make_complex (data, len); |
|
|
768 |
|
|
769 F77_FCN (cffti) (&npts, wsave); |
|
|
770 |
|
|
771 for (int j = 0; j < nsamples; j++) |
|
|
772 F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave); |
|
|
773 |
|
|
774 for (j = 0; j < npts*nsamples; j++) |
|
|
775 tmp_data[j] = tmp_data[j] / (double) npts; |
|
|
776 |
|
|
777 delete [] wsave; |
|
|
778 |
|
|
779 return ComplexMatrix (tmp_data, nr, nc); |
|
|
780 } |
|
|
781 |
|
|
782 DET |
|
|
783 Matrix::determinant (void) const |
|
|
784 { |
|
|
785 int info; |
|
|
786 double rcond; |
|
|
787 return determinant (info, rcond); |
|
|
788 } |
|
|
789 |
|
|
790 DET |
|
|
791 Matrix::determinant (int& info) const |
|
|
792 { |
|
|
793 double rcond; |
|
|
794 return determinant (info, rcond); |
|
|
795 } |
|
|
796 |
|
|
797 DET |
|
|
798 Matrix::determinant (int& info, double& rcond) const |
|
|
799 { |
|
|
800 DET retval; |
|
|
801 |
|
|
802 if (nr == 0 || nc == 0) |
|
|
803 { |
|
|
804 double d[2]; |
|
|
805 d[0] = 1.0; |
|
|
806 d[1] = 0.0; |
|
|
807 return DET (d); |
|
|
808 } |
|
|
809 |
|
|
810 info = 0; |
|
|
811 int *ipvt = new int [nr]; |
|
|
812 |
|
|
813 double *z = new double [nr]; |
|
|
814 double *tmp_data = dup (data, len); |
|
|
815 |
|
|
816 F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
817 |
|
|
818 if (rcond + 1.0 == 1.0) |
|
|
819 { |
|
|
820 info = -1; |
|
|
821 } |
|
|
822 else |
|
|
823 { |
|
|
824 int job = 10; |
|
|
825 double d[2]; |
|
|
826 F77_FCN (dgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job); |
|
|
827 retval = DET (d); |
|
|
828 } |
|
|
829 |
|
|
830 delete [] tmp_data; |
|
|
831 delete [] ipvt; |
|
|
832 delete [] z; |
|
|
833 |
|
|
834 return retval; |
|
|
835 } |
|
|
836 |
|
|
837 Matrix |
|
|
838 Matrix::solve (const Matrix& b) const |
|
|
839 { |
|
|
840 int info; |
|
|
841 double rcond; |
|
|
842 return solve (b, info, rcond); |
|
|
843 } |
|
|
844 |
|
|
845 Matrix |
|
|
846 Matrix::solve (const Matrix& b, int& info) const |
|
|
847 { |
|
|
848 double rcond; |
|
|
849 return solve (b, info, rcond); |
|
|
850 } |
|
|
851 |
|
|
852 Matrix |
|
|
853 Matrix::solve (const Matrix& b, int& info, double& rcond) const |
|
|
854 { |
|
|
855 Matrix retval; |
|
|
856 |
|
|
857 if (nr == 0 || nc == 0 || nr != nc || nr != b.nr) |
|
227
|
858 { |
|
|
859 (*current_liboctave_error_handler) |
|
|
860 ("matrix dimension mismatch solution of linear equations"); |
|
|
861 return Matrix (); |
|
|
862 } |
|
3
|
863 |
|
|
864 info = 0; |
|
|
865 int *ipvt = new int [nr]; |
|
|
866 |
|
|
867 double *z = new double [nr]; |
|
|
868 double *tmp_data = dup (data, len); |
|
|
869 |
|
|
870 F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
871 |
|
|
872 if (rcond + 1.0 == 1.0) |
|
|
873 { |
|
|
874 info = -2; |
|
|
875 } |
|
|
876 else |
|
|
877 { |
|
|
878 int job = 0; |
|
|
879 |
|
|
880 double *result = dup (b.data, b.len); |
|
|
881 |
|
|
882 for (int j = 0; j < b.nc; j++) |
|
|
883 F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job); |
|
|
884 |
|
|
885 retval = Matrix (result, b.nr, b.nc); |
|
|
886 } |
|
|
887 |
|
|
888 delete [] tmp_data; |
|
|
889 delete [] ipvt; |
|
|
890 delete [] z; |
|
|
891 |
|
|
892 return retval; |
|
|
893 } |
|
|
894 |
|
|
895 ComplexMatrix |
|
|
896 Matrix::solve (const ComplexMatrix& b) const |
|
|
897 { |
|
|
898 ComplexMatrix tmp (*this); |
|
|
899 return tmp.solve (b); |
|
|
900 } |
|
|
901 |
|
|
902 ComplexMatrix |
|
|
903 Matrix::solve (const ComplexMatrix& b, int& info) const |
|
|
904 { |
|
|
905 ComplexMatrix tmp (*this); |
|
|
906 return tmp.solve (b, info); |
|
|
907 } |
|
|
908 |
|
|
909 ComplexMatrix |
|
|
910 Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const |
|
|
911 { |
|
|
912 ComplexMatrix tmp (*this); |
|
|
913 return tmp.solve (b, info, rcond); |
|
|
914 } |
|
|
915 |
|
|
916 ColumnVector |
|
|
917 Matrix::solve (const ColumnVector& b) const |
|
|
918 { |
|
|
919 int info; |
|
|
920 double rcond; |
|
|
921 return solve (b, info, rcond); |
|
|
922 } |
|
|
923 |
|
|
924 ColumnVector |
|
|
925 Matrix::solve (const ColumnVector& b, int& info) const |
|
|
926 { |
|
|
927 double rcond; |
|
|
928 return solve (b, info, rcond); |
|
|
929 } |
|
|
930 |
|
|
931 ColumnVector |
|
|
932 Matrix::solve (const ColumnVector& b, int& info, double& rcond) const |
|
|
933 { |
|
|
934 ColumnVector retval; |
|
|
935 |
|
|
936 if (nr == 0 || nc == 0 || nr != nc || nr != b.len) |
|
227
|
937 { |
|
|
938 (*current_liboctave_error_handler) |
|
|
939 ("matrix dimension mismatch solution of linear equations"); |
|
|
940 return ColumnVector (); |
|
|
941 } |
|
3
|
942 |
|
|
943 info = 0; |
|
|
944 int *ipvt = new int [nr]; |
|
|
945 |
|
|
946 double *z = new double [nr]; |
|
|
947 double *tmp_data = dup (data, len); |
|
|
948 |
|
|
949 F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
950 |
|
|
951 if (rcond + 1.0 == 1.0) |
|
|
952 { |
|
|
953 info = -2; |
|
|
954 } |
|
|
955 else |
|
|
956 { |
|
|
957 int job = 0; |
|
|
958 |
|
|
959 double *result = dup (b.data, b.len); |
|
|
960 |
|
|
961 F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, result, &job); |
|
|
962 |
|
|
963 retval = ColumnVector (result, b.len); |
|
|
964 } |
|
|
965 |
|
|
966 delete [] tmp_data; |
|
|
967 delete [] ipvt; |
|
|
968 delete [] z; |
|
|
969 |
|
|
970 return retval; |
|
|
971 } |
|
|
972 |
|
|
973 ComplexColumnVector |
|
|
974 Matrix::solve (const ComplexColumnVector& b) const |
|
|
975 { |
|
|
976 ComplexMatrix tmp (*this); |
|
|
977 return tmp.solve (b); |
|
|
978 } |
|
|
979 |
|
|
980 ComplexColumnVector |
|
|
981 Matrix::solve (const ComplexColumnVector& b, int& info) const |
|
|
982 { |
|
|
983 ComplexMatrix tmp (*this); |
|
|
984 return tmp.solve (b, info); |
|
|
985 } |
|
|
986 |
|
|
987 ComplexColumnVector |
|
|
988 Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const |
|
|
989 { |
|
|
990 ComplexMatrix tmp (*this); |
|
|
991 return tmp.solve (b, info, rcond); |
|
|
992 } |
|
|
993 |
|
|
994 Matrix |
|
|
995 Matrix::lssolve (const Matrix& b) const |
|
|
996 { |
|
|
997 int info; |
|
|
998 int rank; |
|
|
999 return lssolve (b, info, rank); |
|
|
1000 } |
|
|
1001 |
|
|
1002 Matrix |
|
|
1003 Matrix::lssolve (const Matrix& b, int& info) const |
|
|
1004 { |
|
|
1005 int rank; |
|
|
1006 return lssolve (b, info, rank); |
|
|
1007 } |
|
|
1008 |
|
|
1009 Matrix |
|
|
1010 Matrix::lssolve (const Matrix& b, int& info, int& rank) const |
|
|
1011 { |
|
|
1012 int nrhs = b.nc; |
|
|
1013 |
|
|
1014 int m = nr; |
|
|
1015 int n = nc; |
|
|
1016 |
|
|
1017 if (m == 0 || n == 0 || m != b.nr) |
|
227
|
1018 { |
|
|
1019 (*current_liboctave_error_handler) |
|
|
1020 ("matrix dimension mismatch in solution of least squares problem"); |
|
|
1021 return Matrix (); |
|
|
1022 } |
|
3
|
1023 |
|
|
1024 double *tmp_data = dup (data, len); |
|
|
1025 |
|
|
1026 int nrr = m > n ? m : n; |
|
|
1027 Matrix result (nrr, nrhs); |
|
|
1028 |
|
|
1029 int i, j; |
|
|
1030 for (j = 0; j < nrhs; j++) |
|
|
1031 for (i = 0; i < m; i++) |
|
|
1032 result.elem (i, j) = b.elem (i, j); |
|
|
1033 |
|
|
1034 double *presult = result.fortran_vec (); |
|
|
1035 |
|
|
1036 int len_s = m < n ? m : n; |
|
|
1037 double *s = new double [len_s]; |
|
|
1038 double rcond = -1.0; |
|
|
1039 int lwork; |
|
|
1040 if (m < n) |
|
|
1041 lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n)); |
|
|
1042 else |
|
|
1043 lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m)); |
|
|
1044 |
|
|
1045 double *work = new double [lwork]; |
|
|
1046 |
|
|
1047 F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, |
|
|
1048 &rcond, &rank, work, &lwork, &info); |
|
|
1049 |
|
|
1050 Matrix retval (n, nrhs); |
|
|
1051 for (j = 0; j < nrhs; j++) |
|
|
1052 for (i = 0; i < n; i++) |
|
|
1053 retval.elem (i, j) = result.elem (i, j); |
|
|
1054 |
|
|
1055 delete [] tmp_data; |
|
|
1056 delete [] s; |
|
|
1057 delete [] work; |
|
|
1058 |
|
|
1059 return retval; |
|
|
1060 } |
|
|
1061 |
|
|
1062 ComplexMatrix |
|
|
1063 Matrix::lssolve (const ComplexMatrix& b) const |
|
|
1064 { |
|
|
1065 ComplexMatrix tmp (*this); |
|
|
1066 return tmp.lssolve (b); |
|
|
1067 } |
|
|
1068 |
|
|
1069 ComplexMatrix |
|
|
1070 Matrix::lssolve (const ComplexMatrix& b, int& info) const |
|
|
1071 { |
|
|
1072 ComplexMatrix tmp (*this); |
|
|
1073 return tmp.lssolve (b, info); |
|
|
1074 } |
|
|
1075 |
|
|
1076 ComplexMatrix |
|
|
1077 Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const |
|
|
1078 { |
|
|
1079 ComplexMatrix tmp (*this); |
|
|
1080 return tmp.lssolve (b, info, rank); |
|
|
1081 } |
|
|
1082 |
|
|
1083 ColumnVector |
|
|
1084 Matrix::lssolve (const ColumnVector& b) const |
|
|
1085 { |
|
|
1086 int info; |
|
|
1087 int rank; |
|
|
1088 return lssolve (b, info, rank); |
|
|
1089 } |
|
|
1090 |
|
|
1091 ColumnVector |
|
|
1092 Matrix::lssolve (const ColumnVector& b, int& info) const |
|
|
1093 { |
|
|
1094 int rank; |
|
|
1095 return lssolve (b, info, rank); |
|
|
1096 } |
|
|
1097 |
|
|
1098 ColumnVector |
|
|
1099 Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const |
|
|
1100 { |
|
|
1101 int nrhs = 1; |
|
|
1102 |
|
|
1103 int m = nr; |
|
|
1104 int n = nc; |
|
|
1105 |
|
|
1106 if (m == 0 || n == 0 || m != b.len) |
|
227
|
1107 { |
|
|
1108 (*current_liboctave_error_handler) |
|
|
1109 ("matrix dimension mismatch in solution of least squares problem"); |
|
|
1110 return ColumnVector (); |
|
|
1111 } |
|
3
|
1112 |
|
|
1113 double *tmp_data = dup (data, len); |
|
|
1114 |
|
|
1115 int nrr = m > n ? m : n; |
|
|
1116 ColumnVector result (nrr); |
|
|
1117 |
|
|
1118 int i; |
|
|
1119 for (i = 0; i < m; i++) |
|
|
1120 result.elem (i) = b.elem (i); |
|
|
1121 |
|
|
1122 double *presult = result.fortran_vec (); |
|
|
1123 |
|
|
1124 int len_s = m < n ? m : n; |
|
|
1125 double *s = new double [len_s]; |
|
|
1126 double rcond = -1.0; |
|
|
1127 int lwork; |
|
|
1128 if (m < n) |
|
|
1129 lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n)); |
|
|
1130 else |
|
|
1131 lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m)); |
|
|
1132 |
|
|
1133 double *work = new double [lwork]; |
|
|
1134 |
|
|
1135 F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, |
|
|
1136 &rcond, &rank, work, &lwork, &info); |
|
|
1137 |
|
|
1138 ColumnVector retval (n); |
|
|
1139 for (i = 0; i < n; i++) |
|
|
1140 retval.elem (i) = result.elem (i); |
|
|
1141 |
|
|
1142 delete [] tmp_data; |
|
|
1143 delete [] s; |
|
|
1144 delete [] work; |
|
|
1145 |
|
|
1146 return retval; |
|
|
1147 } |
|
|
1148 |
|
|
1149 ComplexColumnVector |
|
|
1150 Matrix::lssolve (const ComplexColumnVector& b) const |
|
|
1151 { |
|
|
1152 ComplexMatrix tmp (*this); |
|
|
1153 return tmp.lssolve (b); |
|
|
1154 } |
|
|
1155 |
|
|
1156 ComplexColumnVector |
|
|
1157 Matrix::lssolve (const ComplexColumnVector& b, int& info) const |
|
|
1158 { |
|
|
1159 ComplexMatrix tmp (*this); |
|
|
1160 return tmp.lssolve (b, info); |
|
|
1161 } |
|
|
1162 |
|
|
1163 ComplexColumnVector |
|
|
1164 Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const |
|
|
1165 { |
|
|
1166 ComplexMatrix tmp (*this); |
|
|
1167 return tmp.lssolve (b, info, rank); |
|
|
1168 } |
|
|
1169 |
|
|
1170 // matrix by scalar -> matrix operations. |
|
|
1171 |
|
|
1172 Matrix |
|
|
1173 Matrix::operator + (double s) const |
|
|
1174 { |
|
|
1175 return Matrix (add (data, len, s), nr, nc); |
|
|
1176 } |
|
|
1177 |
|
|
1178 Matrix |
|
|
1179 Matrix::operator - (double s) const |
|
|
1180 { |
|
|
1181 return Matrix (subtract (data, len, s), nr, nc); |
|
|
1182 } |
|
|
1183 |
|
|
1184 Matrix |
|
|
1185 Matrix::operator * (double s) const |
|
|
1186 { |
|
|
1187 return Matrix (multiply (data, len, s), nr, nc); |
|
|
1188 } |
|
|
1189 |
|
|
1190 Matrix |
|
|
1191 Matrix::operator / (double s) const |
|
|
1192 { |
|
|
1193 return Matrix (divide (data, len, s), nr, nc); |
|
|
1194 } |
|
|
1195 |
|
|
1196 ComplexMatrix |
|
161
|
1197 Matrix::operator + (const Complex& s) const |
|
3
|
1198 { |
|
|
1199 return ComplexMatrix (add (data, len, s), nr, nc); |
|
|
1200 } |
|
|
1201 |
|
|
1202 ComplexMatrix |
|
161
|
1203 Matrix::operator - (const Complex& s) const |
|
3
|
1204 { |
|
|
1205 return ComplexMatrix (subtract (data, len, s), nr, nc); |
|
|
1206 } |
|
|
1207 |
|
|
1208 ComplexMatrix |
|
161
|
1209 Matrix::operator * (const Complex& s) const |
|
3
|
1210 { |
|
|
1211 return ComplexMatrix (multiply (data, len, s), nr, nc); |
|
|
1212 } |
|
|
1213 |
|
|
1214 ComplexMatrix |
|
161
|
1215 Matrix::operator / (const Complex& s) const |
|
3
|
1216 { |
|
|
1217 return ComplexMatrix (divide (data, len, s), nr, nc); |
|
|
1218 } |
|
|
1219 |
|
|
1220 // scalar by matrix -> matrix operations |
|
|
1221 |
|
|
1222 Matrix |
|
|
1223 operator + (double s, const Matrix& a) |
|
|
1224 { |
|
|
1225 return Matrix (add (a.data, a.len, s), a.nr, a.nc); |
|
|
1226 } |
|
|
1227 |
|
|
1228 Matrix |
|
|
1229 operator - (double s, const Matrix& a) |
|
|
1230 { |
|
|
1231 return Matrix (subtract (s, a.data, a.len), a.nr, a.nc); |
|
|
1232 } |
|
|
1233 |
|
|
1234 Matrix |
|
|
1235 operator * (double s, const Matrix& a) |
|
|
1236 { |
|
|
1237 return Matrix (multiply (a.data, a.len, s), a.nr, a.nc); |
|
|
1238 } |
|
|
1239 |
|
|
1240 Matrix |
|
|
1241 operator / (double s, const Matrix& a) |
|
|
1242 { |
|
|
1243 return Matrix (divide (s, a.data, a.len), a.nr, a.nc); |
|
|
1244 } |
|
|
1245 |
|
|
1246 // matrix by column vector -> column vector operations |
|
|
1247 |
|
|
1248 ColumnVector |
|
|
1249 Matrix::operator * (const ColumnVector& a) const |
|
|
1250 { |
|
|
1251 if (nc != a.len) |
|
227
|
1252 { |
|
|
1253 (*current_liboctave_error_handler) |
|
|
1254 ("nonconformant matrix multiplication attempted"); |
|
|
1255 return ColumnVector (); |
|
|
1256 } |
|
3
|
1257 |
|
|
1258 if (nr == 0 || nc == 0) |
|
|
1259 return ColumnVector (0); |
|
|
1260 |
|
|
1261 char trans = 'N'; |
|
|
1262 int ld = nr; |
|
|
1263 double alpha = 1.0; |
|
|
1264 double beta = 0.0; |
|
|
1265 int i_one = 1; |
|
|
1266 |
|
125
|
1267 double *y = new double [nr]; |
|
3
|
1268 |
|
|
1269 F77_FCN (dgemv) (&trans, &nr, &nc, &alpha, data, &ld, a.data, |
|
|
1270 &i_one, &beta, y, &i_one, 1L); |
|
|
1271 |
|
124
|
1272 return ColumnVector (y, nr); |
|
3
|
1273 } |
|
|
1274 |
|
|
1275 ComplexColumnVector |
|
|
1276 Matrix::operator * (const ComplexColumnVector& a) const |
|
|
1277 { |
|
|
1278 ComplexMatrix tmp (*this); |
|
|
1279 return tmp * a; |
|
|
1280 } |
|
|
1281 |
|
|
1282 // matrix by diagonal matrix -> matrix operations |
|
|
1283 |
|
|
1284 Matrix |
|
|
1285 Matrix::operator + (const DiagMatrix& a) const |
|
|
1286 { |
|
|
1287 if (nr != a.nr || nc != a.nc) |
|
227
|
1288 { |
|
|
1289 (*current_liboctave_error_handler) |
|
|
1290 ("nonconformant matrix addition attempted"); |
|
|
1291 return Matrix (); |
|
|
1292 } |
|
3
|
1293 |
|
|
1294 if (nr == 0 || nc == 0) |
|
|
1295 return Matrix (nr, nc); |
|
|
1296 |
|
|
1297 Matrix result (*this); |
|
|
1298 for (int i = 0; i < a.len; i++) |
|
|
1299 result.elem (i, i) += a.data[i]; |
|
|
1300 |
|
|
1301 return result; |
|
|
1302 } |
|
|
1303 |
|
|
1304 Matrix |
|
|
1305 Matrix::operator - (const DiagMatrix& a) const |
|
|
1306 { |
|
|
1307 if (nr != a.nr || nc != a.nc) |
|
227
|
1308 { |
|
|
1309 (*current_liboctave_error_handler) |
|
|
1310 ("nonconformant matrix subtraction attempted"); |
|
|
1311 return Matrix (); |
|
|
1312 } |
|
3
|
1313 |
|
|
1314 if (nr == 0 || nc == 0) |
|
|
1315 return Matrix (nr, nc); |
|
|
1316 |
|
|
1317 Matrix result (*this); |
|
|
1318 for (int i = 0; i < a.len; i++) |
|
|
1319 result.elem (i, i) -= a.data[i]; |
|
|
1320 |
|
|
1321 return result; |
|
|
1322 } |
|
|
1323 |
|
|
1324 Matrix |
|
|
1325 Matrix::operator * (const DiagMatrix& a) const |
|
|
1326 { |
|
|
1327 if (nc != a.nr) |
|
227
|
1328 { |
|
|
1329 (*current_liboctave_error_handler) |
|
|
1330 ("nonconformant matrix multiplication attempted"); |
|
|
1331 return Matrix (); |
|
|
1332 } |
|
3
|
1333 |
|
|
1334 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
1335 return Matrix (nr, a.nc, 0.0); |
|
|
1336 |
|
|
1337 double *c = new double [nr*a.nc]; |
|
|
1338 double *ctmp = (double *) NULL; |
|
|
1339 |
|
|
1340 for (int j = 0; j < a.len; j++) |
|
|
1341 { |
|
|
1342 int idx = j * nr; |
|
|
1343 ctmp = c + idx; |
|
|
1344 if (a.data[j] == 1.0) |
|
|
1345 { |
|
|
1346 for (int i = 0; i < nr; i++) |
|
|
1347 ctmp[i] = elem (i, j); |
|
|
1348 } |
|
|
1349 else if (a.data[j] == 0.0) |
|
|
1350 { |
|
|
1351 for (int i = 0; i < nr; i++) |
|
|
1352 ctmp[i] = 0.0; |
|
|
1353 } |
|
|
1354 else |
|
|
1355 { |
|
|
1356 for (int i = 0; i < nr; i++) |
|
|
1357 ctmp[i] = a.data[j] * elem (i, j); |
|
|
1358 } |
|
|
1359 } |
|
|
1360 |
|
|
1361 if (a.nr < a.nc) |
|
|
1362 { |
|
|
1363 for (int i = nr * nc; i < nr * a.nc; i++) |
|
|
1364 ctmp[i] = 0.0; |
|
|
1365 } |
|
|
1366 |
|
|
1367 return Matrix (c, nr, a.nc); |
|
|
1368 } |
|
|
1369 |
|
|
1370 ComplexMatrix |
|
|
1371 Matrix::operator + (const ComplexDiagMatrix& a) const |
|
|
1372 { |
|
|
1373 if (nr != a.nr || nc != a.nc) |
|
227
|
1374 { |
|
|
1375 (*current_liboctave_error_handler) |
|
|
1376 ("nonconformant matrix addition attempted"); |
|
|
1377 return ComplexMatrix (); |
|
|
1378 } |
|
3
|
1379 |
|
|
1380 if (nr == 0 || nc == 0) |
|
|
1381 return ComplexMatrix (nr, nc); |
|
|
1382 |
|
|
1383 ComplexMatrix result (*this); |
|
|
1384 for (int i = 0; i < a.len; i++) |
|
|
1385 result.elem (i, i) += a.data[i]; |
|
|
1386 |
|
|
1387 return result; |
|
|
1388 } |
|
|
1389 |
|
|
1390 ComplexMatrix |
|
|
1391 Matrix::operator - (const ComplexDiagMatrix& a) const |
|
|
1392 { |
|
|
1393 if (nr != a.nr || nc != a.nc) |
|
227
|
1394 { |
|
|
1395 (*current_liboctave_error_handler) |
|
|
1396 ("nonconformant matrix subtraction attempted"); |
|
|
1397 return ComplexMatrix (); |
|
|
1398 } |
|
3
|
1399 |
|
|
1400 if (nr == 0 || nc == 0) |
|
|
1401 return ComplexMatrix (nr, nc); |
|
|
1402 |
|
|
1403 ComplexMatrix result (*this); |
|
|
1404 for (int i = 0; i < a.len; i++) |
|
|
1405 result.elem (i, i) -= a.data[i]; |
|
|
1406 |
|
|
1407 return result; |
|
|
1408 } |
|
|
1409 |
|
|
1410 ComplexMatrix |
|
|
1411 Matrix::operator * (const ComplexDiagMatrix& a) const |
|
|
1412 { |
|
|
1413 if (nc != a.nr) |
|
227
|
1414 { |
|
|
1415 (*current_liboctave_error_handler) |
|
|
1416 ("nonconformant matrix multiplication attempted"); |
|
|
1417 return ComplexMatrix (); |
|
|
1418 } |
|
3
|
1419 |
|
|
1420 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
1421 return ComplexMatrix (nr, a.nc, 0.0); |
|
|
1422 |
|
|
1423 Complex *c = new Complex [nr*a.nc]; |
|
|
1424 Complex *ctmp = (Complex *) NULL; |
|
|
1425 |
|
|
1426 for (int j = 0; j < a.len; j++) |
|
|
1427 { |
|
|
1428 int idx = j * nr; |
|
|
1429 ctmp = c + idx; |
|
|
1430 if (a.data[j] == 1.0) |
|
|
1431 { |
|
|
1432 for (int i = 0; i < nr; i++) |
|
|
1433 ctmp[i] = elem (i, j); |
|
|
1434 } |
|
|
1435 else if (a.data[j] == 0.0) |
|
|
1436 { |
|
|
1437 for (int i = 0; i < nr; i++) |
|
|
1438 ctmp[i] = 0.0; |
|
|
1439 } |
|
|
1440 else |
|
|
1441 { |
|
|
1442 for (int i = 0; i < nr; i++) |
|
|
1443 ctmp[i] = a.data[j] * elem (i, j); |
|
|
1444 } |
|
|
1445 } |
|
|
1446 |
|
|
1447 if (a.nr < a.nc) |
|
|
1448 { |
|
|
1449 for (int i = nr * nc; i < nr * a.nc; i++) |
|
|
1450 ctmp[i] = 0.0; |
|
|
1451 } |
|
|
1452 |
|
|
1453 return ComplexMatrix (c, nr, a.nc); |
|
|
1454 } |
|
|
1455 |
|
|
1456 Matrix& |
|
|
1457 Matrix::operator += (const DiagMatrix& a) |
|
|
1458 { |
|
|
1459 if (nr != a.nr || nc != a.nc) |
|
227
|
1460 { |
|
|
1461 (*current_liboctave_error_handler) |
|
|
1462 ("nonconformant matrix += operation attempted"); |
|
|
1463 return *this; |
|
|
1464 } |
|
3
|
1465 |
|
|
1466 for (int i = 0; i < a.len; i++) |
|
|
1467 elem (i, i) += a.data[i]; |
|
|
1468 |
|
|
1469 return *this; |
|
|
1470 } |
|
|
1471 |
|
|
1472 Matrix& |
|
|
1473 Matrix::operator -= (const DiagMatrix& a) |
|
|
1474 { |
|
|
1475 if (nr != a.nr || nc != a.nc) |
|
227
|
1476 { |
|
|
1477 (*current_liboctave_error_handler) |
|
|
1478 ("nonconformant matrix += operation attempted"); |
|
|
1479 return *this; |
|
|
1480 } |
|
3
|
1481 |
|
|
1482 for (int i = 0; i < a.len; i++) |
|
|
1483 elem (i, i) -= a.data[i]; |
|
|
1484 |
|
|
1485 return *this; |
|
|
1486 } |
|
|
1487 |
|
|
1488 // matrix by matrix -> matrix operations |
|
|
1489 |
|
|
1490 Matrix |
|
|
1491 Matrix::operator + (const Matrix& a) const |
|
|
1492 { |
|
|
1493 if (nr != a.nr || nc != a.nc) |
|
227
|
1494 { |
|
|
1495 (*current_liboctave_error_handler) |
|
|
1496 ("nonconformant matrix addition attempted"); |
|
|
1497 return Matrix (); |
|
|
1498 } |
|
3
|
1499 |
|
|
1500 if (nr == 0 || nc == 0) |
|
|
1501 return Matrix (nr, nc); |
|
|
1502 |
|
|
1503 return Matrix (add (data, a.data, len), nr, nc); |
|
|
1504 } |
|
|
1505 |
|
|
1506 Matrix |
|
|
1507 Matrix::operator - (const Matrix& a) const |
|
|
1508 { |
|
|
1509 if (nr != a.nr || nc != a.nc) |
|
227
|
1510 { |
|
|
1511 (*current_liboctave_error_handler) |
|
|
1512 ("nonconformant matrix subtraction attempted"); |
|
|
1513 return Matrix (); |
|
|
1514 } |
|
3
|
1515 |
|
|
1516 if (nr == 0 || nc == 0) |
|
|
1517 return Matrix (nr, nc); |
|
|
1518 |
|
|
1519 return Matrix (subtract (data, a.data, len), nr, nc); |
|
|
1520 } |
|
|
1521 |
|
|
1522 Matrix |
|
|
1523 Matrix::operator * (const Matrix& a) const |
|
|
1524 { |
|
|
1525 if (nc != a.nr) |
|
227
|
1526 { |
|
|
1527 (*current_liboctave_error_handler) |
|
|
1528 ("nonconformant matrix multiplication attempted"); |
|
|
1529 return Matrix (); |
|
|
1530 } |
|
3
|
1531 |
|
|
1532 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
1533 return Matrix (nr, a.nc, 0.0); |
|
|
1534 |
|
|
1535 char trans = 'N'; |
|
|
1536 char transa = 'N'; |
|
|
1537 |
|
|
1538 int ld = nr; |
|
|
1539 int lda = a.nr; |
|
|
1540 |
|
|
1541 double alpha = 1.0; |
|
|
1542 double beta = 0.0; |
|
|
1543 int anc = a.nc; |
|
|
1544 |
|
|
1545 double *c = new double [nr*a.nc]; |
|
|
1546 |
|
|
1547 F77_FCN (dgemm) (&trans, &transa, &nr, &anc, &nc, &alpha, data, &ld, |
|
|
1548 a.data, &lda, &beta, c, &nr, 1L, 1L); |
|
|
1549 |
|
|
1550 return Matrix (c, nr, a.nc); |
|
|
1551 } |
|
|
1552 |
|
|
1553 ComplexMatrix |
|
|
1554 Matrix::operator + (const ComplexMatrix& a) const |
|
|
1555 { |
|
|
1556 if (nr != a.nr || nc != a.nc) |
|
227
|
1557 { |
|
|
1558 (*current_liboctave_error_handler) |
|
|
1559 ("nonconformant matrix addition attempted"); |
|
|
1560 return ComplexMatrix (); |
|
|
1561 } |
|
3
|
1562 |
|
|
1563 return ComplexMatrix (add (data, a.data, len), nr, nc); |
|
|
1564 } |
|
|
1565 |
|
|
1566 ComplexMatrix |
|
|
1567 Matrix::operator - (const ComplexMatrix& a) const |
|
|
1568 { |
|
|
1569 if (nr != a.nr || nc != a.nc) |
|
227
|
1570 { |
|
|
1571 (*current_liboctave_error_handler) |
|
|
1572 ("nonconformant matrix subtraction attempted"); |
|
|
1573 return ComplexMatrix (); |
|
|
1574 } |
|
3
|
1575 |
|
|
1576 if (nr == 0 || nc == 0) |
|
|
1577 return ComplexMatrix (nr, nc); |
|
|
1578 |
|
|
1579 return ComplexMatrix (subtract (data, a.data, len), nr, nc); |
|
|
1580 } |
|
|
1581 |
|
|
1582 ComplexMatrix |
|
|
1583 Matrix::operator * (const ComplexMatrix& a) const |
|
|
1584 { |
|
|
1585 ComplexMatrix tmp (*this); |
|
|
1586 return tmp * a; |
|
|
1587 } |
|
|
1588 |
|
|
1589 Matrix |
|
|
1590 Matrix::product (const Matrix& a) const |
|
|
1591 { |
|
|
1592 if (nr != a.nr || nc != a.nc) |
|
227
|
1593 { |
|
|
1594 (*current_liboctave_error_handler) |
|
|
1595 ("nonconformant matrix product attempted"); |
|
|
1596 return Matrix (); |
|
|
1597 } |
|
3
|
1598 |
|
|
1599 if (nr == 0 || nc == 0) |
|
|
1600 return Matrix (nr, nc); |
|
|
1601 |
|
|
1602 return Matrix (multiply (data, a.data, len), nr, nc); |
|
|
1603 } |
|
|
1604 |
|
|
1605 Matrix |
|
|
1606 Matrix::quotient (const Matrix& a) const |
|
|
1607 { |
|
|
1608 if (nr != a.nr || nc != a.nc) |
|
227
|
1609 { |
|
|
1610 (*current_liboctave_error_handler) |
|
|
1611 ("nonconformant matrix quotient attempted"); |
|
|
1612 return Matrix (); |
|
|
1613 } |
|
3
|
1614 |
|
|
1615 if (nr == 0 || nc == 0) |
|
|
1616 return Matrix (nr, nc); |
|
|
1617 |
|
|
1618 return Matrix (divide (data, a.data, len), nr, nc); |
|
|
1619 } |
|
|
1620 |
|
|
1621 ComplexMatrix |
|
|
1622 Matrix::product (const ComplexMatrix& a) const |
|
|
1623 { |
|
|
1624 if (nr != a.nr || nc != a.nc) |
|
227
|
1625 { |
|
|
1626 (*current_liboctave_error_handler) |
|
|
1627 ("nonconformant matrix product attempted"); |
|
|
1628 return ComplexMatrix (); |
|
|
1629 } |
|
3
|
1630 |
|
|
1631 if (nr == 0 || nc == 0) |
|
|
1632 return ComplexMatrix (nr, nc); |
|
|
1633 |
|
|
1634 return ComplexMatrix (multiply (data, a.data, len), nr, nc); |
|
|
1635 } |
|
|
1636 |
|
|
1637 ComplexMatrix |
|
|
1638 Matrix::quotient (const ComplexMatrix& a) const |
|
|
1639 { |
|
|
1640 if (nr != a.nr || nc != a.nc) |
|
227
|
1641 { |
|
|
1642 (*current_liboctave_error_handler) |
|
|
1643 ("nonconformant matrix quotient attempted"); |
|
|
1644 return ComplexMatrix (); |
|
|
1645 } |
|
3
|
1646 |
|
|
1647 if (nr == 0 || nc == 0) |
|
|
1648 return ComplexMatrix (nr, nc); |
|
|
1649 |
|
|
1650 return ComplexMatrix (divide (data, a.data, len), nr, nc); |
|
|
1651 } |
|
|
1652 |
|
|
1653 Matrix& |
|
|
1654 Matrix::operator += (const Matrix& a) |
|
|
1655 { |
|
|
1656 if (nr != a.nr || nc != a.nc) |
|
227
|
1657 { |
|
|
1658 (*current_liboctave_error_handler) |
|
|
1659 ("nonconformant matrix += operation attempted"); |
|
|
1660 return *this; |
|
|
1661 } |
|
3
|
1662 |
|
|
1663 if (nr == 0 || nc == 0) |
|
|
1664 return *this; |
|
|
1665 |
|
|
1666 add2 (data, a.data, len); |
|
|
1667 return *this; |
|
|
1668 } |
|
|
1669 |
|
|
1670 Matrix& |
|
|
1671 Matrix::operator -= (const Matrix& a) |
|
|
1672 { |
|
|
1673 if (nr != a.nr || nc != a.nc) |
|
227
|
1674 { |
|
|
1675 (*current_liboctave_error_handler) |
|
|
1676 ("nonconformant matrix -= operation attempted"); |
|
|
1677 return *this; |
|
|
1678 } |
|
3
|
1679 |
|
|
1680 if (nr == 0 || nc == 0) |
|
|
1681 return *this; |
|
|
1682 |
|
|
1683 subtract2 (data, a.data, len); |
|
|
1684 return *this; |
|
|
1685 } |
|
|
1686 |
|
|
1687 // other operations. |
|
|
1688 |
|
|
1689 Matrix |
|
|
1690 map (d_d_Mapper f, const Matrix& a) |
|
|
1691 { |
|
|
1692 Matrix b (a); |
|
|
1693 b.map (f); |
|
|
1694 return b; |
|
|
1695 } |
|
|
1696 |
|
|
1697 void |
|
|
1698 Matrix::map (d_d_Mapper f) |
|
|
1699 { |
|
|
1700 for (int i = 0; i < len; i++) |
|
|
1701 data[i] = f (data[i]); |
|
|
1702 } |
|
|
1703 |
|
|
1704 // XXX FIXME XXX Do these really belong here? They should maybe be |
|
|
1705 // cleaned up a bit, no? What about corresponding functions for the |
|
|
1706 // Vectors? |
|
|
1707 |
|
|
1708 Matrix |
|
|
1709 Matrix::all (void) const |
|
|
1710 { |
|
|
1711 Matrix retval; |
|
|
1712 if (nr > 0 && nc > 0) |
|
|
1713 { |
|
|
1714 if (nr == 1) |
|
|
1715 { |
|
|
1716 retval.resize (1, 1); |
|
|
1717 retval.elem (0, 0) = 1.0; |
|
|
1718 for (int j = 0; j < nc; j++) |
|
|
1719 { |
|
|
1720 if (elem (0, j) == 0.0) |
|
|
1721 { |
|
|
1722 retval.elem (0, 0) = 0.0; |
|
|
1723 break; |
|
|
1724 } |
|
|
1725 } |
|
|
1726 } |
|
|
1727 else if (nc == 1) |
|
|
1728 { |
|
|
1729 retval.resize (1, 1); |
|
|
1730 retval.elem (0, 0) = 1.0; |
|
|
1731 for (int i = 0; i < nr; i++) |
|
|
1732 { |
|
|
1733 if (elem (i, 0) == 0.0) |
|
|
1734 { |
|
|
1735 retval.elem (0, 0) = 0.0; |
|
|
1736 break; |
|
|
1737 } |
|
|
1738 } |
|
|
1739 } |
|
|
1740 else |
|
|
1741 { |
|
|
1742 retval.resize (1, nc); |
|
|
1743 for (int j = 0; j < nc; j++) |
|
|
1744 { |
|
|
1745 retval.elem (0, j) = 1.0; |
|
|
1746 for (int i = 0; i < nr; i++) |
|
|
1747 { |
|
|
1748 if (elem (i, j) == 0.0) |
|
|
1749 { |
|
|
1750 retval.elem (0, j) = 0.0; |
|
|
1751 break; |
|
|
1752 } |
|
|
1753 } |
|
|
1754 } |
|
|
1755 } |
|
|
1756 } |
|
|
1757 return retval; |
|
|
1758 } |
|
|
1759 |
|
|
1760 Matrix |
|
|
1761 Matrix::any (void) const |
|
|
1762 { |
|
|
1763 Matrix retval; |
|
|
1764 if (nr > 0 && nc > 0) |
|
|
1765 { |
|
|
1766 if (nr == 1) |
|
|
1767 { |
|
|
1768 retval.resize (1, 1); |
|
|
1769 retval.elem (0, 0) = 0.0; |
|
|
1770 for (int j = 0; j < nc; j++) |
|
|
1771 { |
|
|
1772 if (elem (0, j) != 0.0) |
|
|
1773 { |
|
|
1774 retval.elem (0, 0) = 1.0; |
|
|
1775 break; |
|
|
1776 } |
|
|
1777 } |
|
|
1778 } |
|
|
1779 else if (nc == 1) |
|
|
1780 { |
|
|
1781 retval.resize (1, 1); |
|
|
1782 retval.elem (0, 0) = 0.0; |
|
|
1783 for (int i = 0; i < nr; i++) |
|
|
1784 { |
|
|
1785 if (elem (i, 0) != 0.0) |
|
|
1786 { |
|
|
1787 retval.elem (0, 0) = 1.0; |
|
|
1788 break; |
|
|
1789 } |
|
|
1790 } |
|
|
1791 } |
|
|
1792 else |
|
|
1793 { |
|
|
1794 retval.resize (1, nc); |
|
|
1795 for (int j = 0; j < nc; j++) |
|
|
1796 { |
|
|
1797 retval.elem (0, j) = 0.0; |
|
|
1798 for (int i = 0; i < nr; i++) |
|
|
1799 { |
|
|
1800 if (elem (i, j) != 0.0) |
|
|
1801 { |
|
|
1802 retval.elem (0, j) = 1.0; |
|
|
1803 break; |
|
|
1804 } |
|
|
1805 } |
|
|
1806 } |
|
|
1807 } |
|
|
1808 } |
|
|
1809 return retval; |
|
|
1810 } |
|
|
1811 |
|
|
1812 Matrix |
|
|
1813 Matrix::cumprod (void) const |
|
|
1814 { |
|
|
1815 Matrix retval; |
|
|
1816 if (nr == 1) |
|
|
1817 { |
|
|
1818 retval.resize (1, nc); |
|
|
1819 if (nc > 0) |
|
|
1820 { |
|
|
1821 double prod = elem (0, 0); |
|
|
1822 for (int j = 0; j < nc; j++) |
|
|
1823 { |
|
|
1824 retval.elem (0, j) = prod; |
|
|
1825 if (j < nc - 1) |
|
|
1826 prod *= elem (0, j+1); |
|
|
1827 } |
|
|
1828 } |
|
|
1829 } |
|
|
1830 else if (nc == 1) |
|
|
1831 { |
|
|
1832 retval.resize (nr, 1); |
|
|
1833 if (nr > 0) |
|
|
1834 { |
|
|
1835 double prod = elem (0, 0); |
|
|
1836 for (int i = 0; i < nr; i++) |
|
|
1837 { |
|
|
1838 retval.elem (i, 0) = prod; |
|
|
1839 if (i < nr - 1) |
|
|
1840 prod *= elem (i+1, 0); |
|
|
1841 } |
|
|
1842 } |
|
|
1843 } |
|
|
1844 else |
|
|
1845 { |
|
|
1846 retval.resize (nr, nc); |
|
|
1847 if (nr > 0 && nc > 0) |
|
|
1848 { |
|
|
1849 for (int j = 0; j < nc; j++) |
|
|
1850 { |
|
|
1851 double prod = elem (0, j); |
|
|
1852 for (int i = 0; i < nr; i++) |
|
|
1853 { |
|
|
1854 retval.elem (i, j) = prod; |
|
|
1855 if (i < nr - 1) |
|
|
1856 prod *= elem (i+1, j); |
|
|
1857 } |
|
|
1858 } |
|
|
1859 } |
|
|
1860 } |
|
|
1861 return retval; |
|
|
1862 } |
|
|
1863 |
|
|
1864 Matrix |
|
|
1865 Matrix::cumsum (void) const |
|
|
1866 { |
|
|
1867 Matrix retval; |
|
|
1868 if (nr == 1) |
|
|
1869 { |
|
|
1870 retval.resize (1, nc); |
|
|
1871 if (nc > 0) |
|
|
1872 { |
|
|
1873 double sum = elem (0, 0); |
|
|
1874 for (int j = 0; j < nc; j++) |
|
|
1875 { |
|
|
1876 retval.elem (0, j) = sum; |
|
|
1877 if (j < nc - 1) |
|
|
1878 sum += elem (0, j+1); |
|
|
1879 } |
|
|
1880 } |
|
|
1881 } |
|
|
1882 else if (nc == 1) |
|
|
1883 { |
|
|
1884 retval.resize (nr, 1); |
|
|
1885 if (nr > 0) |
|
|
1886 { |
|
|
1887 double sum = elem (0, 0); |
|
|
1888 for (int i = 0; i < nr; i++) |
|
|
1889 { |
|
|
1890 retval.elem (i, 0) = sum; |
|
|
1891 if (i < nr - 1) |
|
|
1892 sum += elem (i+1, 0); |
|
|
1893 } |
|
|
1894 } |
|
|
1895 } |
|
|
1896 else |
|
|
1897 { |
|
|
1898 retval.resize (nr, nc); |
|
|
1899 if (nr > 0 && nc > 0) |
|
|
1900 { |
|
|
1901 for (int j = 0; j < nc; j++) |
|
|
1902 { |
|
|
1903 double sum = elem (0, j); |
|
|
1904 for (int i = 0; i < nr; i++) |
|
|
1905 { |
|
|
1906 retval.elem (i, j) = sum; |
|
|
1907 if (i < nr - 1) |
|
|
1908 sum += elem (i+1, j); |
|
|
1909 } |
|
|
1910 } |
|
|
1911 } |
|
|
1912 } |
|
|
1913 return retval; |
|
|
1914 } |
|
|
1915 |
|
|
1916 Matrix |
|
|
1917 Matrix::prod (void) const |
|
|
1918 { |
|
|
1919 Matrix retval; |
|
|
1920 if (nr == 1) |
|
|
1921 { |
|
|
1922 retval.resize (1, 1); |
|
|
1923 retval.elem (0, 0) = 1.0; |
|
|
1924 for (int j = 0; j < nc; j++) |
|
|
1925 retval.elem (0, 0) *= elem (0, j); |
|
|
1926 } |
|
|
1927 else if (nc == 1) |
|
|
1928 { |
|
|
1929 retval.resize (1, 1); |
|
|
1930 retval.elem (0, 0) = 1.0; |
|
|
1931 for (int i = 0; i < nr; i++) |
|
|
1932 retval.elem (0, 0) *= elem (i, 0); |
|
|
1933 } |
|
|
1934 else |
|
|
1935 { |
|
|
1936 if (nc == 0) |
|
|
1937 { |
|
|
1938 retval.resize (1, 1); |
|
|
1939 retval.elem (0, 0) = 1.0; |
|
|
1940 } |
|
|
1941 else |
|
|
1942 retval.resize (1, nc); |
|
|
1943 |
|
|
1944 for (int j = 0; j < nc; j++) |
|
|
1945 { |
|
|
1946 retval.elem (0, j) = 1.0; |
|
|
1947 for (int i = 0; i < nr; i++) |
|
|
1948 retval.elem (0, j) *= elem (i, j); |
|
|
1949 } |
|
|
1950 } |
|
|
1951 return retval; |
|
|
1952 } |
|
|
1953 |
|
|
1954 Matrix |
|
|
1955 Matrix::sum (void) const |
|
|
1956 { |
|
|
1957 Matrix retval; |
|
|
1958 if (nr == 1) |
|
|
1959 { |
|
|
1960 retval.resize (1, 1); |
|
|
1961 retval.elem (0, 0) = 0.0; |
|
|
1962 for (int j = 0; j < nc; j++) |
|
|
1963 retval.elem (0, 0) += elem (0, j); |
|
|
1964 } |
|
|
1965 else if (nc == 1) |
|
|
1966 { |
|
|
1967 retval.resize (1, 1); |
|
|
1968 retval.elem (0, 0) = 0.0; |
|
|
1969 for (int i = 0; i < nr; i++) |
|
|
1970 retval.elem (0, 0) += elem (i, 0); |
|
|
1971 } |
|
|
1972 else |
|
|
1973 { |
|
|
1974 if (nc == 0) |
|
|
1975 { |
|
|
1976 retval.resize (1, 1); |
|
|
1977 retval.elem (0, 0) = 0.0; |
|
|
1978 } |
|
|
1979 else |
|
|
1980 retval.resize (1, nc); |
|
|
1981 |
|
|
1982 for (int j = 0; j < nc; j++) |
|
|
1983 { |
|
|
1984 retval.elem (0, j) = 0.0; |
|
|
1985 for (int i = 0; i < nr; i++) |
|
|
1986 retval.elem (0, j) += elem (i, j); |
|
|
1987 } |
|
|
1988 } |
|
|
1989 return retval; |
|
|
1990 } |
|
|
1991 |
|
|
1992 Matrix |
|
|
1993 Matrix::sumsq (void) const |
|
|
1994 { |
|
|
1995 Matrix retval; |
|
|
1996 if (nr == 1) |
|
|
1997 { |
|
|
1998 retval.resize (1, 1); |
|
|
1999 retval.elem (0, 0) = 0.0; |
|
|
2000 for (int j = 0; j < nc; j++) |
|
|
2001 { |
|
|
2002 double d = elem (0, j); |
|
|
2003 retval.elem (0, 0) += d * d; |
|
|
2004 } |
|
|
2005 } |
|
|
2006 else if (nc == 1) |
|
|
2007 { |
|
|
2008 retval.resize (1, 1); |
|
|
2009 retval.elem (0, 0) = 0.0; |
|
|
2010 for (int i = 0; i < nr; i++) |
|
|
2011 { |
|
|
2012 double d = elem (i, 0); |
|
|
2013 retval.elem (0, 0) += d * d; |
|
|
2014 } |
|
|
2015 } |
|
|
2016 else |
|
|
2017 { |
|
|
2018 retval.resize (1, nc); |
|
|
2019 for (int j = 0; j < nc; j++) |
|
|
2020 { |
|
|
2021 retval.elem (0, j) = 0.0; |
|
|
2022 for (int i = 0; i < nr; i++) |
|
|
2023 { |
|
|
2024 double d = elem (i, j); |
|
|
2025 retval.elem (0, j) += d * d; |
|
|
2026 } |
|
|
2027 } |
|
|
2028 } |
|
|
2029 return retval; |
|
|
2030 } |
|
|
2031 |
|
|
2032 ColumnVector |
|
|
2033 Matrix::diag (void) const |
|
|
2034 { |
|
|
2035 return diag (0); |
|
|
2036 } |
|
|
2037 |
|
|
2038 ColumnVector |
|
|
2039 Matrix::diag (int k) const |
|
|
2040 { |
|
|
2041 int nnr = nr; |
|
|
2042 int nnc = nc; |
|
|
2043 if (k > 0) |
|
|
2044 nnc -= k; |
|
|
2045 else if (k < 0) |
|
|
2046 nnr += k; |
|
|
2047 |
|
|
2048 ColumnVector d; |
|
|
2049 |
|
|
2050 if (nnr > 0 && nnc > 0) |
|
|
2051 { |
|
|
2052 int ndiag = (nnr < nnc) ? nnr : nnc; |
|
|
2053 |
|
|
2054 d.resize (ndiag); |
|
|
2055 |
|
|
2056 if (k > 0) |
|
|
2057 { |
|
|
2058 for (int i = 0; i < ndiag; i++) |
|
|
2059 d.elem (i) = elem (i, i+k); |
|
|
2060 } |
|
|
2061 else if ( k < 0) |
|
|
2062 { |
|
|
2063 for (int i = 0; i < ndiag; i++) |
|
|
2064 d.elem (i) = elem (i-k, i); |
|
|
2065 } |
|
|
2066 else |
|
|
2067 { |
|
|
2068 for (int i = 0; i < ndiag; i++) |
|
|
2069 d.elem (i) = elem (i, i); |
|
|
2070 } |
|
|
2071 } |
|
|
2072 else |
|
|
2073 cerr << "diag: requested diagonal out of range\n"; |
|
|
2074 |
|
|
2075 return d; |
|
|
2076 } |
|
|
2077 |
|
|
2078 // unary operations |
|
|
2079 |
|
|
2080 Matrix |
|
|
2081 Matrix::operator - (void) const |
|
|
2082 { |
|
|
2083 return Matrix (negate (data, len), nr, nc); |
|
|
2084 } |
|
|
2085 |
|
|
2086 Matrix |
|
|
2087 Matrix::operator ! (void) const |
|
|
2088 { |
|
|
2089 Matrix b (nr, nc); |
|
|
2090 |
|
|
2091 for (int j = 0; j < nc; j++) |
|
|
2092 for (int i = 0; i < nr; i++) |
|
|
2093 b.elem (i, j) = ! elem (i, j); |
|
|
2094 |
|
|
2095 return b; |
|
|
2096 } |
|
|
2097 |
|
|
2098 ColumnVector |
|
|
2099 Matrix::row_min (void) const |
|
|
2100 { |
|
|
2101 ColumnVector result; |
|
|
2102 |
|
|
2103 if (nr > 0 && nc > 0) |
|
|
2104 { |
|
|
2105 result.resize (nr); |
|
|
2106 |
|
|
2107 for (int i = 0; i < nr; i++) |
|
|
2108 { |
|
|
2109 double res = elem (i, 0); |
|
|
2110 for (int j = 1; j < nc; j++) |
|
|
2111 if (elem (i, j) < res) |
|
|
2112 res = elem (i, j); |
|
|
2113 result.elem (i) = res; |
|
|
2114 } |
|
|
2115 } |
|
|
2116 |
|
|
2117 return result; |
|
|
2118 } |
|
|
2119 |
|
|
2120 ColumnVector |
|
210
|
2121 Matrix::row_min_loc (void) const |
|
|
2122 { |
|
|
2123 ColumnVector result; |
|
|
2124 |
|
|
2125 if (nr > 0 && nc > 0) |
|
|
2126 { |
|
|
2127 result.resize (nr); |
|
|
2128 |
|
|
2129 for (int i = 0; i < nr; i++) |
|
|
2130 { |
|
|
2131 int res = 0; |
|
|
2132 for (int j = 0; j < nc; j++) |
|
|
2133 if (elem (i, j) < elem (i, res)) |
|
|
2134 res = j; |
|
|
2135 result.elem (i) = (double) (res + 1); |
|
|
2136 } |
|
|
2137 } |
|
|
2138 |
|
|
2139 return result; |
|
|
2140 } |
|
|
2141 |
|
|
2142 ColumnVector |
|
3
|
2143 Matrix::row_max (void) const |
|
|
2144 { |
|
|
2145 ColumnVector result; |
|
|
2146 |
|
|
2147 if (nr > 0 && nc > 0) |
|
|
2148 { |
|
|
2149 result.resize (nr); |
|
|
2150 |
|
|
2151 for (int i = 0; i < nr; i++) |
|
|
2152 { |
|
|
2153 double res = elem (i, 0); |
|
|
2154 for (int j = 1; j < nc; j++) |
|
|
2155 if (elem (i, j) > res) |
|
|
2156 res = elem (i, j); |
|
|
2157 result.elem (i) = res; |
|
|
2158 } |
|
|
2159 } |
|
|
2160 |
|
|
2161 return result; |
|
|
2162 } |
|
|
2163 |
|
210
|
2164 ColumnVector |
|
|
2165 Matrix::row_max_loc (void) const |
|
|
2166 { |
|
|
2167 ColumnVector result; |
|
|
2168 |
|
|
2169 if (nr > 0 && nc > 0) |
|
|
2170 { |
|
|
2171 result.resize (nr); |
|
|
2172 |
|
|
2173 for (int i = 0; i < nr; i++) |
|
|
2174 { |
|
|
2175 int res = 0; |
|
|
2176 for (int j = 0; j < nc; j++) |
|
|
2177 if (elem (i, j) > elem (i, res)) |
|
|
2178 res = j; |
|
|
2179 result.elem (i) = (double) (res + 1); |
|
|
2180 } |
|
|
2181 } |
|
|
2182 |
|
|
2183 return result; |
|
|
2184 } |
|
|
2185 |
|
3
|
2186 RowVector |
|
|
2187 Matrix::column_min (void) const |
|
|
2188 { |
|
|
2189 RowVector result; |
|
|
2190 |
|
|
2191 if (nr > 0 && nc > 0) |
|
|
2192 { |
|
|
2193 result.resize (nc); |
|
|
2194 |
|
|
2195 for (int j = 0; j < nc; j++) |
|
|
2196 { |
|
|
2197 double res = elem (0, j); |
|
|
2198 for (int i = 1; i < nr; i++) |
|
|
2199 if (elem (i, j) < res) |
|
|
2200 res = elem (i, j); |
|
|
2201 result.elem (j) = res; |
|
|
2202 } |
|
|
2203 } |
|
|
2204 |
|
|
2205 return result; |
|
|
2206 } |
|
210
|
2207 RowVector |
|
|
2208 Matrix::column_min_loc (void) const |
|
|
2209 { |
|
|
2210 RowVector result; |
|
|
2211 |
|
|
2212 if (nr > 0 && nc > 0) |
|
|
2213 { |
|
|
2214 result.resize (nc); |
|
|
2215 |
|
|
2216 for (int j = 0; j < nc; j++) |
|
|
2217 { |
|
227
|
2218 int res = 0; |
|
210
|
2219 for (int i = 0; i < nr; i++) |
|
|
2220 if (elem (i, j) < elem (res, j)) |
|
|
2221 res = i; |
|
|
2222 result.elem (j) = (double) (res + 1); |
|
|
2223 } |
|
|
2224 } |
|
|
2225 |
|
|
2226 return result; |
|
|
2227 } |
|
|
2228 |
|
3
|
2229 |
|
|
2230 RowVector |
|
|
2231 Matrix::column_max (void) const |
|
|
2232 { |
|
|
2233 RowVector result; |
|
|
2234 |
|
|
2235 if (nr > 0 && nc > 0) |
|
|
2236 { |
|
|
2237 result.resize (nc); |
|
|
2238 |
|
|
2239 for (int j = 0; j < nc; j++) |
|
|
2240 { |
|
|
2241 double res = elem (0, j); |
|
|
2242 for (int i = 1; i < nr; i++) |
|
|
2243 if (elem (i, j) > res) |
|
|
2244 res = elem (i, j); |
|
|
2245 result.elem (j) = res; |
|
|
2246 } |
|
|
2247 } |
|
|
2248 |
|
|
2249 return result; |
|
|
2250 } |
|
|
2251 |
|
210
|
2252 RowVector |
|
|
2253 Matrix::column_max_loc (void) const |
|
|
2254 { |
|
|
2255 RowVector result; |
|
|
2256 |
|
|
2257 if (nr > 0 && nc > 0) |
|
|
2258 { |
|
|
2259 result.resize (nc); |
|
|
2260 |
|
|
2261 for (int j = 0; j < nc; j++) |
|
|
2262 { |
|
|
2263 int res = 0; |
|
|
2264 for (int i = 0; i < nr; i++) |
|
|
2265 if (elem (i, j) > elem (res, j)) |
|
|
2266 res = i; |
|
|
2267 result.elem (j) = (double) (res + 1); |
|
|
2268 } |
|
|
2269 } |
|
|
2270 |
|
|
2271 return result; |
|
|
2272 } |
|
|
2273 |
|
3
|
2274 ostream& |
|
|
2275 operator << (ostream& os, const Matrix& a) |
|
|
2276 { |
|
|
2277 // int field_width = os.precision () + 7; |
|
|
2278 for (int i = 0; i < a.nr; i++) |
|
|
2279 { |
|
|
2280 for (int j = 0; j < a.nc; j++) |
|
|
2281 os << " " /* setw (field_width) */ << a.elem (i, j); |
|
|
2282 os << "\n"; |
|
|
2283 } |
|
|
2284 return os; |
|
|
2285 } |
|
|
2286 |
|
|
2287 istream& |
|
|
2288 operator >> (istream& is, Matrix& a) |
|
|
2289 { |
|
|
2290 int nr = a.rows (); |
|
|
2291 int nc = a.columns (); |
|
|
2292 |
|
|
2293 if (nr < 1 || nc < 1) |
|
|
2294 is.clear (ios::badbit); |
|
|
2295 else |
|
|
2296 { |
|
|
2297 double tmp; |
|
|
2298 for (int i = 0; i < nr; i++) |
|
|
2299 for (int j = 0; j < nc; j++) |
|
|
2300 { |
|
|
2301 is >> tmp; |
|
|
2302 if (is) |
|
|
2303 a.elem (i, j) = tmp; |
|
|
2304 else |
|
|
2305 break; |
|
|
2306 } |
|
|
2307 } |
|
|
2308 |
|
|
2309 return is; |
|
|
2310 } |
|
|
2311 |
|
|
2312 /* |
|
|
2313 * Complex Matrix class |
|
|
2314 */ |
|
|
2315 |
|
|
2316 ComplexMatrix::ComplexMatrix (int r, int c) |
|
|
2317 { |
|
|
2318 if (r < 0 || c < 0) |
|
227
|
2319 { |
|
|
2320 (*current_liboctave_error_handler) |
|
|
2321 ("can't construct matrix with negative dimensions"); |
|
|
2322 nr = 0; |
|
|
2323 nc = 0; |
|
|
2324 len = 0; |
|
|
2325 data = (Complex *) NULL; |
|
|
2326 return; |
|
|
2327 } |
|
3
|
2328 |
|
|
2329 nr = r; |
|
|
2330 nc = c; |
|
|
2331 len = nr * nc; |
|
|
2332 if (len > 0) |
|
|
2333 data = new Complex [len]; |
|
|
2334 else |
|
|
2335 data = (Complex *) NULL; |
|
|
2336 } |
|
|
2337 |
|
|
2338 ComplexMatrix::ComplexMatrix (int r, int c, double val) |
|
|
2339 { |
|
|
2340 if (r < 0 || c < 0) |
|
227
|
2341 { |
|
|
2342 (*current_liboctave_error_handler) |
|
|
2343 ("can't construct matrix with negative dimensions"); |
|
|
2344 nr = 0; |
|
|
2345 nc = 0; |
|
|
2346 len = 0; |
|
|
2347 data = (Complex *) NULL; |
|
|
2348 return; |
|
|
2349 } |
|
3
|
2350 |
|
|
2351 nr = r; |
|
|
2352 nc = c; |
|
|
2353 len = nr * nc; |
|
|
2354 if (len > 0) |
|
|
2355 { |
|
|
2356 data = new Complex [len]; |
|
|
2357 copy (data, len, val); |
|
|
2358 } |
|
|
2359 else |
|
|
2360 data = (Complex *) NULL; |
|
|
2361 } |
|
|
2362 |
|
161
|
2363 ComplexMatrix::ComplexMatrix (int r, int c, const Complex& val) |
|
3
|
2364 { |
|
|
2365 if (r < 0 || c < 0) |
|
227
|
2366 { |
|
|
2367 (*current_liboctave_error_handler) |
|
|
2368 ("can't construct matrix with negative dimensions"); |
|
|
2369 nr = 0; |
|
|
2370 nc = 0; |
|
|
2371 len = 0; |
|
|
2372 data = (Complex *) NULL; |
|
|
2373 return; |
|
|
2374 } |
|
3
|
2375 |
|
|
2376 nr = r; |
|
|
2377 nc = c; |
|
|
2378 len = nr * nc; |
|
|
2379 if (len > 0) |
|
|
2380 { |
|
|
2381 data = new Complex [len]; |
|
|
2382 copy (data, len, val); |
|
|
2383 } |
|
|
2384 else |
|
|
2385 data = (Complex *) NULL; |
|
|
2386 } |
|
|
2387 |
|
|
2388 ComplexMatrix::ComplexMatrix (const Matrix& a) |
|
|
2389 { |
|
|
2390 nr = a.nr; |
|
|
2391 nc = a.nc; |
|
|
2392 len = a.len; |
|
|
2393 if (len > 0) |
|
|
2394 { |
|
|
2395 data = new Complex [len]; |
|
|
2396 copy (data, a.data, len); |
|
|
2397 } |
|
|
2398 else |
|
|
2399 data = (Complex *) NULL; |
|
|
2400 } |
|
|
2401 |
|
|
2402 ComplexMatrix::ComplexMatrix (const ComplexMatrix& a) |
|
|
2403 { |
|
|
2404 nr = a.nr; |
|
|
2405 nc = a.nc; |
|
|
2406 len = a.len; |
|
|
2407 if (len > 0) |
|
|
2408 { |
|
|
2409 data = new Complex [len]; |
|
|
2410 copy (data, a.data, len); |
|
|
2411 } |
|
|
2412 else |
|
|
2413 data = (Complex *) NULL; |
|
|
2414 } |
|
|
2415 |
|
|
2416 ComplexMatrix::ComplexMatrix (const DiagMatrix& a) |
|
|
2417 { |
|
|
2418 nr = a.nr; |
|
|
2419 nc = a.nc; |
|
|
2420 len = nr * nc; |
|
|
2421 if (len > 0) |
|
|
2422 { |
|
|
2423 data = new Complex [len]; |
|
|
2424 copy (data, len, 0.0); |
|
|
2425 for (int i = 0; i < a.len; i++) |
|
|
2426 data[nr*i+i] = a.data[i]; |
|
|
2427 } |
|
|
2428 else |
|
|
2429 data = (Complex *) NULL; |
|
|
2430 } |
|
|
2431 |
|
|
2432 ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a) |
|
|
2433 { |
|
|
2434 nr = a.nr; |
|
|
2435 nc = a.nc; |
|
|
2436 len = nr * nc; |
|
|
2437 if (len > 0) |
|
|
2438 { |
|
|
2439 data = new Complex [len]; |
|
|
2440 copy (data, len, 0.0); |
|
|
2441 for (int i = 0; i < a.len; i++) |
|
|
2442 data[nr*i+i] = a.data[i]; |
|
|
2443 } |
|
|
2444 else |
|
|
2445 data = (Complex *) NULL; |
|
|
2446 } |
|
|
2447 |
|
|
2448 ComplexMatrix::ComplexMatrix (double a) |
|
|
2449 { |
|
|
2450 nr = 1; |
|
|
2451 nc = 1; |
|
|
2452 len = 1; |
|
|
2453 data = new Complex [1]; |
|
|
2454 data[0] = a; |
|
|
2455 } |
|
|
2456 |
|
161
|
2457 ComplexMatrix::ComplexMatrix (const Complex& a) |
|
3
|
2458 { |
|
|
2459 nr = 1; |
|
|
2460 nc = 1; |
|
|
2461 len = 1; |
|
|
2462 data = new Complex [1]; |
|
|
2463 data[0] = Complex (a); |
|
|
2464 } |
|
|
2465 |
|
|
2466 ComplexMatrix& |
|
|
2467 ComplexMatrix::operator = (const Matrix& a) |
|
|
2468 { |
|
|
2469 delete [] data; |
|
|
2470 nr = a.nr; |
|
|
2471 nc = a.nc; |
|
|
2472 len = a.len; |
|
|
2473 if (len > 0) |
|
|
2474 { |
|
|
2475 data = new Complex [len]; |
|
|
2476 copy (data, a.data, len); |
|
|
2477 } |
|
|
2478 else |
|
|
2479 data = (Complex *) NULL; |
|
|
2480 return *this; |
|
|
2481 } |
|
|
2482 |
|
|
2483 ComplexMatrix& |
|
|
2484 ComplexMatrix::operator = (const ComplexMatrix& a) |
|
|
2485 { |
|
|
2486 if (this != &a) |
|
|
2487 { |
|
|
2488 delete [] data; |
|
|
2489 nr = a.nr; |
|
|
2490 nc = a.nc; |
|
|
2491 len = a.len; |
|
|
2492 if (len > 0) |
|
|
2493 { |
|
|
2494 data = new Complex [len]; |
|
|
2495 copy (data, a.data, len); |
|
|
2496 } |
|
|
2497 else |
|
|
2498 data = (Complex *) NULL; |
|
|
2499 } |
|
|
2500 return *this; |
|
|
2501 } |
|
|
2502 |
|
227
|
2503 Complex& |
|
|
2504 ComplexMatrix::checkelem (int r, int c) |
|
|
2505 { |
|
|
2506 #ifndef NO_RANGE_CHECK |
|
|
2507 if (r < 0 || r >= nr || c < 0 || c >= nc) |
|
|
2508 { |
|
|
2509 (*current_liboctave_error_handler) ("range error"); |
|
|
2510 static Complex foo (0.0); |
|
|
2511 return foo; |
|
|
2512 } |
|
|
2513 #endif |
|
|
2514 |
|
|
2515 return elem (r, c); |
|
|
2516 } |
|
|
2517 |
|
|
2518 Complex |
|
|
2519 ComplexMatrix::checkelem (int r, int c) const |
|
|
2520 { |
|
|
2521 #ifndef NO_RANGE_CHECK |
|
|
2522 if (r < 0 || r >= nr || c < 0 || c >= nc) |
|
|
2523 { |
|
|
2524 (*current_liboctave_error_handler) ("range error"); |
|
|
2525 return Complex (0.0); |
|
|
2526 } |
|
|
2527 #endif |
|
|
2528 |
|
|
2529 return elem (r, c); |
|
|
2530 } |
|
|
2531 |
|
3
|
2532 ComplexMatrix& |
|
|
2533 ComplexMatrix::resize (int r, int c) |
|
|
2534 { |
|
|
2535 if (r < 0 || c < 0) |
|
227
|
2536 { |
|
|
2537 (*current_liboctave_error_handler) |
|
|
2538 ("can't resize to negative dimensions"); |
|
|
2539 return *this; |
|
|
2540 } |
|
3
|
2541 |
|
|
2542 int new_len = r * c; |
|
|
2543 Complex* new_data = (Complex *) NULL; |
|
|
2544 if (new_len > 0) |
|
|
2545 { |
|
|
2546 new_data = new Complex [new_len]; |
|
|
2547 |
|
|
2548 int min_r = nr < r ? nr : r; |
|
|
2549 int min_c = nc < c ? nc : c; |
|
|
2550 |
|
|
2551 for (int j = 0; j < min_c; j++) |
|
|
2552 for (int i = 0; i < min_r; i++) |
|
|
2553 new_data[r*j+i] = elem (i, j); |
|
|
2554 } |
|
|
2555 |
|
|
2556 delete [] data; |
|
|
2557 nr = r; |
|
|
2558 nc = c; |
|
|
2559 len = new_len; |
|
|
2560 data = new_data; |
|
|
2561 |
|
|
2562 return *this; |
|
|
2563 } |
|
|
2564 |
|
|
2565 ComplexMatrix& |
|
|
2566 ComplexMatrix::resize (int r, int c, double val) |
|
|
2567 { |
|
|
2568 if (r < 0 || c < 0) |
|
227
|
2569 { |
|
|
2570 (*current_liboctave_error_handler) |
|
|
2571 ("can't resize to negative dimensions"); |
|
|
2572 return *this; |
|
|
2573 } |
|
3
|
2574 |
|
|
2575 int new_len = r * c; |
|
|
2576 Complex *new_data = (Complex *) NULL; |
|
|
2577 if (new_len > 0) |
|
|
2578 { |
|
|
2579 new_data = new Complex [new_len]; |
|
|
2580 |
|
|
2581 // There may be faster or cleaner ways to do this. |
|
|
2582 |
|
|
2583 if (r > nr || c > nc) |
|
|
2584 copy (new_data, new_len, val); |
|
|
2585 |
|
|
2586 int min_r = nr < r ? nr : r; |
|
|
2587 int min_c = nc < c ? nc : c; |
|
|
2588 |
|
|
2589 for (int j = 0; j < min_c; j++) |
|
|
2590 for (int i = 0; i < min_r; i++) |
|
|
2591 new_data[r*j+i] = elem (i, j); |
|
|
2592 } |
|
|
2593 |
|
|
2594 delete [] data; |
|
|
2595 nr = r; |
|
|
2596 nc = c; |
|
|
2597 len = new_len; |
|
|
2598 data = new_data; |
|
|
2599 |
|
|
2600 return *this; |
|
|
2601 } |
|
|
2602 |
|
|
2603 ComplexMatrix& |
|
161
|
2604 ComplexMatrix::resize (int r, int c, const Complex& val) |
|
3
|
2605 { |
|
|
2606 if (r < 0 || c < 0) |
|
227
|
2607 { |
|
|
2608 (*current_liboctave_error_handler) |
|
|
2609 ("can't resize to negative dimensions"); |
|
|
2610 return *this; |
|
|
2611 } |
|
3
|
2612 |
|
|
2613 int new_len = r * c; |
|
|
2614 Complex *new_data = (Complex *) NULL; |
|
|
2615 if (new_len > 0) |
|
|
2616 { |
|
|
2617 new_data = new Complex [new_len]; |
|
|
2618 |
|
|
2619 // There may be faster or cleaner ways to do this. |
|
|
2620 |
|
|
2621 if (r > nr || c > nc) |
|
|
2622 copy (new_data, new_len, val); |
|
|
2623 |
|
|
2624 int min_r = nr < r ? nr : r; |
|
|
2625 int min_c = nc < c ? nc : c; |
|
|
2626 |
|
|
2627 for (int j = 0; j < min_c; j++) |
|
|
2628 for (int i = 0; i < min_r; i++) |
|
|
2629 new_data[r*j+i] = elem (i, j); |
|
|
2630 } |
|
|
2631 |
|
|
2632 delete [] data; |
|
|
2633 nr = r; |
|
|
2634 nc = c; |
|
|
2635 len = new_len; |
|
|
2636 data = new_data; |
|
|
2637 |
|
|
2638 return *this; |
|
|
2639 } |
|
|
2640 |
|
|
2641 int |
|
|
2642 ComplexMatrix::operator == (const ComplexMatrix& a) const |
|
|
2643 { |
|
|
2644 if (nr != a.nr || nc != a.nc) |
|
|
2645 return 0; |
|
|
2646 |
|
|
2647 return equal (data, a.data, len); |
|
|
2648 } |
|
|
2649 |
|
|
2650 int |
|
|
2651 ComplexMatrix::operator != (const ComplexMatrix& a) const |
|
|
2652 { |
|
|
2653 return !(*this == a); |
|
|
2654 } |
|
|
2655 |
|
|
2656 // destructive insert/delete/reorder operations |
|
|
2657 |
|
|
2658 ComplexMatrix& |
|
|
2659 ComplexMatrix::insert (const Matrix& a, int r, int c) |
|
|
2660 { |
|
|
2661 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
2662 { |
|
|
2663 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2664 return *this; |
|
|
2665 } |
|
3
|
2666 |
|
|
2667 for (int j = 0; j < a.nc; j++) |
|
|
2668 for (int i = 0; i < a.nr; i++) |
|
|
2669 elem (r+i, c+j) = a.elem (i, j); |
|
|
2670 |
|
|
2671 return *this; |
|
|
2672 } |
|
|
2673 |
|
|
2674 ComplexMatrix& |
|
|
2675 ComplexMatrix::insert (const RowVector& a, int r, int c) |
|
|
2676 { |
|
|
2677 if (r < 0 || r >= nr || c < 0 || c + a.len - 1 > nc) |
|
227
|
2678 { |
|
|
2679 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2680 return *this; |
|
|
2681 } |
|
3
|
2682 |
|
|
2683 for (int i = 0; i < a.len; i++) |
|
|
2684 elem (r, c+i) = a.data[i]; |
|
|
2685 |
|
|
2686 return *this; |
|
|
2687 } |
|
|
2688 |
|
|
2689 ComplexMatrix& |
|
|
2690 ComplexMatrix::insert (const ColumnVector& a, int r, int c) |
|
|
2691 { |
|
|
2692 if (r < 0 || r + a.len - 1 > nr || c < 0 || c >= nc) |
|
227
|
2693 { |
|
|
2694 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2695 return *this; |
|
|
2696 } |
|
3
|
2697 |
|
|
2698 for (int i = 0; i < a.len; i++) |
|
|
2699 elem (r+i, c) = a.data[i]; |
|
|
2700 |
|
|
2701 return *this; |
|
|
2702 } |
|
|
2703 |
|
|
2704 ComplexMatrix& |
|
|
2705 ComplexMatrix::insert (const DiagMatrix& a, int r, int c) |
|
|
2706 { |
|
|
2707 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
2708 { |
|
|
2709 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2710 return *this; |
|
|
2711 } |
|
3
|
2712 |
|
|
2713 for (int i = 0; i < a.len; i++) |
|
|
2714 elem (r+i, c+i) = a.data[i]; |
|
|
2715 |
|
|
2716 return *this; |
|
|
2717 } |
|
|
2718 |
|
|
2719 ComplexMatrix& |
|
|
2720 ComplexMatrix::insert (const ComplexMatrix& a, int r, int c) |
|
|
2721 { |
|
|
2722 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
2723 { |
|
|
2724 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2725 return *this; |
|
|
2726 } |
|
3
|
2727 |
|
|
2728 for (int j = 0; j < a.nc; j++) |
|
|
2729 for (int i = 0; i < a.nr; i++) |
|
|
2730 elem (r+i, c+j) = a.elem (i, j); |
|
|
2731 |
|
|
2732 return *this; |
|
|
2733 } |
|
|
2734 |
|
|
2735 ComplexMatrix& |
|
|
2736 ComplexMatrix::insert (const ComplexRowVector& a, int r, int c) |
|
|
2737 { |
|
|
2738 if (r < 0 || r >= nr || c < 0 || c + a.len - 1 > nc) |
|
227
|
2739 { |
|
|
2740 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2741 return *this; |
|
|
2742 } |
|
3
|
2743 |
|
|
2744 for (int i = 0; i < a.len; i++) |
|
|
2745 elem (r, c+i) = a.data[i]; |
|
|
2746 |
|
|
2747 return *this; |
|
|
2748 } |
|
|
2749 |
|
|
2750 ComplexMatrix& |
|
|
2751 ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c) |
|
|
2752 { |
|
|
2753 if (r < 0 || r + a.len - 1 > nr || c < 0 || c >= nc) |
|
227
|
2754 { |
|
|
2755 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2756 return *this; |
|
|
2757 } |
|
3
|
2758 |
|
|
2759 for (int i = 0; i < a.len; i++) |
|
|
2760 elem (r+i, c) = a.data[i]; |
|
|
2761 |
|
|
2762 return *this; |
|
|
2763 } |
|
|
2764 |
|
|
2765 ComplexMatrix& |
|
|
2766 ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c) |
|
|
2767 { |
|
|
2768 if (r < 0 || r + a.nr - 1 > nr || c < 0 || c + a.nc - 1 > nc) |
|
227
|
2769 { |
|
|
2770 (*current_liboctave_error_handler) ("range error for insert"); |
|
|
2771 return *this; |
|
|
2772 } |
|
3
|
2773 |
|
|
2774 for (int i = 0; i < a.len; i++) |
|
|
2775 elem (r+i, c+i) = a.data[i]; |
|
|
2776 |
|
|
2777 return *this; |
|
|
2778 } |
|
|
2779 |
|
|
2780 ComplexMatrix& |
|
|
2781 ComplexMatrix::fill (double val) |
|
|
2782 { |
|
|
2783 if (nr > 0 && nc > 0) |
|
|
2784 copy (data, len, val); |
|
|
2785 return *this; |
|
|
2786 } |
|
|
2787 |
|
|
2788 ComplexMatrix& |
|
161
|
2789 ComplexMatrix::fill (const Complex& val) |
|
3
|
2790 { |
|
|
2791 if (nr > 0 && nc > 0) |
|
|
2792 copy (data, len, val); |
|
|
2793 return *this; |
|
|
2794 } |
|
|
2795 |
|
|
2796 ComplexMatrix& |
|
|
2797 ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2) |
|
|
2798 { |
|
|
2799 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 |
|
|
2800 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) |
|
227
|
2801 { |
|
|
2802 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
2803 return *this; |
|
|
2804 } |
|
3
|
2805 |
|
|
2806 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
2807 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
2808 |
|
|
2809 for (int j = c1; j <= c2; j++) |
|
|
2810 for (int i = r1; i <= r2; i++) |
|
|
2811 elem (i, j) = val; |
|
|
2812 |
|
|
2813 return *this; |
|
|
2814 } |
|
|
2815 |
|
|
2816 ComplexMatrix& |
|
161
|
2817 ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2) |
|
3
|
2818 { |
|
|
2819 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 |
|
|
2820 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) |
|
227
|
2821 { |
|
|
2822 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
2823 return *this; |
|
|
2824 } |
|
3
|
2825 |
|
|
2826 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
2827 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
2828 |
|
|
2829 for (int j = c1; j <= c2; j++) |
|
|
2830 for (int i = r1; i <= r2; i++) |
|
|
2831 elem (i, j) = val; |
|
|
2832 |
|
|
2833 return *this; |
|
|
2834 } |
|
|
2835 |
|
|
2836 ComplexMatrix |
|
|
2837 ComplexMatrix::append (const Matrix& a) const |
|
|
2838 { |
|
|
2839 if (nr != a.nr) |
|
227
|
2840 { |
|
|
2841 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2842 return *this; |
|
|
2843 } |
|
3
|
2844 |
|
|
2845 int nc_insert = nc; |
|
|
2846 ComplexMatrix retval (nr, nc + a.nc); |
|
|
2847 retval.insert (*this, 0, 0); |
|
|
2848 retval.insert (a, 0, nc_insert); |
|
|
2849 return retval; |
|
|
2850 } |
|
|
2851 |
|
|
2852 ComplexMatrix |
|
|
2853 ComplexMatrix::append (const RowVector& a) const |
|
|
2854 { |
|
|
2855 if (nr != 1) |
|
227
|
2856 { |
|
|
2857 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2858 return *this; |
|
|
2859 } |
|
3
|
2860 |
|
|
2861 int nc_insert = nc; |
|
|
2862 ComplexMatrix retval (nr, nc + a.len); |
|
|
2863 retval.insert (*this, 0, 0); |
|
|
2864 retval.insert (a, 0, nc_insert); |
|
|
2865 return retval; |
|
|
2866 } |
|
|
2867 |
|
|
2868 ComplexMatrix |
|
|
2869 ComplexMatrix::append (const ColumnVector& a) const |
|
|
2870 { |
|
|
2871 if (nr != a.len) |
|
227
|
2872 { |
|
|
2873 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2874 return *this; |
|
|
2875 } |
|
3
|
2876 |
|
|
2877 int nc_insert = nc; |
|
|
2878 ComplexMatrix retval (nr, nc + 1); |
|
|
2879 retval.insert (*this, 0, 0); |
|
|
2880 retval.insert (a, 0, nc_insert); |
|
|
2881 return retval; |
|
|
2882 } |
|
|
2883 |
|
|
2884 ComplexMatrix |
|
|
2885 ComplexMatrix::append (const DiagMatrix& a) const |
|
|
2886 { |
|
|
2887 if (nr != a.nr) |
|
227
|
2888 { |
|
|
2889 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2890 return *this; |
|
|
2891 } |
|
3
|
2892 |
|
|
2893 int nc_insert = nc; |
|
|
2894 ComplexMatrix retval (nr, nc + a.nc); |
|
|
2895 retval.insert (*this, 0, 0); |
|
|
2896 retval.insert (a, 0, nc_insert); |
|
|
2897 return retval; |
|
|
2898 } |
|
|
2899 |
|
|
2900 ComplexMatrix |
|
|
2901 ComplexMatrix::append (const ComplexMatrix& a) const |
|
|
2902 { |
|
|
2903 if (nr != a.nr) |
|
227
|
2904 { |
|
|
2905 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2906 return *this; |
|
|
2907 } |
|
3
|
2908 |
|
|
2909 int nc_insert = nc; |
|
|
2910 ComplexMatrix retval (nr, nc + a.nc); |
|
|
2911 retval.insert (*this, 0, 0); |
|
|
2912 retval.insert (a, 0, nc_insert); |
|
|
2913 return retval; |
|
|
2914 } |
|
|
2915 |
|
|
2916 ComplexMatrix |
|
|
2917 ComplexMatrix::append (const ComplexRowVector& a) const |
|
|
2918 { |
|
|
2919 if (nr != 1) |
|
227
|
2920 { |
|
|
2921 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2922 return *this; |
|
|
2923 } |
|
3
|
2924 |
|
|
2925 int nc_insert = nc; |
|
|
2926 ComplexMatrix retval (nr, nc + a.len); |
|
|
2927 retval.insert (*this, 0, 0); |
|
|
2928 retval.insert (a, 0, nc_insert); |
|
|
2929 return retval; |
|
|
2930 } |
|
|
2931 |
|
|
2932 ComplexMatrix |
|
|
2933 ComplexMatrix::append (const ComplexColumnVector& a) const |
|
|
2934 { |
|
|
2935 if (nr != a.len) |
|
227
|
2936 { |
|
|
2937 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2938 return *this; |
|
|
2939 } |
|
3
|
2940 |
|
|
2941 int nc_insert = nc; |
|
|
2942 ComplexMatrix retval (nr, nc + 1); |
|
|
2943 retval.insert (*this, 0, 0); |
|
|
2944 retval.insert (a, 0, nc_insert); |
|
|
2945 return retval; |
|
|
2946 } |
|
|
2947 |
|
|
2948 ComplexMatrix |
|
|
2949 ComplexMatrix::append (const ComplexDiagMatrix& a) const |
|
|
2950 { |
|
|
2951 if (nr != a.nr) |
|
227
|
2952 { |
|
|
2953 (*current_liboctave_error_handler) ("row dimension mismatch for append"); |
|
|
2954 return *this; |
|
|
2955 } |
|
3
|
2956 |
|
|
2957 int nc_insert = nc; |
|
|
2958 ComplexMatrix retval (nr, nc + a.nc); |
|
|
2959 retval.insert (*this, 0, 0); |
|
|
2960 retval.insert (a, 0, nc_insert); |
|
|
2961 return retval; |
|
|
2962 } |
|
|
2963 |
|
|
2964 ComplexMatrix |
|
|
2965 ComplexMatrix::stack (const Matrix& a) const |
|
|
2966 { |
|
|
2967 if (nc != a.nc) |
|
227
|
2968 { |
|
|
2969 (*current_liboctave_error_handler) |
|
|
2970 ("column dimension mismatch for stack"); |
|
|
2971 return *this; |
|
|
2972 } |
|
3
|
2973 |
|
|
2974 int nr_insert = nr; |
|
|
2975 ComplexMatrix retval (nr + a.nr, nc); |
|
|
2976 retval.insert (*this, 0, 0); |
|
|
2977 retval.insert (a, nr_insert, 0); |
|
|
2978 return retval; |
|
|
2979 } |
|
|
2980 |
|
|
2981 ComplexMatrix |
|
|
2982 ComplexMatrix::stack (const RowVector& a) const |
|
|
2983 { |
|
|
2984 if (nc != a.len) |
|
227
|
2985 { |
|
|
2986 (*current_liboctave_error_handler) |
|
|
2987 ("column dimension mismatch for stack"); |
|
|
2988 return *this; |
|
|
2989 } |
|
3
|
2990 |
|
|
2991 int nr_insert = nr; |
|
|
2992 ComplexMatrix retval (nr + 1, nc); |
|
|
2993 retval.insert (*this, 0, 0); |
|
|
2994 retval.insert (a, nr_insert, 0); |
|
|
2995 return retval; |
|
|
2996 } |
|
|
2997 |
|
|
2998 ComplexMatrix |
|
|
2999 ComplexMatrix::stack (const ColumnVector& a) const |
|
|
3000 { |
|
|
3001 if (nc != 1) |
|
227
|
3002 { |
|
|
3003 (*current_liboctave_error_handler) |
|
|
3004 ("column dimension mismatch for stack"); |
|
|
3005 return *this; |
|
|
3006 } |
|
3
|
3007 |
|
|
3008 int nr_insert = nr; |
|
|
3009 ComplexMatrix retval (nr + a.len, nc); |
|
|
3010 retval.insert (*this, 0, 0); |
|
|
3011 retval.insert (a, nr_insert, 0); |
|
|
3012 return retval; |
|
|
3013 } |
|
|
3014 |
|
|
3015 ComplexMatrix |
|
|
3016 ComplexMatrix::stack (const DiagMatrix& a) const |
|
|
3017 { |
|
|
3018 if (nc != a.nc) |
|
227
|
3019 { |
|
|
3020 (*current_liboctave_error_handler) |
|
|
3021 ("column dimension mismatch for stack"); |
|
|
3022 return *this; |
|
|
3023 } |
|
3
|
3024 |
|
|
3025 int nr_insert = nr; |
|
|
3026 ComplexMatrix retval (nr + a.nr, nc); |
|
|
3027 retval.insert (*this, 0, 0); |
|
|
3028 retval.insert (a, nr_insert, 0); |
|
|
3029 return retval; |
|
|
3030 } |
|
|
3031 |
|
|
3032 ComplexMatrix |
|
|
3033 ComplexMatrix::stack (const ComplexMatrix& a) const |
|
|
3034 { |
|
|
3035 if (nc != a.nc) |
|
227
|
3036 { |
|
|
3037 (*current_liboctave_error_handler) |
|
|
3038 ("column dimension mismatch for stack"); |
|
|
3039 return *this; |
|
|
3040 } |
|
3
|
3041 |
|
|
3042 int nr_insert = nr; |
|
|
3043 ComplexMatrix retval (nr + a.nr, nc); |
|
|
3044 retval.insert (*this, 0, 0); |
|
|
3045 retval.insert (a, nr_insert, 0); |
|
|
3046 return retval; |
|
|
3047 } |
|
|
3048 |
|
|
3049 ComplexMatrix |
|
|
3050 ComplexMatrix::stack (const ComplexRowVector& a) const |
|
|
3051 { |
|
|
3052 if (nc != a.len) |
|
227
|
3053 { |
|
|
3054 (*current_liboctave_error_handler) |
|
|
3055 ("column dimension mismatch for stack"); |
|
|
3056 return *this; |
|
|
3057 } |
|
3
|
3058 |
|
|
3059 int nr_insert = nr; |
|
|
3060 ComplexMatrix retval (nr + 1, nc); |
|
|
3061 retval.insert (*this, 0, 0); |
|
|
3062 retval.insert (a, nr_insert, 0); |
|
|
3063 return retval; |
|
|
3064 } |
|
|
3065 |
|
|
3066 ComplexMatrix |
|
|
3067 ComplexMatrix::stack (const ComplexColumnVector& a) const |
|
|
3068 { |
|
|
3069 if (nc != 1) |
|
227
|
3070 { |
|
|
3071 (*current_liboctave_error_handler) |
|
|
3072 ("column dimension mismatch for stack"); |
|
|
3073 return *this; |
|
|
3074 } |
|
3
|
3075 |
|
|
3076 int nr_insert = nr; |
|
|
3077 ComplexMatrix retval (nr + a.len, nc); |
|
|
3078 retval.insert (*this, 0, 0); |
|
|
3079 retval.insert (a, nr_insert, 0); |
|
|
3080 return retval; |
|
|
3081 } |
|
|
3082 |
|
|
3083 ComplexMatrix |
|
|
3084 ComplexMatrix::stack (const ComplexDiagMatrix& a) const |
|
|
3085 { |
|
|
3086 if (nc != a.nc) |
|
227
|
3087 { |
|
|
3088 (*current_liboctave_error_handler) |
|
|
3089 ("column dimension mismatch for stack"); |
|
|
3090 return *this; |
|
|
3091 } |
|
3
|
3092 |
|
|
3093 int nr_insert = nr; |
|
|
3094 ComplexMatrix retval (nr + a.nr, nc); |
|
|
3095 retval.insert (*this, 0, 0); |
|
|
3096 retval.insert (a, nr_insert, 0); |
|
|
3097 return retval; |
|
|
3098 } |
|
|
3099 |
|
|
3100 ComplexMatrix |
|
|
3101 ComplexMatrix::hermitian (void) const |
|
|
3102 { |
|
|
3103 ComplexMatrix result; |
|
|
3104 if (len > 0) |
|
|
3105 { |
|
|
3106 result.resize (nc, nr); |
|
|
3107 for (int j = 0; j < nc; j++) |
|
|
3108 for (int i = 0; i < nr; i++) |
|
|
3109 result.data[nc*i+j] = conj (data[nr*j+i]); |
|
|
3110 } |
|
|
3111 return result; |
|
|
3112 } |
|
|
3113 |
|
|
3114 ComplexMatrix |
|
|
3115 ComplexMatrix::transpose (void) const |
|
|
3116 { |
|
88
|
3117 ComplexMatrix result (nc, nr); |
|
3
|
3118 if (len > 0) |
|
|
3119 { |
|
|
3120 for (int j = 0; j < nc; j++) |
|
|
3121 for (int i = 0; i < nr; i++) |
|
|
3122 result.data[nc*i+j] = data[nr*j+i]; |
|
|
3123 } |
|
|
3124 return result; |
|
|
3125 } |
|
|
3126 |
|
|
3127 Matrix |
|
|
3128 real (const ComplexMatrix& a) |
|
|
3129 { |
|
|
3130 Matrix retval; |
|
|
3131 if (a.len > 0) |
|
|
3132 retval = Matrix (real_dup (a.data, a.len), a.nr, a.nc); |
|
|
3133 return retval; |
|
|
3134 } |
|
|
3135 |
|
|
3136 Matrix |
|
|
3137 imag (const ComplexMatrix& a) |
|
|
3138 { |
|
|
3139 Matrix retval; |
|
|
3140 if (a.len > 0) |
|
|
3141 retval = Matrix (imag_dup (a.data, a.len), a.nr, a.nc); |
|
|
3142 return retval; |
|
|
3143 } |
|
|
3144 |
|
|
3145 ComplexMatrix |
|
|
3146 conj (const ComplexMatrix& a) |
|
|
3147 { |
|
|
3148 ComplexMatrix retval; |
|
|
3149 if (a.len > 0) |
|
|
3150 retval = ComplexMatrix (conj_dup (a.data, a.len), a.nr, a.nc); |
|
|
3151 return retval; |
|
|
3152 } |
|
|
3153 |
|
|
3154 // resize is the destructive equivalent for this one |
|
|
3155 |
|
|
3156 ComplexMatrix |
|
|
3157 ComplexMatrix::extract (int r1, int c1, int r2, int c2) const |
|
|
3158 { |
|
|
3159 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
3160 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
3161 |
|
|
3162 int new_r = r2 - r1 + 1; |
|
|
3163 int new_c = c2 - c1 + 1; |
|
|
3164 |
|
|
3165 ComplexMatrix result (new_r, new_c); |
|
|
3166 |
|
|
3167 for (int j = 0; j < new_c; j++) |
|
|
3168 for (int i = 0; i < new_r; i++) |
|
|
3169 result.data[new_r*j+i] = elem (r1+i, c1+j); |
|
|
3170 |
|
|
3171 return result; |
|
|
3172 } |
|
|
3173 |
|
|
3174 // extract row or column i. |
|
|
3175 |
|
|
3176 ComplexRowVector |
|
|
3177 ComplexMatrix::row (int i) const |
|
|
3178 { |
|
|
3179 if (i < 0 || i >= nr) |
|
227
|
3180 { |
|
|
3181 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
3182 return ComplexRowVector (); |
|
|
3183 } |
|
3
|
3184 |
|
|
3185 ComplexRowVector retval (nc); |
|
|
3186 for (int j = 0; j < nc; j++) |
|
|
3187 retval.elem (j) = elem (i, j); |
|
|
3188 |
|
|
3189 return retval; |
|
|
3190 } |
|
|
3191 |
|
|
3192 ComplexRowVector |
|
|
3193 ComplexMatrix::row (char *s) const |
|
|
3194 { |
|
|
3195 if (s == (char *) NULL) |
|
227
|
3196 { |
|
|
3197 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
3198 return ComplexRowVector (); |
|
|
3199 } |
|
3
|
3200 |
|
|
3201 char c = *s; |
|
|
3202 if (c == 'f' || c == 'F') |
|
|
3203 return row (0); |
|
|
3204 else if (c == 'l' || c == 'L') |
|
|
3205 return row (nr - 1); |
|
|
3206 else |
|
227
|
3207 { |
|
|
3208 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
3209 return ComplexRowVector (); |
|
|
3210 } |
|
3
|
3211 } |
|
|
3212 |
|
|
3213 ComplexColumnVector |
|
|
3214 ComplexMatrix::column (int i) const |
|
|
3215 { |
|
|
3216 if (i < 0 || i >= nc) |
|
227
|
3217 { |
|
|
3218 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
3219 return ComplexColumnVector (); |
|
|
3220 } |
|
3
|
3221 |
|
|
3222 ComplexColumnVector retval (nr); |
|
|
3223 for (int j = 0; j < nr; j++) |
|
|
3224 retval.elem (j) = elem (j, i); |
|
|
3225 |
|
|
3226 return retval; |
|
|
3227 } |
|
|
3228 |
|
|
3229 ComplexColumnVector |
|
|
3230 ComplexMatrix::column (char *s) const |
|
|
3231 { |
|
|
3232 if (s == (char *) NULL) |
|
227
|
3233 { |
|
|
3234 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
3235 return ComplexColumnVector (); |
|
|
3236 } |
|
3
|
3237 |
|
|
3238 char c = *s; |
|
|
3239 if (c == 'f' || c == 'F') |
|
|
3240 return column (0); |
|
|
3241 else if (c == 'l' || c == 'L') |
|
|
3242 return column (nc - 1); |
|
|
3243 else |
|
227
|
3244 { |
|
|
3245 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
3246 return ComplexColumnVector (); |
|
|
3247 } |
|
3
|
3248 } |
|
|
3249 |
|
|
3250 ComplexMatrix |
|
|
3251 ComplexMatrix::inverse (int& info, double& rcond) const |
|
|
3252 { |
|
|
3253 if (nr != nc) |
|
227
|
3254 { |
|
|
3255 (*current_liboctave_error_handler) ("inverse requires square matrix"); |
|
|
3256 return ComplexMatrix (); |
|
|
3257 } |
|
3
|
3258 |
|
|
3259 info = 0; |
|
|
3260 |
|
|
3261 int *ipvt = new int [nr]; |
|
|
3262 Complex *z = new Complex [nr]; |
|
|
3263 Complex *tmp_data = dup (data, len); |
|
|
3264 |
|
|
3265 F77_FCN (zgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z); |
|
|
3266 |
|
|
3267 if (rcond + 1.0 == 1.0) |
|
|
3268 { |
|
|
3269 info = -1; |
|
|
3270 copy (tmp_data, data, len); // Restore contents. |
|
|
3271 } |
|
|
3272 else |
|
|
3273 { |
|
|
3274 int job = 1; |
|
|
3275 Complex dummy; |
|
|
3276 |
|
|
3277 F77_FCN (zgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job); |
|
|
3278 } |
|
|
3279 |
|
|
3280 delete [] ipvt; |
|
|
3281 delete [] z; |
|
|
3282 |
|
|
3283 return ComplexMatrix (tmp_data, nr, nc); |
|
|
3284 } |
|
|
3285 |
|
|
3286 ComplexMatrix |
|
|
3287 ComplexMatrix::inverse (int& info) const |
|
|
3288 { |
|
|
3289 double rcond; |
|
|
3290 return inverse (info, rcond); |
|
|
3291 } |
|
|
3292 |
|
|
3293 ComplexMatrix |
|
|
3294 ComplexMatrix::inverse (void) const |
|
|
3295 { |
|
|
3296 int info; |
|
|
3297 double rcond; |
|
|
3298 return inverse (info, rcond); |
|
|
3299 } |
|
|
3300 |
|
|
3301 ComplexMatrix |
|
|
3302 ComplexMatrix::fourier (void) const |
|
|
3303 { |
|
|
3304 int npts, nsamples; |
|
|
3305 if (nr == 1 || nc == 1) |
|
|
3306 { |
|
|
3307 npts = nr > nc ? nr : nc; |
|
|
3308 nsamples = 1; |
|
|
3309 } |
|
|
3310 else |
|
|
3311 { |
|
|
3312 npts = nr; |
|
|
3313 nsamples = nc; |
|
|
3314 } |
|
|
3315 |
|
|
3316 int nn = 4*npts+15; |
|
|
3317 Complex *wsave = new Complex [nn]; |
|
|
3318 Complex *tmp_data = dup (data, len); |
|
|
3319 |
|
|
3320 F77_FCN (cffti) (&npts, wsave); |
|
|
3321 |
|
|
3322 for (int j = 0; j < nsamples; j++) |
|
|
3323 F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave); |
|
|
3324 |
|
|
3325 delete [] wsave; |
|
|
3326 |
|
|
3327 return ComplexMatrix (tmp_data, nr, nc); |
|
|
3328 } |
|
|
3329 |
|
|
3330 ComplexMatrix |
|
|
3331 ComplexMatrix::ifourier (void) const |
|
|
3332 { |
|
|
3333 int npts, nsamples; |
|
|
3334 if (nr == 1 || nc == 1) |
|
|
3335 { |
|
|
3336 npts = nr > nc ? nr : nc; |
|
|
3337 nsamples = 1; |
|
|
3338 } |
|
|
3339 else |
|
|
3340 { |
|
|
3341 npts = nr; |
|
|
3342 nsamples = nc; |
|
|
3343 } |
|
|
3344 |
|
|
3345 int nn = 4*npts+15; |
|
|
3346 Complex *wsave = new Complex [nn]; |
|
|
3347 Complex *tmp_data = dup (data, len); |
|
|
3348 |
|
|
3349 F77_FCN (cffti) (&npts, wsave); |
|
|
3350 |
|
|
3351 for (int j = 0; j < nsamples; j++) |
|
|
3352 F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave); |
|
|
3353 |
|
|
3354 for (j = 0; j < npts*nsamples; j++) |
|
|
3355 tmp_data[j] = tmp_data[j] / (double) npts; |
|
|
3356 |
|
|
3357 delete [] wsave; |
|
|
3358 |
|
|
3359 return ComplexMatrix (tmp_data, nr, nc); |
|
|
3360 } |
|
|
3361 |
|
|
3362 ComplexDET |
|
|
3363 ComplexMatrix::determinant (void) const |
|
|
3364 { |
|
|
3365 int info; |
|
|
3366 double rcond; |
|
|
3367 return determinant (info, rcond); |
|
|
3368 } |
|
|
3369 |
|
|
3370 ComplexDET |
|
|
3371 ComplexMatrix::determinant (int& info) const |
|
|
3372 { |
|
|
3373 double rcond; |
|
|
3374 return determinant (info, rcond); |
|
|
3375 } |
|
|
3376 |
|
|
3377 ComplexDET |
|
|
3378 ComplexMatrix::determinant (int& info, double& rcond) const |
|
|
3379 { |
|
|
3380 ComplexDET retval; |
|
|
3381 |
|
|
3382 if (nr == 0 || nc == 0) |
|
|
3383 { |
|
|
3384 Complex d[2]; |
|
|
3385 d[0] = 1.0; |
|
|
3386 d[1] = 0.0; |
|
|
3387 return ComplexDET (d); |
|
|
3388 } |
|
|
3389 |
|
|
3390 info = 0; |
|
|
3391 int *ipvt = new int [nr]; |
|
|
3392 |
|
|
3393 Complex *z = new Complex [nr]; |
|
|
3394 Complex *tmp_data = dup (data, len); |
|
|
3395 |
|
|
3396 F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
3397 |
|
|
3398 if (rcond + 1.0 == 1.0) |
|
|
3399 { |
|
|
3400 info = -1; |
|
|
3401 } |
|
|
3402 else |
|
|
3403 { |
|
|
3404 int job = 10; |
|
|
3405 Complex d[2]; |
|
|
3406 F77_FCN (zgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job); |
|
|
3407 retval = ComplexDET (d); |
|
|
3408 } |
|
|
3409 |
|
|
3410 delete [] tmp_data; |
|
|
3411 delete [] ipvt; |
|
|
3412 delete [] z; |
|
|
3413 |
|
|
3414 return retval; |
|
|
3415 } |
|
|
3416 |
|
|
3417 ComplexMatrix |
|
|
3418 ComplexMatrix::solve (const Matrix& b) const |
|
|
3419 { |
|
|
3420 int info; |
|
|
3421 double rcond; |
|
|
3422 return solve (b, info, rcond); |
|
|
3423 } |
|
|
3424 |
|
|
3425 ComplexMatrix |
|
|
3426 ComplexMatrix::solve (const Matrix& b, int& info) const |
|
|
3427 { |
|
|
3428 double rcond; |
|
|
3429 return solve (b, info, rcond); |
|
|
3430 } |
|
|
3431 |
|
|
3432 ComplexMatrix |
|
|
3433 ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const |
|
|
3434 { |
|
|
3435 ComplexMatrix tmp (b); |
|
|
3436 return solve (tmp, info, rcond); |
|
|
3437 } |
|
|
3438 |
|
|
3439 ComplexMatrix |
|
|
3440 ComplexMatrix::solve (const ComplexMatrix& b) const |
|
|
3441 { |
|
|
3442 int info; |
|
|
3443 double rcond; |
|
|
3444 return solve (b, info, rcond); |
|
|
3445 } |
|
|
3446 |
|
|
3447 ComplexMatrix |
|
|
3448 ComplexMatrix::solve (const ComplexMatrix& b, int& info) const |
|
|
3449 { |
|
|
3450 double rcond; |
|
|
3451 return solve (b, info, rcond); |
|
|
3452 } |
|
|
3453 |
|
|
3454 ComplexMatrix |
|
|
3455 ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const |
|
|
3456 { |
|
|
3457 ComplexMatrix retval; |
|
|
3458 |
|
|
3459 if (nr == 0 || nc == 0 || nr != nc || nr != b.nr) |
|
227
|
3460 { |
|
|
3461 (*current_liboctave_error_handler) |
|
|
3462 ("matrix dimension mismatch in solution of linear equations"); |
|
|
3463 return ComplexMatrix (); |
|
|
3464 } |
|
3
|
3465 |
|
|
3466 info = 0; |
|
|
3467 int *ipvt = new int [nr]; |
|
|
3468 |
|
|
3469 Complex *z = new Complex [nr]; |
|
|
3470 Complex *tmp_data = dup (data, len); |
|
|
3471 |
|
|
3472 F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
3473 |
|
|
3474 if (rcond + 1.0 == 1.0) |
|
|
3475 { |
|
|
3476 info = -2; |
|
|
3477 } |
|
|
3478 else |
|
|
3479 { |
|
|
3480 int job = 0; |
|
|
3481 |
|
|
3482 Complex *result = dup (b.data, b.len); |
|
|
3483 |
|
|
3484 for (int j = 0; j < b.nc; j++) |
|
|
3485 F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job); |
|
|
3486 |
|
|
3487 retval = ComplexMatrix (result, b.nr, b.nc); |
|
|
3488 } |
|
|
3489 |
|
|
3490 delete [] tmp_data; |
|
|
3491 delete [] ipvt; |
|
|
3492 delete [] z; |
|
|
3493 |
|
|
3494 return retval; |
|
|
3495 } |
|
|
3496 |
|
|
3497 ComplexColumnVector |
|
|
3498 ComplexMatrix::solve (const ColumnVector& b) const |
|
|
3499 { |
|
|
3500 int info; |
|
|
3501 double rcond; |
|
|
3502 return solve (b, info, rcond); |
|
|
3503 } |
|
|
3504 |
|
|
3505 ComplexColumnVector |
|
|
3506 ComplexMatrix::solve (const ColumnVector& b, int& info) const |
|
|
3507 { |
|
|
3508 double rcond; |
|
|
3509 return solve (b, info, rcond); |
|
|
3510 } |
|
|
3511 |
|
|
3512 ComplexColumnVector |
|
|
3513 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond) const |
|
|
3514 { |
|
|
3515 ComplexColumnVector tmp (b); |
|
|
3516 return solve (tmp, info, rcond); |
|
|
3517 } |
|
|
3518 |
|
|
3519 ComplexColumnVector |
|
|
3520 ComplexMatrix::solve (const ComplexColumnVector& b) const |
|
|
3521 { |
|
|
3522 int info; |
|
|
3523 double rcond; |
|
|
3524 return solve (b, info, rcond); |
|
|
3525 } |
|
|
3526 |
|
|
3527 ComplexColumnVector |
|
|
3528 ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const |
|
|
3529 { |
|
|
3530 double rcond; |
|
|
3531 return solve (b, info, rcond); |
|
|
3532 } |
|
|
3533 |
|
|
3534 ComplexColumnVector |
|
|
3535 ComplexMatrix::solve (const ComplexColumnVector& b, int& info, |
|
|
3536 double& rcond) const |
|
|
3537 { |
|
|
3538 ComplexColumnVector retval; |
|
|
3539 |
|
|
3540 if (nr == 0 || nc == 0 || nr != nc || nr != b.len) |
|
227
|
3541 { |
|
|
3542 (*current_liboctave_error_handler) |
|
|
3543 ("matrix dimension mismatch in solution of linear equations"); |
|
|
3544 return ComplexColumnVector (); |
|
|
3545 } |
|
3
|
3546 |
|
|
3547 info = 0; |
|
|
3548 int *ipvt = new int [nr]; |
|
|
3549 |
|
|
3550 Complex *z = new Complex [nr]; |
|
|
3551 Complex *tmp_data = dup (data, len); |
|
|
3552 |
|
|
3553 F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); |
|
|
3554 |
|
|
3555 if (rcond + 1.0 == 1.0) |
|
|
3556 { |
|
|
3557 info = -2; |
|
|
3558 } |
|
|
3559 else |
|
|
3560 { |
|
|
3561 int job = 0; |
|
|
3562 |
|
|
3563 Complex *result = dup (b.data, b.len); |
|
|
3564 |
|
|
3565 F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, result, &job); |
|
|
3566 |
|
|
3567 retval = ComplexColumnVector (result, b.len); |
|
|
3568 } |
|
|
3569 |
|
|
3570 delete [] tmp_data; |
|
|
3571 delete [] ipvt; |
|
|
3572 delete [] z; |
|
|
3573 |
|
|
3574 return retval; |
|
|
3575 } |
|
|
3576 |
|
|
3577 ComplexMatrix |
|
|
3578 ComplexMatrix::lssolve (const Matrix& b) const |
|
|
3579 { |
|
|
3580 int info; |
|
|
3581 int rank; |
|
|
3582 return lssolve (b, info, rank); |
|
|
3583 } |
|
|
3584 |
|
|
3585 ComplexMatrix |
|
|
3586 ComplexMatrix::lssolve (const Matrix& b, int& info) const |
|
|
3587 { |
|
|
3588 int rank; |
|
|
3589 return lssolve (b, info, rank); |
|
|
3590 } |
|
|
3591 |
|
|
3592 ComplexMatrix |
|
|
3593 ComplexMatrix::lssolve (const Matrix& b, int& info, int& rank) const |
|
|
3594 { |
|
|
3595 ComplexMatrix tmp (b); |
|
|
3596 return lssolve (tmp, info, rank); |
|
|
3597 } |
|
|
3598 |
|
|
3599 ComplexMatrix |
|
|
3600 ComplexMatrix::lssolve (const ComplexMatrix& b) const |
|
|
3601 { |
|
|
3602 int info; |
|
|
3603 int rank; |
|
|
3604 return lssolve (b, info, rank); |
|
|
3605 } |
|
|
3606 |
|
|
3607 ComplexMatrix |
|
|
3608 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const |
|
|
3609 { |
|
|
3610 int rank; |
|
|
3611 return lssolve (b, info, rank); |
|
|
3612 } |
|
|
3613 |
|
|
3614 ComplexMatrix |
|
|
3615 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const |
|
|
3616 { |
|
|
3617 int nrhs = b.nc; |
|
|
3618 |
|
|
3619 int m = nr; |
|
|
3620 int n = nc; |
|
|
3621 |
|
|
3622 if (m == 0 || n == 0 || m != b.nr) |
|
227
|
3623 { |
|
|
3624 (*current_liboctave_error_handler) |
|
|
3625 ("matrix dimension mismatch solution of linear equations"); |
|
|
3626 return Matrix (); |
|
|
3627 } |
|
3
|
3628 |
|
|
3629 Complex *tmp_data = dup (data, len); |
|
|
3630 |
|
|
3631 int nrr = m > n ? m : n; |
|
|
3632 ComplexMatrix result (nrr, nrhs); |
|
|
3633 |
|
|
3634 int i, j; |
|
|
3635 for (j = 0; j < nrhs; j++) |
|
|
3636 for (i = 0; i < m; i++) |
|
|
3637 result.elem (i, j) = b.elem (i, j); |
|
|
3638 |
|
|
3639 Complex *presult = result.fortran_vec (); |
|
|
3640 |
|
|
3641 int len_s = m < n ? m : n; |
|
|
3642 double *s = new double [len_s]; |
|
|
3643 double rcond = -1.0; |
|
|
3644 int lwork; |
|
|
3645 if (m < n) |
|
|
3646 lwork = 2*m + (nrhs > n ? nrhs : n); |
|
|
3647 else |
|
|
3648 lwork = 2*n + (nrhs > m ? nrhs : m); |
|
|
3649 |
|
|
3650 Complex *work = new Complex [lwork]; |
|
|
3651 |
|
|
3652 int lrwork = (5 * (m < n ? m : n)) - 4; |
|
|
3653 lrwork = lrwork > 1 ? lrwork : 1; |
|
|
3654 double *rwork = new double [lrwork]; |
|
|
3655 |
|
|
3656 F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, |
|
|
3657 &rcond, &rank, work, &lwork, rwork, &info); |
|
|
3658 |
|
|
3659 ComplexMatrix retval (n, nrhs); |
|
|
3660 for (j = 0; j < nrhs; j++) |
|
|
3661 for (i = 0; i < n; i++) |
|
|
3662 retval.elem (i, j) = result.elem (i, j); |
|
|
3663 |
|
|
3664 delete [] tmp_data; |
|
|
3665 delete [] s; |
|
|
3666 delete [] work; |
|
|
3667 delete [] rwork; |
|
|
3668 |
|
|
3669 return retval; |
|
|
3670 } |
|
|
3671 |
|
|
3672 ComplexColumnVector |
|
|
3673 ComplexMatrix::lssolve (const ColumnVector& b) const |
|
|
3674 { |
|
|
3675 int info; |
|
|
3676 int rank; |
|
|
3677 return lssolve (b, info, rank); |
|
|
3678 } |
|
|
3679 |
|
|
3680 ComplexColumnVector |
|
|
3681 ComplexMatrix::lssolve (const ColumnVector& b, int& info) const |
|
|
3682 { |
|
|
3683 int rank; |
|
|
3684 return lssolve (b, info, rank); |
|
|
3685 } |
|
|
3686 |
|
|
3687 ComplexColumnVector |
|
|
3688 ComplexMatrix::lssolve (const ColumnVector& b, int& info, int& rank) const |
|
|
3689 { |
|
|
3690 ComplexColumnVector tmp (b); |
|
|
3691 return lssolve (tmp, info, rank); |
|
|
3692 } |
|
|
3693 |
|
|
3694 ComplexColumnVector |
|
|
3695 ComplexMatrix::lssolve (const ComplexColumnVector& b) const |
|
|
3696 { |
|
|
3697 int info; |
|
|
3698 int rank; |
|
|
3699 return lssolve (b, info, rank); |
|
|
3700 } |
|
|
3701 |
|
|
3702 ComplexColumnVector |
|
|
3703 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const |
|
|
3704 { |
|
|
3705 int rank; |
|
|
3706 return lssolve (b, info, rank); |
|
|
3707 } |
|
|
3708 |
|
|
3709 ComplexColumnVector |
|
|
3710 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info, |
|
|
3711 int& rank) const |
|
|
3712 { |
|
|
3713 int nrhs = 1; |
|
|
3714 |
|
|
3715 int m = nr; |
|
|
3716 int n = nc; |
|
|
3717 |
|
|
3718 if (m == 0 || n == 0 || m != b.len) |
|
227
|
3719 { |
|
|
3720 (*current_liboctave_error_handler) |
|
|
3721 ("matrix dimension mismatch solution of least squares problem"); |
|
|
3722 return ComplexColumnVector (); |
|
|
3723 } |
|
3
|
3724 |
|
|
3725 Complex *tmp_data = dup (data, len); |
|
|
3726 |
|
|
3727 int nrr = m > n ? m : n; |
|
|
3728 ComplexColumnVector result (nrr); |
|
|
3729 |
|
|
3730 int i; |
|
|
3731 for (i = 0; i < m; i++) |
|
|
3732 result.elem (i) = b.elem (i); |
|
|
3733 |
|
|
3734 Complex *presult = result.fortran_vec (); |
|
|
3735 |
|
|
3736 int len_s = m < n ? m : n; |
|
|
3737 double *s = new double [len_s]; |
|
|
3738 double rcond = -1.0; |
|
|
3739 int lwork; |
|
|
3740 if (m < n) |
|
|
3741 lwork = 2*m + (nrhs > n ? nrhs : n); |
|
|
3742 else |
|
|
3743 lwork = 2*n + (nrhs > m ? nrhs : m); |
|
|
3744 |
|
|
3745 Complex *work = new Complex [lwork]; |
|
|
3746 |
|
|
3747 int lrwork = (5 * (m < n ? m : n)) - 4; |
|
|
3748 lrwork = lrwork > 1 ? lrwork : 1; |
|
|
3749 double *rwork = new double [lrwork]; |
|
|
3750 |
|
|
3751 F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, |
|
|
3752 &rcond, &rank, work, &lwork, rwork, &info); |
|
|
3753 |
|
|
3754 ComplexColumnVector retval (n); |
|
|
3755 for (i = 0; i < n; i++) |
|
|
3756 retval.elem (i) = result.elem (i); |
|
|
3757 |
|
|
3758 delete [] tmp_data; |
|
|
3759 delete [] s; |
|
|
3760 delete [] work; |
|
|
3761 delete [] rwork; |
|
|
3762 |
|
|
3763 return retval; |
|
|
3764 } |
|
|
3765 |
|
|
3766 // matrix by scalar -> matrix operations |
|
|
3767 |
|
|
3768 ComplexMatrix |
|
|
3769 ComplexMatrix::operator + (double s) const |
|
|
3770 { |
|
|
3771 return ComplexMatrix (add (data, len, s), nr, nc); |
|
|
3772 } |
|
|
3773 |
|
|
3774 ComplexMatrix |
|
|
3775 ComplexMatrix::operator - (double s) const |
|
|
3776 { |
|
|
3777 return ComplexMatrix (subtract (data, len, s), nr, nc); |
|
|
3778 } |
|
|
3779 |
|
|
3780 ComplexMatrix |
|
|
3781 ComplexMatrix::operator * (double s) const |
|
|
3782 { |
|
|
3783 return ComplexMatrix (multiply (data, len, s), nr, nc); |
|
|
3784 } |
|
|
3785 |
|
|
3786 ComplexMatrix |
|
|
3787 ComplexMatrix::operator / (double s) const |
|
|
3788 { |
|
|
3789 return ComplexMatrix (divide (data, len, s), nr, nc); |
|
|
3790 } |
|
|
3791 |
|
|
3792 ComplexMatrix |
|
161
|
3793 ComplexMatrix::operator + (const Complex& s) const |
|
3
|
3794 { |
|
|
3795 return ComplexMatrix (add (data, len, s), nr, nc); |
|
|
3796 } |
|
|
3797 |
|
|
3798 ComplexMatrix |
|
161
|
3799 ComplexMatrix::operator - (const Complex& s) const |
|
3
|
3800 { |
|
|
3801 return ComplexMatrix (subtract (data, len, s), nr, nc); |
|
|
3802 } |
|
|
3803 |
|
|
3804 ComplexMatrix |
|
161
|
3805 ComplexMatrix::operator * (const Complex& s) const |
|
3
|
3806 { |
|
|
3807 return ComplexMatrix (multiply (data, len, s), nr, nc); |
|
|
3808 } |
|
|
3809 |
|
|
3810 ComplexMatrix |
|
161
|
3811 ComplexMatrix::operator / (const Complex& s) const |
|
3
|
3812 { |
|
|
3813 return ComplexMatrix (divide (data, len, s), nr, nc); |
|
|
3814 } |
|
|
3815 |
|
|
3816 // scalar by matrix -> matrix operations |
|
|
3817 |
|
|
3818 ComplexMatrix |
|
|
3819 operator + (double s, const ComplexMatrix& a) |
|
|
3820 { |
|
|
3821 return ComplexMatrix (add (a.data, a.len, s), a.nr, a.nc); |
|
|
3822 } |
|
|
3823 |
|
|
3824 ComplexMatrix |
|
|
3825 operator - (double s, const ComplexMatrix& a) |
|
|
3826 { |
|
|
3827 return ComplexMatrix (subtract (s, a.data, a.len), a.nr, a.nc); |
|
|
3828 } |
|
|
3829 |
|
|
3830 ComplexMatrix |
|
|
3831 operator * (double s, const ComplexMatrix& a) |
|
|
3832 { |
|
|
3833 return ComplexMatrix (multiply (a.data, a.len, s), a.nr, a.nc); |
|
|
3834 } |
|
|
3835 |
|
|
3836 ComplexMatrix |
|
|
3837 operator / (double s, const ComplexMatrix& a) |
|
|
3838 { |
|
|
3839 return ComplexMatrix (divide (s, a.data, a.len), a.nr, a.nc); |
|
|
3840 } |
|
|
3841 |
|
|
3842 ComplexMatrix |
|
161
|
3843 operator + (const Complex& s, const ComplexMatrix& a) |
|
3
|
3844 { |
|
|
3845 return ComplexMatrix (add (s, a.data, a.len), a.nr, a.nc); |
|
|
3846 } |
|
|
3847 |
|
|
3848 ComplexMatrix |
|
161
|
3849 operator - (const Complex& s, const ComplexMatrix& a) |
|
3
|
3850 { |
|
|
3851 return ComplexMatrix (subtract (s, a.data, a.len), a.nr, a.nc); |
|
|
3852 } |
|
|
3853 |
|
|
3854 ComplexMatrix |
|
161
|
3855 operator * (const Complex& s, const ComplexMatrix& a) |
|
3
|
3856 { |
|
|
3857 return ComplexMatrix (multiply (s, a.data, a.len), a.nr, a.nc); |
|
|
3858 } |
|
|
3859 |
|
|
3860 ComplexMatrix |
|
161
|
3861 operator / (const Complex& s, const ComplexMatrix& a) |
|
3
|
3862 { |
|
|
3863 return ComplexMatrix (divide (s, a.data, a.len), a.nr, a.nc); |
|
|
3864 } |
|
|
3865 |
|
|
3866 // matrix by column vector -> column vector operations |
|
|
3867 |
|
|
3868 ComplexColumnVector |
|
|
3869 ComplexMatrix::operator * (const ColumnVector& a) const |
|
|
3870 { |
|
|
3871 ComplexColumnVector tmp (a); |
|
|
3872 return *this * tmp; |
|
|
3873 } |
|
|
3874 |
|
|
3875 ComplexColumnVector |
|
|
3876 ComplexMatrix::operator * (const ComplexColumnVector& a) const |
|
|
3877 { |
|
|
3878 if (nc != a.len) |
|
227
|
3879 { |
|
|
3880 (*current_liboctave_error_handler) |
|
|
3881 ("nonconformant matrix multiplication attempted"); |
|
|
3882 return ComplexColumnVector (); |
|
|
3883 } |
|
3
|
3884 |
|
|
3885 if (nc == 0 || nr == 0) |
|
|
3886 return ComplexColumnVector (0); |
|
|
3887 |
|
|
3888 char trans = 'N'; |
|
|
3889 int ld = nr; |
|
|
3890 Complex alpha (1.0); |
|
|
3891 Complex beta (0.0); |
|
|
3892 int i_one = 1; |
|
|
3893 |
|
125
|
3894 Complex *y = new Complex [nr]; |
|
3
|
3895 |
|
|
3896 F77_FCN (zgemv) (&trans, &nr, &nc, &alpha, data, &ld, a.data, |
|
|
3897 &i_one, &beta, y, &i_one, 1L); |
|
|
3898 |
|
124
|
3899 return ComplexColumnVector (y, nr); |
|
3
|
3900 } |
|
|
3901 |
|
|
3902 // matrix by diagonal matrix -> matrix operations |
|
|
3903 |
|
|
3904 ComplexMatrix |
|
|
3905 ComplexMatrix::operator + (const DiagMatrix& a) const |
|
|
3906 { |
|
|
3907 if (nr != a.nr || nc != a.nc) |
|
227
|
3908 { |
|
|
3909 (*current_liboctave_error_handler) |
|
|
3910 ("nonconformant matrix addition attempted"); |
|
|
3911 return ComplexMatrix (); |
|
|
3912 } |
|
3
|
3913 |
|
|
3914 if (nr == 0 || nc == 0) |
|
|
3915 return ComplexMatrix (nr, nc); |
|
|
3916 |
|
|
3917 ComplexMatrix result (*this); |
|
|
3918 for (int i = 0; i < a.len; i++) |
|
|
3919 result.elem (i, i) += a.data[i]; |
|
|
3920 |
|
|
3921 return result; |
|
|
3922 } |
|
|
3923 |
|
|
3924 ComplexMatrix |
|
|
3925 ComplexMatrix::operator - (const DiagMatrix& a) const |
|
|
3926 { |
|
|
3927 if (nr != a.nr || nc != a.nc) |
|
227
|
3928 { |
|
|
3929 (*current_liboctave_error_handler) |
|
|
3930 ("nonconformant matrix subtraction attempted"); |
|
|
3931 return ComplexMatrix (); |
|
|
3932 } |
|
3
|
3933 |
|
|
3934 if (nr == 0 || nc == 0) |
|
|
3935 return ComplexMatrix (nr, nc); |
|
|
3936 |
|
|
3937 ComplexMatrix result (*this); |
|
|
3938 for (int i = 0; i < a.len; i++) |
|
|
3939 result.elem (i, i) -= a.data[i]; |
|
|
3940 |
|
|
3941 return result; |
|
|
3942 } |
|
|
3943 |
|
|
3944 ComplexMatrix |
|
|
3945 ComplexMatrix::operator * (const DiagMatrix& a) const |
|
|
3946 { |
|
|
3947 if (nc != a.nr) |
|
227
|
3948 { |
|
|
3949 (*current_liboctave_error_handler) |
|
|
3950 ("nonconformant matrix multiplication attempted"); |
|
|
3951 return ComplexMatrix (); |
|
|
3952 } |
|
3
|
3953 |
|
|
3954 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
3955 return ComplexMatrix (nr, nc, 0.0); |
|
|
3956 |
|
|
3957 Complex *c = new Complex [nr*a.nc]; |
|
|
3958 Complex *ctmp = (Complex *) NULL; |
|
|
3959 |
|
|
3960 for (int j = 0; j < a.len; j++) |
|
|
3961 { |
|
|
3962 int idx = j * nr; |
|
|
3963 ctmp = c + idx; |
|
|
3964 if (a.data[j] == 1.0) |
|
|
3965 { |
|
|
3966 for (int i = 0; i < nr; i++) |
|
|
3967 ctmp[i] = elem (i, j); |
|
|
3968 } |
|
|
3969 else if (a.data[j] == 0.0) |
|
|
3970 { |
|
|
3971 for (int i = 0; i < nr; i++) |
|
|
3972 ctmp[i] = 0.0; |
|
|
3973 } |
|
|
3974 else |
|
|
3975 { |
|
|
3976 for (int i = 0; i < nr; i++) |
|
|
3977 ctmp[i] = a.data[j] * elem (i, j); |
|
|
3978 } |
|
|
3979 } |
|
|
3980 |
|
|
3981 if (a.nr < a.nc) |
|
|
3982 { |
|
|
3983 for (int i = nr * nc; i < nr * a.nc; i++) |
|
|
3984 ctmp[i] = 0.0; |
|
|
3985 } |
|
|
3986 |
|
|
3987 return ComplexMatrix (c, nr, a.nc); |
|
|
3988 } |
|
|
3989 |
|
|
3990 ComplexMatrix |
|
|
3991 ComplexMatrix::operator + (const ComplexDiagMatrix& a) const |
|
|
3992 { |
|
|
3993 if (nr != a.nr || nc != a.nc) |
|
227
|
3994 { |
|
|
3995 (*current_liboctave_error_handler) |
|
|
3996 ("nonconformant matrix addition attempted"); |
|
|
3997 return ComplexMatrix (); |
|
|
3998 } |
|
3
|
3999 |
|
|
4000 if (nr == 0 || nc == 0) |
|
|
4001 return ComplexMatrix (nr, nc); |
|
|
4002 |
|
|
4003 ComplexMatrix result (*this); |
|
|
4004 for (int i = 0; i < a.len; i++) |
|
|
4005 result.elem (i, i) += a.data[i]; |
|
|
4006 |
|
|
4007 return result; |
|
|
4008 } |
|
|
4009 |
|
|
4010 ComplexMatrix |
|
|
4011 ComplexMatrix::operator - (const ComplexDiagMatrix& a) const |
|
|
4012 { |
|
|
4013 if (nr != a.nr || nc != a.nc) |
|
227
|
4014 { |
|
|
4015 (*current_liboctave_error_handler) |
|
|
4016 ("nonconformant matrix subtraction attempted"); |
|
|
4017 return ComplexMatrix (); |
|
|
4018 } |
|
3
|
4019 |
|
|
4020 if (nr == 0 || nc == 0) |
|
|
4021 return ComplexMatrix (nr, nc); |
|
|
4022 |
|
|
4023 ComplexMatrix result (*this); |
|
|
4024 for (int i = 0; i < a.len; i++) |
|
|
4025 result.elem (i, i) -= a.data[i]; |
|
|
4026 |
|
|
4027 return result; |
|
|
4028 } |
|
|
4029 |
|
|
4030 ComplexMatrix |
|
|
4031 ComplexMatrix::operator * (const ComplexDiagMatrix& a) const |
|
|
4032 { |
|
|
4033 if (nc != a.nr) |
|
227
|
4034 { |
|
|
4035 (*current_liboctave_error_handler) |
|
|
4036 ("nonconformant matrix multiplication attempted"); |
|
|
4037 return ComplexMatrix (); |
|
|
4038 } |
|
3
|
4039 |
|
|
4040 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
4041 return ComplexMatrix (nr, nc, 0.0); |
|
|
4042 |
|
|
4043 Complex *c = new Complex [nr*a.nc]; |
|
|
4044 Complex *ctmp = (Complex *) NULL; |
|
|
4045 |
|
|
4046 for (int j = 0; j < a.len; j++) |
|
|
4047 { |
|
|
4048 int idx = j * nr; |
|
|
4049 ctmp = c + idx; |
|
|
4050 if (a.data[j] == 1.0) |
|
|
4051 { |
|
|
4052 for (int i = 0; i < nr; i++) |
|
|
4053 ctmp[i] = elem (i, j); |
|
|
4054 } |
|
|
4055 else if (a.data[j] == 0.0) |
|
|
4056 { |
|
|
4057 for (int i = 0; i < nr; i++) |
|
|
4058 ctmp[i] = 0.0; |
|
|
4059 } |
|
|
4060 else |
|
|
4061 { |
|
|
4062 for (int i = 0; i < nr; i++) |
|
|
4063 ctmp[i] = a.data[j] * elem (i, j); |
|
|
4064 } |
|
|
4065 } |
|
|
4066 |
|
|
4067 if (a.nr < a.nc) |
|
|
4068 { |
|
|
4069 for (int i = nr * nc; i < nr * a.nc; i++) |
|
|
4070 ctmp[i] = 0.0; |
|
|
4071 } |
|
|
4072 |
|
|
4073 return ComplexMatrix (c, nr, a.nc); |
|
|
4074 } |
|
|
4075 |
|
|
4076 ComplexMatrix& |
|
|
4077 ComplexMatrix::operator += (const DiagMatrix& a) |
|
|
4078 { |
|
|
4079 if (nr != a.nr || nc != a.nc) |
|
227
|
4080 { |
|
|
4081 (*current_liboctave_error_handler) |
|
|
4082 ("nonconformant matrix += operation attempted"); |
|
|
4083 return ComplexMatrix (); |
|
|
4084 } |
|
3
|
4085 |
|
|
4086 for (int i = 0; i < a.len; i++) |
|
|
4087 elem (i, i) += a.data[i]; |
|
|
4088 |
|
|
4089 return *this; |
|
|
4090 } |
|
|
4091 |
|
|
4092 ComplexMatrix& |
|
|
4093 ComplexMatrix::operator -= (const DiagMatrix& a) |
|
|
4094 { |
|
|
4095 if (nr != a.nr || nc != a.nc) |
|
227
|
4096 { |
|
|
4097 (*current_liboctave_error_handler) |
|
|
4098 ("nonconformant matrix -= operation attempted"); |
|
|
4099 return ComplexMatrix (); |
|
|
4100 } |
|
3
|
4101 |
|
|
4102 for (int i = 0; i < a.len; i++) |
|
|
4103 elem (i, i) -= a.data[i]; |
|
|
4104 |
|
|
4105 return *this; |
|
|
4106 } |
|
|
4107 |
|
|
4108 ComplexMatrix& |
|
|
4109 ComplexMatrix::operator += (const ComplexDiagMatrix& a) |
|
|
4110 { |
|
|
4111 if (nr != a.nr || nc != a.nc) |
|
227
|
4112 { |
|
|
4113 (*current_liboctave_error_handler) |
|
|
4114 ("nonconformant matrix += operation attempted"); |
|
|
4115 return ComplexMatrix (); |
|
|
4116 } |
|
3
|
4117 |
|
|
4118 for (int i = 0; i < a.len; i++) |
|
|
4119 elem (i, i) += a.data[i]; |
|
|
4120 |
|
|
4121 return *this; |
|
|
4122 } |
|
|
4123 |
|
|
4124 ComplexMatrix& |
|
|
4125 ComplexMatrix::operator -= (const ComplexDiagMatrix& a) |
|
|
4126 { |
|
|
4127 if (nr != a.nr || nc != a.nc) |
|
227
|
4128 { |
|
|
4129 (*current_liboctave_error_handler) |
|
|
4130 ("nonconformant matrix -= operation attempted"); |
|
|
4131 return ComplexMatrix (); |
|
|
4132 } |
|
3
|
4133 |
|
|
4134 for (int i = 0; i < a.len; i++) |
|
|
4135 elem (i, i) -= a.data[i]; |
|
|
4136 |
|
|
4137 return *this; |
|
|
4138 } |
|
|
4139 |
|
|
4140 // matrix by matrix -> matrix operations |
|
|
4141 |
|
|
4142 ComplexMatrix |
|
|
4143 ComplexMatrix::operator + (const Matrix& a) const |
|
|
4144 { |
|
|
4145 if (nr != a.nr || nc != a.nc) |
|
227
|
4146 { |
|
|
4147 (*current_liboctave_error_handler) |
|
|
4148 ("nonconformant matrix addition attempted"); |
|
|
4149 return ComplexMatrix (); |
|
|
4150 } |
|
3
|
4151 |
|
|
4152 if (nr == 0 || nc == 0) |
|
|
4153 return ComplexMatrix (nr, nc); |
|
|
4154 |
|
|
4155 return ComplexMatrix (add (data, a.data, len), nr, nc); |
|
|
4156 } |
|
|
4157 |
|
|
4158 ComplexMatrix |
|
|
4159 ComplexMatrix::operator - (const Matrix& a) const |
|
|
4160 { |
|
|
4161 if (nr != a.nr || nc != a.nc) |
|
227
|
4162 { |
|
|
4163 (*current_liboctave_error_handler) |
|
|
4164 ("nonconformant matrix subtraction attempted"); |
|
|
4165 return ComplexMatrix (); |
|
|
4166 } |
|
3
|
4167 |
|
|
4168 if (nr == 0 || nc == 0) |
|
|
4169 return ComplexMatrix (nr, nc); |
|
|
4170 |
|
|
4171 return ComplexMatrix (subtract (data, a.data, len), nr, nc); |
|
|
4172 } |
|
|
4173 |
|
|
4174 ComplexMatrix |
|
|
4175 ComplexMatrix::operator * (const Matrix& a) const |
|
|
4176 { |
|
|
4177 ComplexMatrix tmp (a); |
|
|
4178 return *this * tmp; |
|
|
4179 } |
|
|
4180 |
|
|
4181 ComplexMatrix |
|
|
4182 ComplexMatrix::operator + (const ComplexMatrix& a) const |
|
|
4183 { |
|
|
4184 if (nr != a.nr || nc != a.nc) |
|
227
|
4185 { |
|
|
4186 (*current_liboctave_error_handler) |
|
|
4187 ("nonconformant matrix addition attempted"); |
|
|
4188 return ComplexMatrix (); |
|
|
4189 } |
|
3
|
4190 |
|
|
4191 if (nr == 0 || nc == 0) |
|
|
4192 return ComplexMatrix (nr, nc); |
|
|
4193 |
|
|
4194 return ComplexMatrix (add (data, a.data, len), nr, nc); |
|
|
4195 } |
|
|
4196 |
|
|
4197 ComplexMatrix |
|
|
4198 ComplexMatrix::operator - (const ComplexMatrix& a) const |
|
|
4199 { |
|
|
4200 if (nr != a.nr || nc != a.nc) |
|
227
|
4201 { |
|
|
4202 (*current_liboctave_error_handler) |
|
|
4203 ("nonconformant matrix subtraction attempted"); |
|
|
4204 return ComplexMatrix (); |
|
|
4205 } |
|
3
|
4206 |
|
|
4207 if (nr == 0 || nc == 0) |
|
|
4208 return ComplexMatrix (nr, nc); |
|
|
4209 |
|
|
4210 return ComplexMatrix (subtract (data, a.data, len), nr, nc); |
|
|
4211 } |
|
|
4212 |
|
|
4213 ComplexMatrix |
|
|
4214 ComplexMatrix::operator * (const ComplexMatrix& a) const |
|
|
4215 { |
|
|
4216 if (nc != a.nr) |
|
227
|
4217 { |
|
|
4218 (*current_liboctave_error_handler) |
|
|
4219 ("nonconformant matrix multiplication attempted"); |
|
|
4220 return ComplexMatrix (); |
|
|
4221 } |
|
3
|
4222 |
|
|
4223 if (nr == 0 || nc == 0 || a.nc == 0) |
|
|
4224 return ComplexMatrix (nr, nc, 0.0); |
|
|
4225 |
|
|
4226 char trans = 'N'; |
|
|
4227 char transa = 'N'; |
|
|
4228 |
|
|
4229 int ld = nr; |
|
|
4230 int lda = a.nr; |
|
|
4231 |
|
|
4232 Complex alpha (1.0); |
|
|
4233 Complex beta (0.0); |
|
|
4234 int anc = a.nc; |
|
|
4235 |
|
|
4236 Complex *c = new Complex [nr*a.nc]; |
|
|
4237 |
|
|
4238 F77_FCN (zgemm) (&trans, &transa, &nr, &anc, &nc, &alpha, data, &ld, |
|
|
4239 a.data, &lda, &beta, c, &nr, 1L, 1L); |
|
|
4240 |
|
|
4241 return ComplexMatrix (c, nr, a.nc); |
|
|
4242 } |
|
|
4243 |
|
|
4244 ComplexMatrix |
|
|
4245 ComplexMatrix::product (const Matrix& a) const |
|
|
4246 { |
|
|
4247 if (nr != a.nr || nc != a.nc) |
|
227
|
4248 { |
|
|
4249 (*current_liboctave_error_handler) |
|
|
4250 ("nonconformant matrix product attempted"); |
|
|
4251 return ComplexMatrix (); |
|
|
4252 } |
|
3
|
4253 |
|
|
4254 if (nr == 0 || nc == 0) |
|
|
4255 return ComplexMatrix (nr, nc); |
|
|
4256 |
|
|
4257 return ComplexMatrix (multiply (data, a.data, len), nr, nc); |
|
|
4258 } |
|
|
4259 |
|
|
4260 ComplexMatrix |
|
|
4261 ComplexMatrix::quotient (const Matrix& a) const |
|
|
4262 { |
|
|
4263 if (nr != a.nr || nc != a.nc) |
|
227
|
4264 { |
|
|
4265 (*current_liboctave_error_handler) |
|
|
4266 ("nonconformant matrix quotient attempted"); |
|
|
4267 return ComplexMatrix (); |
|
|
4268 } |
|
3
|
4269 |
|
|
4270 if (nr == 0 || nc == 0) |
|
|
4271 return ComplexMatrix (nr, nc); |
|
|
4272 |
|
|
4273 return ComplexMatrix (divide (data, a.data, len), nr, nc); |
|
|
4274 } |
|
|
4275 |
|
|
4276 ComplexMatrix |
|
|
4277 ComplexMatrix::product (const ComplexMatrix& a) const |
|
|
4278 { |
|
|
4279 if (nr != a.nr || nc != a.nc) |
|
227
|
4280 { |
|
|
4281 (*current_liboctave_error_handler) |
|
|
4282 ("nonconformant matrix product attempted"); |
|
|
4283 return ComplexMatrix (); |
|
|
4284 } |
|
3
|
4285 |
|
|
4286 if (nr == 0 || nc == 0) |
|
|
4287 return ComplexMatrix (nr, nc); |
|
|
4288 |
|
|
4289 return ComplexMatrix (multiply (data, a.data, len), nr, nc); |
|
|
4290 } |
|
|
4291 |
|
|
4292 ComplexMatrix |
|
|
4293 ComplexMatrix::quotient (const ComplexMatrix& a) const |
|
|
4294 { |
|
|
4295 if (nr != a.nr || nc != a.nc) |
|
227
|
4296 { |
|
|
4297 (*current_liboctave_error_handler) |
|
|
4298 ("nonconformant matrix quotient attempted"); |
|
|
4299 return ComplexMatrix (); |
|
|
4300 } |
|
3
|
4301 |
|
|
4302 if (nr == 0 || nc == 0) |
|
|
4303 return ComplexMatrix (nr, nc); |
|
|
4304 |
|
|
4305 return ComplexMatrix (divide (data, a.data, len), nr, nc); |
|
|
4306 } |
|
|
4307 |
|
|
4308 ComplexMatrix& |
|
|
4309 ComplexMatrix::operator += (const Matrix& a) |
|
|
4310 { |
|
|
4311 if (nr != a.nr || nc != a.nc) |
|
227
|
4312 { |
|
|
4313 (*current_liboctave_error_handler) |
|
|
4314 ("nonconformant matrix += operation attempted"); |
|
|
4315 return *this; |
|
|
4316 } |
|
3
|
4317 |
|
|
4318 if (nr == 0 || nc == 0) |
|
|
4319 return *this; |
|
|
4320 |
|
|
4321 add2 (data, a.data, len); |
|
|
4322 return *this; |
|
|
4323 } |
|
|
4324 |
|
|
4325 ComplexMatrix& |
|
|
4326 ComplexMatrix::operator -= (const Matrix& a) |
|
|
4327 { |
|
|
4328 if (nr != a.nr || nc != a.nc) |
|
227
|
4329 { |
|
|
4330 (*current_liboctave_error_handler) |
|
|
4331 ("nonconformant matrix -= operation attempted"); |
|
|
4332 return *this; |
|
|
4333 } |
|
3
|
4334 |
|
|
4335 if (nr == 0 || nc == 0) |
|
|
4336 return *this; |
|
|
4337 |
|
|
4338 subtract2 (data, a.data, len); |
|
|
4339 return *this; |
|
|
4340 } |
|
|
4341 |
|
|
4342 ComplexMatrix& |
|
|
4343 ComplexMatrix::operator += (const ComplexMatrix& a) |
|
|
4344 { |
|
|
4345 if (nr != a.nr || nc != a.nc) |
|
227
|
4346 { |
|
|
4347 (*current_liboctave_error_handler) |
|
|
4348 ("nonconformant matrix += operation attempted"); |
|
|
4349 return *this; |
|
|
4350 } |
|
3
|
4351 |
|
|
4352 if (nr == 0 || nc == 0) |
|
|
4353 return *this; |
|
|
4354 |
|
|
4355 add2 (data, a.data, len); |
|
|
4356 return *this; |
|
|
4357 } |
|
|
4358 |
|
|
4359 ComplexMatrix& |
|
|
4360 ComplexMatrix::operator -= (const ComplexMatrix& a) |
|
|
4361 { |
|
|
4362 if (nr != a.nr || nc != a.nc) |
|
227
|
4363 { |
|
|
4364 (*current_liboctave_error_handler) |
|
|
4365 ("nonconformant matrix -= operation attempted"); |
|
|
4366 return *this; |
|
|
4367 } |
|
3
|
4368 |
|
|
4369 if (nr == 0 || nc == 0) |
|
|
4370 return *this; |
|
|
4371 |
|
|
4372 subtract2 (data, a.data, len); |
|
|
4373 return *this; |
|
|
4374 } |
|
|
4375 |
|
|
4376 // unary operations |
|
|
4377 |
|
|
4378 ComplexMatrix |
|
|
4379 ComplexMatrix::operator - (void) const |
|
|
4380 { |
|
|
4381 return ComplexMatrix (negate (data, len), nr, nc); |
|
|
4382 } |
|
|
4383 |
|
|
4384 Matrix |
|
|
4385 ComplexMatrix::operator ! (void) const |
|
|
4386 { |
|
|
4387 return Matrix (not (data, len), nr, nc); |
|
|
4388 } |
|
|
4389 |
|
|
4390 // other operations |
|
|
4391 |
|
|
4392 ComplexMatrix |
|
|
4393 map (c_c_Mapper f, const ComplexMatrix& a) |
|
|
4394 { |
|
|
4395 ComplexMatrix b (a); |
|
|
4396 b.map (f); |
|
|
4397 return b; |
|
|
4398 } |
|
|
4399 |
|
|
4400 Matrix |
|
|
4401 map (d_c_Mapper f, const ComplexMatrix& a) |
|
|
4402 { |
|
|
4403 Matrix b (a.nr, a.nc); |
|
|
4404 for (int j = 0; j < a.nc; j++) |
|
|
4405 for (int i = 0; i < a.nr; i++) |
|
|
4406 b.elem (i, j) = f (a.elem (i, j)); |
|
|
4407 return b; |
|
|
4408 } |
|
|
4409 |
|
|
4410 void |
|
|
4411 ComplexMatrix::map (c_c_Mapper f) |
|
|
4412 { |
|
|
4413 for (int i = 0; i < len; i++) |
|
|
4414 data[i] = f (data[i]); |
|
|
4415 } |
|
|
4416 |
|
|
4417 Matrix |
|
|
4418 ComplexMatrix::all (void) const |
|
|
4419 { |
|
|
4420 Matrix retval; |
|
|
4421 if (nr > 0 && nc > 0) |
|
|
4422 { |
|
|
4423 if (nr == 1) |
|
|
4424 { |
|
|
4425 retval.resize (1, 1); |
|
|
4426 retval.elem (0, 0) = 1.0; |
|
|
4427 for (int j = 0; j < nc; j++) |
|
|
4428 { |
|
|
4429 if (elem (0, j) == 0.0) |
|
|
4430 { |
|
|
4431 retval.elem (0, 0) = 0.0; |
|
|
4432 break; |
|
|
4433 } |
|
|
4434 } |
|
|
4435 } |
|
|
4436 else if (nc == 1) |
|
|
4437 { |
|
|
4438 retval.resize (1, 1); |
|
|
4439 retval.elem (0, 0) = 1.0; |
|
|
4440 for (int i = 0; i < nr; i++) |
|
|
4441 { |
|
|
4442 if (elem (i, 0) == 0.0) |
|
|
4443 { |
|
|
4444 retval.elem (0, 0) = 0.0; |
|
|
4445 break; |
|
|
4446 } |
|
|
4447 } |
|
|
4448 } |
|
|
4449 else |
|
|
4450 { |
|
|
4451 retval.resize (1, nc); |
|
|
4452 for (int j = 0; j < nc; j++) |
|
|
4453 { |
|
|
4454 retval.elem (0, j) = 1.0; |
|
|
4455 for (int i = 0; i < nr; i++) |
|
|
4456 { |
|
|
4457 if (elem (i, j) == 0.0) |
|
|
4458 { |
|
|
4459 retval.elem (0, j) = 0.0; |
|
|
4460 break; |
|
|
4461 } |
|
|
4462 } |
|
|
4463 } |
|
|
4464 } |
|
|
4465 } |
|
|
4466 return retval; |
|
|
4467 } |
|
|
4468 |
|
|
4469 Matrix |
|
|
4470 ComplexMatrix::any (void) const |
|
|
4471 { |
|
|
4472 Matrix retval; |
|
|
4473 if (nr > 0 && nc > 0) |
|
|
4474 { |
|
|
4475 if (nr == 1) |
|
|
4476 { |
|
|
4477 retval.resize (1, 1); |
|
|
4478 retval.elem (0, 0) = 0.0; |
|
|
4479 for (int j = 0; j < nc; j++) |
|
|
4480 { |
|
|
4481 if (elem (0, j) != 0.0) |
|
|
4482 { |
|
|
4483 retval.elem (0, 0) = 1.0; |
|
|
4484 break; |
|
|
4485 } |
|
|
4486 } |
|
|
4487 } |
|
|
4488 else if (nc == 1) |
|
|
4489 { |
|
|
4490 retval.resize (1, 1); |
|
|
4491 retval.elem (0, 0) = 0.0; |
|
|
4492 for (int i = 0; i < nr; i++) |
|
|
4493 { |
|
|
4494 if (elem (i, 0) != 0.0) |
|
|
4495 { |
|
|
4496 retval.elem (0, 0) = 1.0; |
|
|
4497 break; |
|
|
4498 } |
|
|
4499 } |
|
|
4500 } |
|
|
4501 else |
|
|
4502 { |
|
|
4503 retval.resize (1, nc); |
|
|
4504 for (int j = 0; j < nc; j++) |
|
|
4505 { |
|
|
4506 retval.elem (0, j) = 0.0; |
|
|
4507 for (int i = 0; i < nr; i++) |
|
|
4508 { |
|
|
4509 if (elem (i, j) != 0.0) |
|
|
4510 { |
|
|
4511 retval.elem (0, j) = 1.0; |
|
|
4512 break; |
|
|
4513 } |
|
|
4514 } |
|
|
4515 } |
|
|
4516 } |
|
|
4517 } |
|
|
4518 return retval; |
|
|
4519 } |
|
|
4520 |
|
|
4521 ComplexMatrix |
|
|
4522 ComplexMatrix::cumprod (void) const |
|
|
4523 { |
|
|
4524 ComplexMatrix retval; |
|
|
4525 if (nr > 0 && nc > 0) |
|
|
4526 { |
|
|
4527 if (nr == 1) |
|
|
4528 { |
|
|
4529 retval.resize (1, nc); |
|
|
4530 Complex prod = elem (0, 0); |
|
|
4531 for (int j = 0; j < nc; j++) |
|
|
4532 { |
|
|
4533 retval.elem (0, j) = prod; |
|
|
4534 if (j < nc - 1) |
|
|
4535 prod *= elem (0, j+1); |
|
|
4536 } |
|
|
4537 } |
|
|
4538 else if (nc == 1) |
|
|
4539 { |
|
|
4540 retval.resize (nr, 1); |
|
|
4541 Complex prod = elem (0, 0); |
|
|
4542 for (int i = 0; i < nr; i++) |
|
|
4543 { |
|
|
4544 retval.elem (i, 0) = prod; |
|
|
4545 if (i < nr - 1) |
|
|
4546 prod *= elem (i+1, 0); |
|
|
4547 } |
|
|
4548 } |
|
|
4549 else |
|
|
4550 { |
|
|
4551 retval.resize (nr, nc); |
|
|
4552 for (int j = 0; j < nc; j++) |
|
|
4553 { |
|
|
4554 Complex prod = elem (0, j); |
|
|
4555 for (int i = 0; i < nr; i++) |
|
|
4556 { |
|
|
4557 retval.elem (i, j) = prod; |
|
|
4558 if (i < nr - 1) |
|
|
4559 prod *= elem (i+1, j); |
|
|
4560 } |
|
|
4561 } |
|
|
4562 } |
|
|
4563 } |
|
|
4564 return retval; |
|
|
4565 } |
|
|
4566 |
|
|
4567 ComplexMatrix |
|
|
4568 ComplexMatrix::cumsum (void) const |
|
|
4569 { |
|
|
4570 ComplexMatrix retval; |
|
|
4571 if (nr > 0 && nc > 0) |
|
|
4572 { |
|
|
4573 if (nr == 1) |
|
|
4574 { |
|
|
4575 retval.resize (1, nc); |
|
|
4576 Complex sum = elem (0, 0); |
|
|
4577 for (int j = 0; j < nc; j++) |
|
|
4578 { |
|
|
4579 retval.elem (0, j) = sum; |
|
|
4580 if (j < nc - 1) |
|
|
4581 sum += elem (0, j+1); |
|
|
4582 } |
|
|
4583 } |
|
|
4584 else if (nc == 1) |
|
|
4585 { |
|
|
4586 retval.resize (nr, 1); |
|
|
4587 Complex sum = elem (0, 0); |
|
|
4588 for (int i = 0; i < nr; i++) |
|
|
4589 { |
|
|
4590 retval.elem (i, 0) = sum; |
|
|
4591 if (i < nr - 1) |
|
|
4592 sum += elem (i+1, 0); |
|
|
4593 } |
|
|
4594 } |
|
|
4595 else |
|
|
4596 { |
|
|
4597 retval.resize (nr, nc); |
|
|
4598 for (int j = 0; j < nc; j++) |
|
|
4599 { |
|
|
4600 Complex sum = elem (0, j); |
|
|
4601 for (int i = 0; i < nr; i++) |
|
|
4602 { |
|
|
4603 retval.elem (i, j) = sum; |
|
|
4604 if (i < nr - 1) |
|
|
4605 sum += elem (i+1, j); |
|
|
4606 } |
|
|
4607 } |
|
|
4608 } |
|
|
4609 } |
|
|
4610 return retval; |
|
|
4611 } |
|
|
4612 |
|
|
4613 ComplexMatrix |
|
|
4614 ComplexMatrix::prod (void) const |
|
|
4615 { |
|
|
4616 ComplexMatrix retval; |
|
|
4617 if (nr > 0 && nc > 0) |
|
|
4618 { |
|
|
4619 if (nr == 1) |
|
|
4620 { |
|
|
4621 retval.resize (1, 1); |
|
|
4622 retval.elem (0, 0) = 1.0; |
|
|
4623 for (int j = 0; j < nc; j++) |
|
|
4624 retval.elem (0, 0) *= elem (0, j); |
|
|
4625 } |
|
|
4626 else if (nc == 1) |
|
|
4627 { |
|
|
4628 retval.resize (1, 1); |
|
|
4629 retval.elem (0, 0) = 1.0; |
|
|
4630 for (int i = 0; i < nr; i++) |
|
|
4631 retval.elem (0, 0) *= elem (i, 0); |
|
|
4632 } |
|
|
4633 else |
|
|
4634 { |
|
|
4635 retval.resize (1, nc); |
|
|
4636 for (int j = 0; j < nc; j++) |
|
|
4637 { |
|
|
4638 retval.elem (0, j) = 1.0; |
|
|
4639 for (int i = 0; i < nr; i++) |
|
|
4640 retval.elem (0, j) *= elem (i, j); |
|
|
4641 } |
|
|
4642 } |
|
|
4643 } |
|
|
4644 return retval; |
|
|
4645 } |
|
|
4646 |
|
|
4647 ComplexMatrix |
|
|
4648 ComplexMatrix::sum (void) const |
|
|
4649 { |
|
|
4650 ComplexMatrix retval; |
|
|
4651 if (nr > 0 && nc > 0) |
|
|
4652 { |
|
|
4653 if (nr == 1) |
|
|
4654 { |
|
|
4655 retval.resize (1, 1); |
|
|
4656 retval.elem (0, 0) = 0.0; |
|
|
4657 for (int j = 0; j < nc; j++) |
|
|
4658 retval.elem (0, 0) += elem (0, j); |
|
|
4659 } |
|
|
4660 else if (nc == 1) |
|
|
4661 { |
|
|
4662 retval.resize (1, 1); |
|
|
4663 retval.elem (0, 0) = 0.0; |
|
|
4664 for (int i = 0; i < nr; i++) |
|
|
4665 retval.elem (0, 0) += elem (i, 0); |
|
|
4666 } |
|
|
4667 else |
|
|
4668 { |
|
|
4669 retval.resize (1, nc); |
|
|
4670 for (int j = 0; j < nc; j++) |
|
|
4671 { |
|
|
4672 retval.elem (0, j) = 0.0; |
|
|
4673 for (int i = 0; i < nr; i++) |
|
|
4674 retval.elem (0, j) += elem (i, j); |
|
|
4675 } |
|
|
4676 } |
|
|
4677 } |
|
|
4678 return retval; |
|
|
4679 } |
|
|
4680 |
|
|
4681 ComplexMatrix |
|
|
4682 ComplexMatrix::sumsq (void) const |
|
|
4683 { |
|
|
4684 ComplexMatrix retval; |
|
|
4685 if (nr > 0 && nc > 0) |
|
|
4686 { |
|
|
4687 if (nr == 1) |
|
|
4688 { |
|
|
4689 retval.resize (1, 1); |
|
|
4690 retval.elem (0, 0) = 0.0; |
|
|
4691 for (int j = 0; j < nc; j++) |
|
|
4692 { |
|
|
4693 Complex d = elem (0, j); |
|
|
4694 retval.elem (0, 0) += d * d; |
|
|
4695 } |
|
|
4696 } |
|
|
4697 else if (nc == 1) |
|
|
4698 { |
|
|
4699 retval.resize (1, 1); |
|
|
4700 retval.elem (0, 0) = 0.0; |
|
|
4701 for (int i = 0; i < nr; i++) |
|
|
4702 { |
|
|
4703 Complex d = elem (i, 0); |
|
|
4704 retval.elem (0, 0) += d * d; |
|
|
4705 } |
|
|
4706 } |
|
|
4707 else |
|
|
4708 { |
|
|
4709 retval.resize (1, nc); |
|
|
4710 for (int j = 0; j < nc; j++) |
|
|
4711 { |
|
|
4712 retval.elem (0, j) = 0.0; |
|
|
4713 for (int i = 0; i < nr; i++) |
|
|
4714 { |
|
|
4715 Complex d = elem (i, j); |
|
|
4716 retval.elem (0, j) += d * d; |
|
|
4717 } |
|
|
4718 } |
|
|
4719 } |
|
|
4720 } |
|
|
4721 return retval; |
|
|
4722 } |
|
|
4723 |
|
|
4724 ComplexColumnVector |
|
|
4725 ComplexMatrix::diag (void) const |
|
|
4726 { |
|
|
4727 return diag (0); |
|
|
4728 } |
|
|
4729 |
|
|
4730 ComplexColumnVector |
|
|
4731 ComplexMatrix::diag (int k) const |
|
|
4732 { |
|
|
4733 int nnr = nr; |
|
|
4734 int nnc = nc; |
|
|
4735 if (k > 0) |
|
|
4736 nnc -= k; |
|
|
4737 else if (k < 0) |
|
|
4738 nnr += k; |
|
|
4739 |
|
|
4740 ComplexColumnVector d; |
|
|
4741 |
|
|
4742 if (nnr > 0 && nnc > 0) |
|
|
4743 { |
|
|
4744 int ndiag = (nnr < nnc) ? nnr : nnc; |
|
|
4745 |
|
|
4746 d.resize (ndiag); |
|
|
4747 |
|
|
4748 if (k > 0) |
|
|
4749 { |
|
|
4750 for (int i = 0; i < ndiag; i++) |
|
|
4751 d.elem (i) = elem (i, i+k); |
|
|
4752 } |
|
|
4753 else if ( k < 0) |
|
|
4754 { |
|
|
4755 for (int i = 0; i < ndiag; i++) |
|
|
4756 d.elem (i) = elem (i-k, i); |
|
|
4757 } |
|
|
4758 else |
|
|
4759 { |
|
|
4760 for (int i = 0; i < ndiag; i++) |
|
|
4761 d.elem (i) = elem (i, i); |
|
|
4762 } |
|
|
4763 } |
|
|
4764 else |
|
|
4765 cerr << "diag: requested diagonal out of range\n"; |
|
|
4766 |
|
|
4767 return d; |
|
|
4768 } |
|
|
4769 |
|
|
4770 ComplexColumnVector |
|
|
4771 ComplexMatrix::row_min (void) const |
|
|
4772 { |
|
|
4773 ComplexColumnVector result; |
|
|
4774 |
|
|
4775 if (nr > 0 && nc > 0) |
|
|
4776 { |
|
|
4777 result.resize (nr); |
|
|
4778 |
|
|
4779 for (int i = 0; i < nr; i++) |
|
|
4780 { |
|
|
4781 Complex res = elem (i, 0); |
|
|
4782 double absres = abs (res); |
|
|
4783 for (int j = 1; j < nc; j++) |
|
|
4784 if (abs (elem (i, j)) < absres) |
|
|
4785 { |
|
|
4786 res = elem (i, j); |
|
|
4787 absres = abs (res); |
|
|
4788 } |
|
|
4789 result.elem (i) = res; |
|
|
4790 } |
|
|
4791 } |
|
|
4792 |
|
|
4793 return result; |
|
|
4794 } |
|
|
4795 |
|
|
4796 ComplexColumnVector |
|
210
|
4797 ComplexMatrix::row_min_loc (void) const |
|
|
4798 { |
|
|
4799 ComplexColumnVector result; |
|
|
4800 |
|
|
4801 if (nr > 0 && nc > 0) |
|
|
4802 { |
|
|
4803 result.resize (nr); |
|
|
4804 |
|
|
4805 for (int i = 0; i < nr; i++) |
|
|
4806 { |
|
|
4807 Complex res = 0; |
|
|
4808 double absres = abs (elem (i, 0)); |
|
|
4809 for (int j = 0; j < nc; j++) |
|
|
4810 if (abs (elem (i, j)) < absres) |
|
|
4811 { |
|
|
4812 res = j; |
|
|
4813 absres = abs (elem (i, j)); |
|
|
4814 } |
|
|
4815 result.elem (i) = res + 1; |
|
|
4816 } |
|
|
4817 } |
|
|
4818 |
|
|
4819 return result; |
|
|
4820 } |
|
|
4821 |
|
|
4822 ComplexColumnVector |
|
3
|
4823 ComplexMatrix::row_max (void) const |
|
|
4824 { |
|
|
4825 ComplexColumnVector result; |
|
|
4826 |
|
|
4827 if (nr > 0 && nc > 0) |
|
|
4828 { |
|
|
4829 result.resize (nr); |
|
|
4830 |
|
|
4831 for (int i = 0; i < nr; i++) |
|
|
4832 { |
|
|
4833 Complex res = elem (i, 0); |
|
|
4834 double absres = abs (res); |
|
|
4835 for (int j = 1; j < nc; j++) |
|
|
4836 if (abs (elem (i, j)) > absres) |
|
|
4837 { |
|
|
4838 res = elem (i, j); |
|
|
4839 absres = abs (res); |
|
|
4840 } |
|
|
4841 result.elem (i) = res; |
|
|
4842 } |
|
|
4843 } |
|
|
4844 |
|
|
4845 return result; |
|
|
4846 } |
|
|
4847 |
|
210
|
4848 ComplexColumnVector |
|
|
4849 ComplexMatrix::row_max_loc (void) const |
|
|
4850 { |
|
|
4851 ComplexColumnVector result; |
|
|
4852 |
|
|
4853 if (nr > 0 && nc > 0) |
|
|
4854 { |
|
|
4855 result.resize (nr); |
|
|
4856 |
|
|
4857 for (int i = 0; i < nr; i++) |
|
|
4858 { |
|
|
4859 Complex res = 0; |
|
|
4860 double absres = abs (elem (i, 0)); |
|
|
4861 for (int j = 0; j < nc; j++) |
|
|
4862 if (abs (elem (i, j)) > absres) |
|
|
4863 { |
|
|
4864 res = j; |
|
|
4865 absres = abs (elem (i, j)); |
|
|
4866 } |
|
|
4867 result.elem (i) = res + 1; |
|
|
4868 } |
|
|
4869 } |
|
|
4870 |
|
|
4871 return result; |
|
|
4872 } |
|
|
4873 |
|
3
|
4874 ComplexRowVector |
|
|
4875 ComplexMatrix::column_min (void) const |
|
|
4876 { |
|
|
4877 ComplexRowVector result; |
|
|
4878 |
|
|
4879 if (nr > 0 && nc > 0) |
|
|
4880 { |
|
|
4881 result.resize (nc); |
|
|
4882 |
|
|
4883 for (int j = 0; j < nc; j++) |
|
|
4884 { |
|
|
4885 Complex res = elem (0, j); |
|
|
4886 double absres = abs (res); |
|
|
4887 for (int i = 1; i < nr; i++) |
|
|
4888 if (abs (elem (i, j)) < absres) |
|
|
4889 { |
|
|
4890 res = elem (i, j); |
|
|
4891 absres = abs (res); |
|
|
4892 } |
|
|
4893 result.elem (j) = res; |
|
|
4894 } |
|
|
4895 } |
|
|
4896 |
|
|
4897 return result; |
|
|
4898 } |
|
|
4899 |
|
|
4900 ComplexRowVector |
|
210
|
4901 ComplexMatrix::column_min_loc (void) const |
|
|
4902 { |
|
|
4903 ComplexRowVector result; |
|
|
4904 |
|
|
4905 if (nr > 0 && nc > 0) |
|
|
4906 { |
|
|
4907 result.resize (nc); |
|
|
4908 |
|
|
4909 for (int j = 0; j < nc; j++) |
|
|
4910 { |
|
|
4911 Complex res = 0; |
|
|
4912 double absres = abs (elem (0, j)); |
|
|
4913 for (int i = 0; i < nr; i++) |
|
|
4914 if (abs (elem (i, j)) < absres) |
|
|
4915 { |
|
|
4916 res = i; |
|
|
4917 absres = abs (elem (i, j)); |
|
|
4918 } |
|
|
4919 result.elem (j) = res + 1; |
|
|
4920 } |
|
|
4921 } |
|
|
4922 |
|
|
4923 return result; |
|
|
4924 } |
|
|
4925 |
|
|
4926 ComplexRowVector |
|
3
|
4927 ComplexMatrix::column_max (void) const |
|
|
4928 { |
|
|
4929 ComplexRowVector result; |
|
|
4930 |
|
|
4931 if (nr > 0 && nc > 0) |
|
|
4932 { |
|
|
4933 result.resize (nc); |
|
|
4934 |
|
|
4935 for (int j = 0; j < nc; j++) |
|
|
4936 { |
|
|
4937 Complex res = elem (0, j); |
|
|
4938 double absres = abs (res); |
|
|
4939 for (int i = 1; i < nr; i++) |
|
|
4940 if (abs (elem (i, j)) > absres) |
|
|
4941 { |
|
|
4942 res = elem (i, j); |
|
|
4943 absres = abs (res); |
|
|
4944 } |
|
|
4945 result.elem (j) = res; |
|
|
4946 } |
|
|
4947 } |
|
|
4948 |
|
|
4949 return result; |
|
|
4950 } |
|
|
4951 |
|
210
|
4952 ComplexRowVector |
|
|
4953 ComplexMatrix::column_max_loc (void) const |
|
|
4954 { |
|
|
4955 ComplexRowVector result; |
|
|
4956 |
|
|
4957 if (nr > 0 && nc > 0) |
|
|
4958 { |
|
|
4959 result.resize (nc); |
|
|
4960 |
|
|
4961 for (int j = 0; j < nc; j++) |
|
|
4962 { |
|
|
4963 Complex res = 0; |
|
|
4964 double absres = abs (elem (0, j)); |
|
|
4965 for (int i = 0; i < nr; i++) |
|
|
4966 if (abs (elem (i, j)) > absres) |
|
|
4967 { |
|
|
4968 res = i; |
|
|
4969 absres = abs (elem (i, j)); |
|
|
4970 } |
|
|
4971 result.elem (j) = res + 1; |
|
|
4972 } |
|
|
4973 } |
|
|
4974 |
|
|
4975 return result; |
|
|
4976 } |
|
|
4977 |
|
3
|
4978 // i/o |
|
|
4979 |
|
|
4980 ostream& |
|
|
4981 operator << (ostream& os, const ComplexMatrix& a) |
|
|
4982 { |
|
|
4983 // int field_width = os.precision () + 7; |
|
|
4984 for (int i = 0; i < a.nr; i++) |
|
|
4985 { |
|
|
4986 for (int j = 0; j < a.nc; j++) |
|
|
4987 os << " " /* setw (field_width) */ << a.elem (i, j); |
|
|
4988 os << "\n"; |
|
|
4989 } |
|
|
4990 return os; |
|
|
4991 } |
|
|
4992 |
|
|
4993 istream& |
|
|
4994 operator >> (istream& is, ComplexMatrix& a) |
|
|
4995 { |
|
|
4996 int nr = a.rows (); |
|
|
4997 int nc = a.columns (); |
|
|
4998 |
|
|
4999 if (nr < 1 || nc < 1) |
|
|
5000 is.clear (ios::badbit); |
|
|
5001 else |
|
|
5002 { |
|
|
5003 Complex tmp; |
|
|
5004 for (int i = 0; i < nr; i++) |
|
|
5005 for (int j = 0; j < nc; j++) |
|
|
5006 { |
|
|
5007 is >> tmp; |
|
|
5008 if (is) |
|
|
5009 a.elem (i, j) = tmp; |
|
|
5010 else |
|
|
5011 break; |
|
|
5012 } |
|
|
5013 } |
|
|
5014 |
|
|
5015 return is; |
|
|
5016 } |
|
|
5017 |
|
|
5018 /* |
|
|
5019 ;;; Local Variables: *** |
|
|
5020 ;;; mode: C++ *** |
|
|
5021 ;;; page-delimiter: "^/\\*" *** |
|
|
5022 ;;; End: *** |
|
|
5023 */ |
