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