Mercurial > hg > octave-lyh

comparison liboctave/Matrix.h @ 3:9a4c07481e61

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-08 01:20:23 by jwe] Initial revision
author jwe
date Sun, 08 Aug 1993 01:21:46 +0000
parents
children 1b59f5c6c370
comparison
equal deleted inserted replaced
2:c0190df9885d 3:9a4c07481e61
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 /*
25
26 Should probably say something here about why these classes are not
27 represented by some sort of inheritance tree...
28
29 */
30
31 #if !defined (_Matrix_h)
32 #define _Matrix_h 1
33
34 // I\'m not sure how this is supposed to work if the .h file declares
35 // several classes, each of which is defined in a separate file...
36 //
37 // #ifdef __GNUG__
38 // #pragma interface
39 // #endif
40
41 #include <stdlib.h>
42 #include <stddef.h>
43 #include <math.h>
44 #include <values.h>
45 #include <assert.h>
46 #include <iostream.h>
47 // #include <iomanip.h> // We don\'t use this yet.
48 #include <Complex.h>
49
50 #define FAIL assert(0) /* XXX FIXME XXX */
51
52 #ifndef MAPPER_FCN_TYPEDEFS
53 #define MAPPER_FCN_TYPEDEFS 1
54
55 typedef double (*d_d_Mapper)(double);
56 typedef double (*d_c_Mapper)(const Complex&);
57 typedef Complex (*c_c_Mapper)(const Complex&);
58
59 #endif
60
61 #include "f77-uscore.h"
62
63 // Fortran functions we call.
64
65 extern "C"
66 {
67 int F77_FCN (dgemm) (const char*, const char*, const int*,
68 const int*, const int*, const double*,
69 const double*, const int*, const double*,
70 const int*, const double*, double*, const int*,
71 long, long);
72
73 int F77_FCN (dgemv) (const char*, const int*, const int*,
74 const double*, const double*, const int*,
75 const double*, const int*, const double*,
76 double*, const int*, long);
77
78 int F77_FCN (dgeco) (double*, const int*, const int*, int*, double*,
79 double*);
80
81 int F77_FCN (dgesv) (const int*, const int*, double*, const int*,
82 int*, double*, const int*, int*);
83
84 int F77_FCN (dgeqrf) (const int*, const int*, double*, const int*,
85 double*, double*, const int*, int*);
86
87 int F77_FCN (dorgqr) (const int*, const int*, const int*, double*,
88 const int*, double*, double*, const int*, int*);
89
90 int F77_FCN (dgesl) (const double*, const int*, const int*,
91 const int*, double*, const int*);
92
93 int F77_FCN (dgedi) (double*, const int*, const int*, const int*,
94 double*, double*, const int*);
95
96 double F77_FCN (ddot) (const int*, const double*, const int*,
97 const double*, const int*);
98
99 int F77_FCN (dgeev) (const char*, const char*, const int*, double*,
100 const int*, double*, double*, double*,
101 const int*, double*, const int*, double*,
102 const int*, int*, long, long);
103
104 int F77_FCN (dgeesx) (const char*, const char*, int (*)(), const char*,
105 const int*, double*, const int*, int*, double*,
106 double*, double*, const int*, double*, double*,
107 double*, const int*, int*, const int*, int*,
108 int*, long, long);
109
110 int F77_FCN (dhseqr) (const char*, const char*, const int*,
111 const int*, const int*, double*,
112 const int*, double*, double*,
113 double*, const int*, double*, const int*,
114 int*, long, long);
115
116 int F77_FCN (dgebal) (const char*, const int*, double*,
117 const int*, int*, int*, double*,
118 int*, long, long);
119
120 int F77_FCN (dgebak) (const char*, const char*, const int*, const int*,
121 const int*, double*, const int*, double*, const int*,
122 int*, long, long);
123
124 int F77_FCN (dgehrd) (const int*, const int*, const int*,
125 double*, const int*, double*, double*,
126 const int*, int*, long, long);
127
128 int F77_FCN (dorghr) (const int*, const int*, const int*,
129 double*, const int*, double*, double*,
130 const int*, int*, long, long);
131
132 int F77_FCN (dgesvd) (const char*, const char*, const int*,
133 const int*, double*, const int*, double*,
134 double*, const int*, double*, const int*,
135 double*, const int*, int*, long, long);
136
137 int F77_FCN (dgelss) (const int*, const int*, const int*, double*,
138 const int*, double*, const int*, double*,
139 const double*, int*, double*, const int*,
140 int*);
141
142 /*
143 * f2c translates complex*16 as
144 *
145 * typedef struct { doublereal re, im; } doublecomplex;
146 *
147 * and Complex.h from libg++ uses
148 *
149 * protected:
150 * double re;
151 * double im;
152 *
153 * as the only data members, so this should work (fingers crossed that
154 * things don't change).
155 */
156
157 int F77_FCN (zgemm) (const char*, const char*, const int*,
158 const int*, const int*, const Complex*,
159 const Complex*, const int*, const Complex*,
160 const int*, const Complex*, Complex*, const int*,
161 long, long);
162
163 int F77_FCN (zgemv) (const char*, const int*, const int*,
164 const Complex*, const Complex*, const int*,
165 const Complex*, const int*, const Complex*,
166 Complex*, const int*, long);
167
168 int F77_FCN (zgeco) (Complex*, const int*, const int*, int*,
169 double*, Complex*);
170
171 int F77_FCN (zgesv) (const int*, const int*, Complex*, const int*,
172 int*, Complex*, const int*, int*);
173
174 int F77_FCN (zgeqrf) (const int*, const int*, Complex*, const int*,
175 Complex*, Complex*, const int*, int*);
176
177 int F77_FCN (zgeesx) (const char*, const char*, int (*)(), const char*,
178 const int*, Complex*, const int*, int*,
179 Complex*, Complex*, const int*, double*, double*,
180 Complex*, const int*, double*, int*, int*,
181 long, long);
182
183 int F77_FCN (zhseqr) (const char*, const char*, const int*,
184 const int*, const int*, Complex*, const int*,
185 Complex*, Complex*, const int*, Complex*,
186 const int*, int*, long, long);
187
188 int F77_FCN (zgebal) (const char*, const int*, Complex*, const int*,
189 int*, int*, double*, int*, long, long);
190
191 int F77_FCN (zgebak) (const char*, const char*, const int*, const int*,
192 const int*, double*, const int*, Complex*,
193 const int*, int*, long, long);
194
195 int F77_FCN (zgehrd) (const int*, const int*, const int*, Complex*,
196 const int*, Complex*, Complex*, const int*,
197 int*, long, long);
198
199 int F77_FCN (zunghr) (const int*, const int*, const int*, Complex*,
200 const int*, Complex*, Complex*, const int*,
201 int*, long, long);
202
203 int F77_FCN (zungqr) (const int*, const int*, const int*, Complex*,
204 const int*, Complex*, Complex*, const int*, int*);
205
206 int F77_FCN (zgedi) (Complex*, const int*, const int*, int*,
207 Complex*, Complex*, const int*);
208
209 int F77_FCN (zgesl) (Complex*, const int*, const int*, int*,
210 Complex*, const int*);
211
212 int F77_FCN (zgeev) (const char*, const char*, const int*, Complex*,
213 const int*, Complex*, Complex*, const int*,
214 Complex*, const int*, Complex*, const int*,
215 double*, int*, long, long);
216
217 int F77_FCN (zgesvd) (const char*, const char*, const int*,
218 const int*, Complex*, const int*, double*,
219 Complex*, const int*, Complex*, const int*,
220 Complex*, const int*, double*, int*, long, long);
221
222 int F77_FCN (zgelss) (const int*, const int*, const int*, Complex*,
223 const int*, Complex*, const int*, double*,
224 const double*, int*, Complex*, const int*,
225 double*, int*);
226
227 // Note that the original complex fft routines were not written for
228 // double complex arguments. They have been modified by adding an
229 // implicit double precision (a-h,o-z) statement at the beginning of
230 // each subroutine.
231
232 int F77_FCN (cffti) (const int*, Complex*);
233
234 int F77_FCN (cfftf) (const int*, Complex*, Complex*);
235
236 int F77_FCN (cfftb) (const int*, Complex*, Complex*);
237
238 }
239
240 // Classes we declare.
241
242 class Matrix;
243 class ColumnVector;
244 class RowVector;
245 class DiagMatrix;
246 class ComplexMatrix;
247 class ComplexColumnVector;
248 class ComplexRowVector;
249 class ComplexDiagMatrix;
250 class DET;
251 class ComplexDET;
252 class EIG;
253 class HESS;
254 class ComplexHESS;
255 class SCHUR;
256 class ComplexSCHUR;
257 class SVD;
258 class ComplexSVD;
259 class LU;
260 class ComplexLU;
261 class QR;
262 class ComplexQR;
263
264 /*
265 * Matrix class
266 */
267
268 class Matrix
269 {
270 friend class RowVector;
271 friend class DiagMatrix;
272 friend class ComplexMatrix;
273 friend class ComplexDiagMatrix;
274 friend class EIG;
275 friend class HESS;
276 friend class SCHUR;
277 friend class SVD;
278 friend class LU;
279 friend class QR;
280
281 public:
282 Matrix (void);
283 Matrix (int r, int c);
284 Matrix (int r, int c, double val);
285 Matrix (const Matrix& a);
286 Matrix (const DiagMatrix& a);
287 Matrix (double a);
288 ~Matrix (void);
289
290 #if defined (MDEBUG)
291 void *operator new (size_t size)
292 {
293 Matrix *p = ::new Matrix;
294 cerr << "Matrix::new(): " << p << "\n";
295 return p;
296 }
297
298 void operator delete (void *p, size_t size)
299 {
300 cerr << "Matrix::delete(): " << p << "\n";
301 ::delete p;
302 }
303 #endif
304
305 Matrix& operator = (const Matrix& a);
306
307 int rows (void) const;
308 int cols (void) const;
309 int columns (void) const;
310
311 double& elem (int r, int c);
312 double& checkelem (int r, int c);
313 double& operator () (int r, int c);
314
315 double elem (int r, int c) const; // const access
316 double checkelem (int r, int c) const;
317 double operator () (int r, int c) const;
318
319 Matrix& resize (int r, int c);
320 Matrix& resize (int r, int c, double val);
321
322 int operator == (const Matrix& a) const;
323 int operator != (const Matrix& a) const;
324
325 // destructive insert/delete/reorder operations
326
327 Matrix& insert (const Matrix& a, int r, int c);
328 Matrix& insert (const RowVector& a, int r, int c);
329 Matrix& insert (const ColumnVector& a, int r, int c);
330 Matrix& insert (const DiagMatrix& a, int r, int c);
331
332 Matrix& fill (double val);
333 Matrix& fill (double val, int r1, int c1, int r2, int c2);
334
335 Matrix append (const Matrix& a) const;
336 Matrix append (const RowVector& a) const;
337 Matrix append (const ColumnVector& a) const;
338 Matrix append (const DiagMatrix& a) const;
339
340 Matrix stack (const Matrix& a) const;
341 Matrix stack (const RowVector& a) const;
342 Matrix stack (const ColumnVector& a) const;
343 Matrix stack (const DiagMatrix& a) const;
344
345 Matrix transpose (void) const;
346
347 // resize is the destructive equivalent for this one
348
349 Matrix extract (int r1, int c1, int r2, int c2) const;
350
351 // extract row or column i.
352
353 RowVector row (int i) const;
354 RowVector row (char *s) const;
355
356 ColumnVector column (int i) const;
357 ColumnVector column (char *s) const;
358
359 Matrix inverse (int& info, double& rcond) const;
360 Matrix inverse (int& info) const;
361 Matrix inverse (void) const;
362
363 ComplexMatrix fourier (void) const;
364 ComplexMatrix ifourier (void) const;
365
366 DET determinant (void) const;
367 DET determinant (int& info) const;
368 DET determinant (int& info, double& rcond) const;
369
370 Matrix solve (const Matrix& b) const;
371 Matrix solve (const Matrix& b, int& info) const;
372 Matrix solve (const Matrix& b, int& info, double& rcond) const;
373
374 ComplexMatrix solve (const ComplexMatrix& b) const;
375 ComplexMatrix solve (const ComplexMatrix& b, int& info) const;
376 ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const;
377
378 ColumnVector solve (const ColumnVector& b) const;
379 ColumnVector solve (const ColumnVector& b, int& info) const;
380 ColumnVector solve (const ColumnVector& b, int& info, double& rcond) const;
381
382 ComplexColumnVector solve (const ComplexColumnVector& b) const;
383 ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const;
384 ComplexColumnVector solve (const ComplexColumnVector& b, int& info,
385 double& rcond) const;
386
387 Matrix lssolve (const Matrix& b) const;
388 Matrix lssolve (const Matrix& b, int& info) const;
389 Matrix lssolve (const Matrix& b, int& info, int& rank) const;
390
391 ComplexMatrix lssolve (const ComplexMatrix& b) const;
392 ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const;
393 ComplexMatrix lssolve (const ComplexMatrix& b, int& info,
394 int& rank) const;
395
396 ColumnVector lssolve (const ColumnVector& b) const;
397 ColumnVector lssolve (const ColumnVector& b, int& info) const;
398 ColumnVector lssolve (const ColumnVector& b, int& info, int& rank) const;
399
400 ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
401 ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const;
402 ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
403 int& rank) const;
404
405 // matrix by scalar -> matrix operations
406
407 Matrix operator + (double s) const;
408 Matrix operator - (double s) const;
409 Matrix operator * (double s) const;
410 Matrix operator / (double s) const;
411
412 ComplexMatrix operator + (Complex s) const;
413 ComplexMatrix operator - (Complex s) const;
414 ComplexMatrix operator * (Complex s) const;
415 ComplexMatrix operator / (Complex s) const;
416
417 // scalar by matrix -> matrix operations
418
419 friend Matrix operator + (double s, const Matrix& a);
420 friend Matrix operator - (double s, const Matrix& a);
421 friend Matrix operator * (double s, const Matrix& a);
422 friend Matrix operator / (double s, const Matrix& a);
423
424 // matrix by column vector -> column vector operations
425
426 ColumnVector operator * (const ColumnVector& a) const;
427
428 ComplexColumnVector operator * (const ComplexColumnVector& a) const;
429
430 // matrix by diagonal matrix -> matrix operations
431
432 Matrix operator + (const DiagMatrix& a) const;
433 Matrix operator - (const DiagMatrix& a) const;
434 Matrix operator * (const DiagMatrix& a) const;
435
436 ComplexMatrix operator + (const ComplexDiagMatrix& a) const;
437 ComplexMatrix operator - (const ComplexDiagMatrix& a) const;
438 ComplexMatrix operator * (const ComplexDiagMatrix& a) const;
439
440 Matrix& operator += (const DiagMatrix& a);
441 Matrix& operator -= (const DiagMatrix& a);
442
443 // matrix by matrix -> matrix operations
444
445 Matrix operator + (const Matrix& a) const;
446 Matrix operator - (const Matrix& a) const;
447 Matrix operator * (const Matrix& a) const;
448
449 ComplexMatrix operator + (const ComplexMatrix& a) const;
450 ComplexMatrix operator - (const ComplexMatrix& a) const;
451 ComplexMatrix operator * (const ComplexMatrix& a) const;
452
453 Matrix product (const Matrix& a) const; // element by element
454 Matrix quotient (const Matrix& a) const; // element by element
455
456 ComplexMatrix product (const ComplexMatrix& a) const; // element by element
457 ComplexMatrix quotient (const ComplexMatrix& a) const; // element by element
458
459 Matrix& operator += (const Matrix& a);
460 Matrix& operator -= (const Matrix& a);
461
462 // unary operations
463
464 Matrix operator - (void) const;
465 Matrix operator ! (void) const;
466
467 // other operations
468
469 friend Matrix map (d_d_Mapper f, const Matrix& a);
470 void map (d_d_Mapper f);
471
472 Matrix all (void) const;
473 Matrix any (void) const;
474
475 Matrix cumprod (void) const;
476 Matrix cumsum (void) const;
477 Matrix prod (void) const;
478 Matrix sum (void) const;
479 Matrix sumsq (void) const;
480
481 ColumnVector diag (void) const;
482 ColumnVector diag (int k) const;
483
484 ColumnVector row_min (void) const;
485 ColumnVector row_max (void) const;
486
487 RowVector column_min (void) const;
488 RowVector column_max (void) const;
489
490 // i/o
491
492 friend ostream& operator << (ostream& os, const Matrix& a);
493 friend istream& operator >> (istream& is, Matrix& a);
494
495 // conversions
496
497 double *fortran_vec (void);
498
499 private:
500 int nr;
501 int nc;
502 int len;
503 double *data;
504
505 Matrix (double *d, int r, int c);
506 };
507
508 inline Matrix::Matrix (void) { nr = 0; nc = 0; len = 0; data = 0; }
509
510 inline Matrix::Matrix (double *d, int r, int c)
511 { nr = r; nc = c; len = nr*nc; data = d; }
512
513 inline Matrix::~Matrix (void) { delete [] data; data = 0; }
514
515 inline int Matrix::rows (void) const { return nr; }
516 inline int Matrix::cols (void) const { return nc; }
517 inline int Matrix::columns (void) const { return nc; }
518
519 inline double& Matrix::elem (int r, int c) { return data[nr*c+r]; }
520
521 inline double& Matrix::checkelem (int r, int c)
522 {
523 #ifndef NO_RANGE_CHECK
524 if (r < 0 || r >= nr || c < 0 || c >= nc)
525 FAIL;
526 #endif
527
528 return elem (r, c);
529 }
530
531 inline double& Matrix::operator () (int r, int c)
532 { return checkelem (r, c); }
533
534 inline double Matrix::elem (int r, int c) const { return data[nr*c+r]; }
535
536 inline double Matrix::checkelem (int r, int c) const
537 {
538 #ifndef NO_RANGE_CHECK
539 if (r < 0 || r >= nr || c < 0 || c >= nc)
540 FAIL;
541 #endif
542
543 return elem (r, c);
544 }
545
546 inline double Matrix::operator () (int r, int c) const
547 { return checkelem (r, c); }
548
549 inline double *Matrix::fortran_vec (void) { return data; }
550
551 /*
552 * Column Vector class
553 */
554
555 class ColumnVector
556 {
557 friend class Matrix;
558 friend class RowVector;
559 friend class DiagMatrix;
560 friend class ComplexMatrix;
561 friend class ComplexColumnVector;
562 friend class ComplexDiagMatrix;
563
564 public:
565 ColumnVector (void);
566 ColumnVector (int n);
567 ColumnVector (int n, double val);
568 ColumnVector (const ColumnVector& a);
569 ColumnVector (double a);
570 ~ColumnVector (void);
571
572 ColumnVector& operator = (const ColumnVector& a);
573
574 int capacity (void) const;
575 int length (void) const;
576
577 double& elem (int n);
578 double& checkelem (int n);
579 double& operator () (int n);
580
581 double elem (int n) const; // const access
582 double checkelem (int n) const;
583 double operator () (int n) const;
584
585 ColumnVector& resize (int n);
586 ColumnVector& resize (int n, double val);
587
588 int operator == (const ColumnVector& a) const;
589 int operator != (const ColumnVector& a) const;
590
591 // destructive insert/delete/reorder operations
592
593 ColumnVector& insert (const ColumnVector& a, int r);
594
595 ColumnVector& fill (double val);
596 ColumnVector& fill (double val, int r1, int r2);
597
598 ColumnVector stack (const ColumnVector& a) const;
599
600 RowVector transpose (void) const;
601
602 // resize is the destructive equivalent for this one
603
604 ColumnVector extract (int r1, int r2) const;
605
606 // column vector by scalar -> column vector operations
607
608 ColumnVector operator + (double s) const;
609 ColumnVector operator - (double s) const;
610 ColumnVector operator * (double s) const;
611 ColumnVector operator / (double s) const;
612
613 ComplexColumnVector operator + (Complex s) const;
614 ComplexColumnVector operator - (Complex s) const;
615 ComplexColumnVector operator * (Complex s) const;
616 ComplexColumnVector operator / (Complex s) const;
617
618 // scalar by column vector -> column vector operations
619
620 friend ColumnVector operator + (double s, const ColumnVector& a);
621 friend ColumnVector operator - (double s, const ColumnVector& a);
622 friend ColumnVector operator * (double s, const ColumnVector& a);
623 friend ColumnVector operator / (double s, const ColumnVector& a);
624
625 // column vector by row vector -> matrix operations
626
627 Matrix operator * (const RowVector& a) const;
628
629 ComplexMatrix operator * (const ComplexRowVector& a) const;
630
631 // column vector by column vector -> column vector operations
632
633 ColumnVector operator + (const ColumnVector& a) const;
634 ColumnVector operator - (const ColumnVector& a) const;
635
636 ComplexColumnVector operator + (const ComplexColumnVector& a) const;
637 ComplexColumnVector operator - (const ComplexColumnVector& a) const;
638
639 ColumnVector product (const ColumnVector& a) const; // element by element
640 ColumnVector quotient (const ColumnVector& a) const; // element by element
641
642 ComplexColumnVector product (const ComplexColumnVector& a) const;
643 ComplexColumnVector quotient (const ComplexColumnVector& a) const;
644
645 ColumnVector& operator += (const ColumnVector& a);
646 ColumnVector& operator -= (const ColumnVector& a);
647
648 // unary operations
649
650 ColumnVector operator - (void) const;
651
652 friend ColumnVector map (d_d_Mapper f, const ColumnVector& a);
653 void map (d_d_Mapper f);
654
655 double min (void) const;
656 double max (void) const;
657
658 // i/o
659
660 friend ostream& operator << (ostream& os, const ColumnVector& a);
661
662 // conversions
663
664 double *fortran_vec (void);
665
666 private:
667 int len;
668 double *data;
669
670 ColumnVector (double *d, int l);
671 };
672
673 inline ColumnVector::ColumnVector (void) { len = 0; data = 0; }
674 inline ColumnVector::ColumnVector (double *d, int l) { len = l; data = d; }
675 inline ColumnVector::~ColumnVector (void) { delete [] data; data = 0; }
676
677 inline int ColumnVector::capacity (void) const { return len; }
678 inline int ColumnVector::length (void) const { return len; }
679
680 inline double& ColumnVector::elem (int n) { return data[n]; }
681
682 inline double&
683 ColumnVector::checkelem (int n)
684 {
685 #ifndef NO_RANGE_CHECK
686 if (n < 0 || n >= len)
687 FAIL;
688 #endif
689
690 return elem (n);
691 }
692
693 inline double& ColumnVector::operator () (int n) { return checkelem (n); }
694
695 inline double ColumnVector::elem (int n) const { return data[n]; }
696
697 inline double
698 ColumnVector::checkelem (int n) const
699 {
700 #ifndef NO_RANGE_CHECK
701 if (n < 0 || n >= len)
702 FAIL;
703 #endif
704
705 return elem (n);
706 }
707
708 inline double ColumnVector::operator () (int n) const { return checkelem (n); }
709
710 inline double *ColumnVector::fortran_vec (void) { return data; }
711
712 /*
713 * Row Vector class
714 */
715
716 class RowVector
717 {
718 friend class Matrix;
719 friend class DiagMatrix;
720 friend class ColumnVector;
721 friend class ComplexMatrix;
722 friend class ComplexRowVector;
723 friend class ComplexDiagMatrix;
724
725 public:
726 RowVector (void);
727 RowVector (int n);
728 RowVector (int n, double val);
729 RowVector (const RowVector& a);
730 RowVector (double a);
731 ~RowVector (void);
732
733 RowVector& operator = (const RowVector& a);
734
735 int capacity (void) const;
736 int length (void) const;
737
738 double& elem (int n);
739 double& checkelem (int n);
740 double& operator () (int n);
741
742 double elem (int n) const; // const access
743 double checkelem (int n) const;
744 double operator () (int n) const;
745
746 RowVector& resize (int n);
747 RowVector& resize (int n, double val);
748
749 int operator == (const RowVector& a) const;
750 int operator != (const RowVector& a) const;
751
752 // destructive insert/delete/reorder operations
753
754 RowVector& insert (const RowVector& a, int c);
755
756 RowVector& fill (double val);
757 RowVector& fill (double val, int c1, int c2);
758
759 RowVector append (const RowVector& a) const;
760
761 ColumnVector transpose (void) const;
762
763 // resize is the destructive equivalent for this one
764
765 RowVector extract (int c1, int c2) const;
766
767 // row vector by scalar -> row vector operations
768
769 RowVector operator + (double s) const;
770 RowVector operator - (double s) const;
771 RowVector operator * (double s) const;
772 RowVector operator / (double s) const;
773
774 ComplexRowVector operator + (Complex s) const;
775 ComplexRowVector operator - (Complex s) const;
776 ComplexRowVector operator * (Complex s) const;
777 ComplexRowVector operator / (Complex s) const;
778
779 // scalar by row vector -> row vector operations
780
781 friend RowVector operator + (double s, const RowVector& a);
782 friend RowVector operator - (double s, const RowVector& a);
783 friend RowVector operator * (double s, const RowVector& a);
784 friend RowVector operator / (double s, const RowVector& a);
785
786 // row vector by column vector -> scalar
787
788 double operator * (const ColumnVector& a) const;
789
790 Complex operator * (const ComplexColumnVector& a) const;
791
792 // row vector by matrix -> row vector
793
794 RowVector operator * (const Matrix& a) const;
795
796 ComplexRowVector operator * (const ComplexMatrix& a) const;
797
798 // row vector by row vector -> row vector operations
799
800 RowVector operator + (const RowVector& a) const;
801 RowVector operator - (const RowVector& a) const;
802
803 ComplexRowVector operator + (const ComplexRowVector& a) const;
804 ComplexRowVector operator - (const ComplexRowVector& a) const;
805
806 RowVector product (const RowVector& a) const; // element by element
807 RowVector quotient (const RowVector& a) const; // element by element
808
809 ComplexRowVector product (const ComplexRowVector& a) const; // el by el
810 ComplexRowVector quotient (const ComplexRowVector& a) const; // el by el
811
812 RowVector& operator += (const RowVector& a);
813 RowVector& operator -= (const RowVector& a);
814
815 // unary operations
816
817 RowVector operator - (void) const;
818
819 friend RowVector map (d_d_Mapper f, const RowVector& a);
820 void map (d_d_Mapper f);
821
822 double min (void) const;
823 double max (void) const;
824
825 // i/o
826
827 friend ostream& operator << (ostream& os, const RowVector& a);
828
829 // conversions
830
831 double *fortran_vec (void);
832
833 private:
834 int len;
835 double *data;
836
837 RowVector (double *d, int l);
838 };
839
840 inline RowVector::RowVector (void) { len = 0; data = 0; }
841 inline RowVector::RowVector (double *d, int l) { len = l; data = d; }
842 inline RowVector::~RowVector (void) { delete [] data; data = 0; }
843
844 inline int RowVector::capacity (void) const { return len; }
845 inline int RowVector::length (void) const { return len; }
846
847 inline double& RowVector::elem (int n) { return data[n]; }
848
849 inline double&
850 RowVector::checkelem (int n)
851 {
852 #ifndef NO_RANGE_CHECK
853 if (n < 0 || n >= len)
854 FAIL;
855 #endif
856
857 return elem (n);
858 }
859
860 inline double& RowVector::operator () (int n) { return checkelem (n); }
861
862 inline double RowVector::elem (int n) const { return data[n]; }
863
864 inline double
865 RowVector::checkelem (int n) const
866 {
867 #ifndef NO_RANGE_CHECK
868 if (n < 0 || n >= len)
869 FAIL;
870 #endif
871
872 return elem (n);
873 }
874
875 inline double RowVector::operator () (int n) const { return checkelem (n); }
876
877 inline double *RowVector::fortran_vec (void) { return data; }
878
879 /*
880 * Diagonal Matrix class
881 */
882
883 class DiagMatrix
884 {
885 friend class Matrix;
886 friend class ComplexMatrix;
887 friend class ComplexDiagMatrix;
888
889 public:
890 DiagMatrix (void);
891 DiagMatrix (int n);
892 DiagMatrix (int n, double val);
893 DiagMatrix (int r, int c);
894 DiagMatrix (int r, int c, double val);
895 DiagMatrix (const RowVector& a);
896 DiagMatrix (const ColumnVector& a);
897 DiagMatrix (const DiagMatrix& a);
898 DiagMatrix (double a);
899 ~DiagMatrix (void);
900
901 DiagMatrix& operator = (const DiagMatrix& a);
902
903 int rows (void) const;
904 int cols (void) const;
905 int columns (void) const;
906
907 double& elem (int r, int c);
908 double& checkelem (int r, int c);
909 double& operator () (int r, int c);
910
911 double elem (int r, int c) const; // const access
912 double checkelem (int r, int c) const;
913 double operator () (int r, int c) const;
914
915 DiagMatrix& resize (int r, int c);
916 DiagMatrix& resize (int r, int c, double val);
917
918 int operator == (const DiagMatrix& a) const;
919 int operator != (const DiagMatrix& a) const;
920
921 DiagMatrix& fill (double val);
922 DiagMatrix& fill (double val, int beg, int end);
923 DiagMatrix& fill (const ColumnVector& a);
924 DiagMatrix& fill (const RowVector& a);
925 DiagMatrix& fill (const ColumnVector& a, int beg);
926 DiagMatrix& fill (const RowVector& a, int beg);
927
928 DiagMatrix transpose (void) const;
929
930 // resize is the destructive analog for this one
931
932 Matrix extract (int r1, int c1, int r2, int c2) const;
933
934 // extract row or column i.
935
936 RowVector row (int i) const;
937 RowVector row (char *s) const;
938
939 ColumnVector column (int i) const;
940 ColumnVector column (char *s) const;
941
942 DiagMatrix inverse (int& info) const;
943 DiagMatrix inverse (void) const;
944
945 // diagonal matrix by scalar -> matrix operations
946
947 Matrix operator + (double s) const;
948 Matrix operator - (double s) const;
949
950 ComplexMatrix operator + (Complex s) const;
951 ComplexMatrix operator - (Complex s) const;
952
953 // diagonal matrix by scalar -> diagonal matrix operations
954
955 DiagMatrix operator * (double s) const;
956 DiagMatrix operator / (double s) const;
957
958 ComplexDiagMatrix operator * (Complex s) const;
959 ComplexDiagMatrix operator / (Complex s) const;
960
961 // scalar by diagonal matrix -> matrix operations
962
963 friend Matrix operator + (double s, const DiagMatrix& a);
964 friend Matrix operator - (double s, const DiagMatrix& a);
965
966 // scalar by diagonal matrix -> diagonal matrix operations
967
968 friend DiagMatrix operator * (double s, const DiagMatrix& a);
969 friend DiagMatrix operator / (double s, const DiagMatrix& a);
970
971 // diagonal matrix by column vector -> column vector operations
972
973 ColumnVector operator * (const ColumnVector& a) const;
974
975 ComplexColumnVector operator * (const ComplexColumnVector& a) const;
976
977 // diagonal matrix by diagonal matrix -> diagonal matrix operations
978
979 DiagMatrix operator + (const DiagMatrix& a) const;
980 DiagMatrix operator - (const DiagMatrix& a) const;
981 DiagMatrix operator * (const DiagMatrix& a) const;
982
983 ComplexDiagMatrix operator + (const ComplexDiagMatrix& a) const;
984 ComplexDiagMatrix operator - (const ComplexDiagMatrix& a) const;
985 ComplexDiagMatrix operator * (const ComplexDiagMatrix& a) const;
986
987 DiagMatrix product (const DiagMatrix& a) const; // element by element
988 DiagMatrix quotient (const DiagMatrix& a) const; // element by element
989
990 ComplexDiagMatrix product (const ComplexDiagMatrix& a) const; // el by el
991 ComplexDiagMatrix quotient (const ComplexDiagMatrix& a) const; // el by el
992
993 DiagMatrix& operator += (const DiagMatrix& a);
994 DiagMatrix& operator -= (const DiagMatrix& a);
995
996 // diagonal matrix by matrix -> matrix operations
997
998 Matrix operator + (const Matrix& a) const;
999 Matrix operator - (const Matrix& a) const;
1000 Matrix operator * (const Matrix& a) const;
1001
1002 ComplexMatrix operator + (const ComplexMatrix& a) const;
1003 ComplexMatrix operator - (const ComplexMatrix& a) const;
1004 ComplexMatrix operator * (const ComplexMatrix& a) const;
1005
1006 // unary operations
1007
1008 DiagMatrix operator - (void) const;
1009
1010 ColumnVector diag (void) const;
1011 ColumnVector diag (int k) const;
1012
1013 // i/o
1014
1015 friend ostream& operator << (ostream& os, const DiagMatrix& a);
1016
1017 private:
1018 int nr;
1019 int nc;
1020 int len;
1021 double *data;
1022
1023 DiagMatrix (double *d, int nr, int nc);
1024 };
1025
1026 inline DiagMatrix::DiagMatrix (void)
1027 { nr = 0; nc = 0; len = 0; data = 0; }
1028
1029 inline DiagMatrix::DiagMatrix (double *d, int r, int c)
1030 { nr = r; nc = c; len = nr < nc ? nr : nc; data = d; }
1031
1032 inline DiagMatrix::~DiagMatrix (void) { delete [] data; data = 0; }
1033
1034 inline int DiagMatrix::rows (void) const { return len; }
1035 inline int DiagMatrix::cols (void) const { return len; }
1036 inline int DiagMatrix::columns (void) const { return len; }
1037
1038 // Would be nice to be able to avoid compiler warning and make this
1039 // fail on assignment.
1040 inline double& DiagMatrix::elem (int r, int c)
1041 { return (r == c) ? data[r] : 0; }
1042
1043 inline double& DiagMatrix::checkelem (int r, int c)
1044 {
1045 #ifndef NO_RANGE_CHECK
1046 if (r < 0 || r >= nr || c < 0 || c >= nc)
1047 FAIL;
1048 #endif
1049
1050 return elem (r, c);
1051 }
1052
1053 inline double& DiagMatrix::operator () (int r, int c)
1054 { return checkelem (r, c); }
1055
1056 inline double DiagMatrix::elem (int r, int c) const
1057 { return (r == c) ? data[r] : 0; }
1058
1059 inline double DiagMatrix::checkelem (int r, int c) const
1060 {
1061 #ifndef NO_RANGE_CHECK
1062 if (r < 0 || r >= nr || c < 0 || c >= nc)
1063 FAIL;
1064 #endif
1065
1066 return elem (r, c);
1067 }
1068
1069 inline double DiagMatrix::operator () (int r, int c) const
1070 { return checkelem (r, c); }
1071
1072 /*
1073 * Complex Matrix class
1074 */
1075
1076 class ComplexMatrix
1077 {
1078 friend class Matrix;
1079 friend class DiagMatrix;
1080 friend class ComplexRowVector;
1081 friend class ComplexDiagMatrix;
1082 friend class EIG;
1083 friend class ComplexHESS;
1084 friend class ComplexSVD;
1085 friend class ComplexSCHUR;
1086 friend class ComplexLU;
1087 friend class ComplexQR;
1088
1089 public:
1090 ComplexMatrix (void);
1091 ComplexMatrix (int r, int c);
1092 ComplexMatrix (int r, int c, double val);
1093 ComplexMatrix (int r, int c, Complex val);
1094 ComplexMatrix (const Matrix& a);
1095 ComplexMatrix (const ComplexMatrix& a);
1096 ComplexMatrix (const DiagMatrix& a);
1097 ComplexMatrix (const ComplexDiagMatrix& a);
1098 ComplexMatrix (double a);
1099 ComplexMatrix (Complex a);
1100 ~ComplexMatrix (void);
1101
1102 ComplexMatrix& operator = (const Matrix& a);
1103 ComplexMatrix& operator = (const ComplexMatrix& a);
1104
1105 int rows (void) const;
1106 int cols (void) const;
1107 int columns (void) const;
1108
1109 Complex& elem (int r, int c);
1110 Complex& checkelem (int r, int c);
1111 Complex& operator () (int r, int c);
1112
1113 Complex elem (int r, int c) const; // const access
1114 Complex checkelem (int r, int c) const;
1115 Complex operator () (int r, int c) const;
1116
1117 ComplexMatrix& resize (int r, int c);
1118 ComplexMatrix& resize (int r, int c, double val);
1119 ComplexMatrix& resize (int r, int c, Complex val);
1120
1121 int operator == (const ComplexMatrix& a) const;
1122 int operator != (const ComplexMatrix& a) const;
1123
1124 // destructive insert/delete/reorder operations
1125
1126 ComplexMatrix& insert (const Matrix& a, int r, int c);
1127 ComplexMatrix& insert (const RowVector& a, int r, int c);
1128 ComplexMatrix& insert (const ColumnVector& a, int r, int c);
1129 ComplexMatrix& insert (const DiagMatrix& a, int r, int c);
1130
1131 ComplexMatrix& insert (const ComplexMatrix& a, int r, int c);
1132 ComplexMatrix& insert (const ComplexRowVector& a, int r, int c);
1133 ComplexMatrix& insert (const ComplexColumnVector& a, int r, int c);
1134 ComplexMatrix& insert (const ComplexDiagMatrix& a, int r, int c);
1135
1136 ComplexMatrix& fill (double val);
1137 ComplexMatrix& fill (Complex val);
1138 ComplexMatrix& fill (double val, int r1, int c1, int r2, int c2);
1139 ComplexMatrix& fill (Complex val, int r1, int c1, int r2, int c2);
1140
1141 ComplexMatrix append (const Matrix& a) const;
1142 ComplexMatrix append (const RowVector& a) const;
1143 ComplexMatrix append (const ColumnVector& a) const;
1144 ComplexMatrix append (const DiagMatrix& a) const;
1145
1146 ComplexMatrix append (const ComplexMatrix& a) const;
1147 ComplexMatrix append (const ComplexRowVector& a) const;
1148 ComplexMatrix append (const ComplexColumnVector& a) const;
1149 ComplexMatrix append (const ComplexDiagMatrix& a) const;
1150
1151 ComplexMatrix stack (const Matrix& a) const;
1152 ComplexMatrix stack (const RowVector& a) const;
1153 ComplexMatrix stack (const ColumnVector& a) const;
1154 ComplexMatrix stack (const DiagMatrix& a) const;
1155
1156 ComplexMatrix stack (const ComplexMatrix& a) const;
1157 ComplexMatrix stack (const ComplexRowVector& a) const;
1158 ComplexMatrix stack (const ComplexColumnVector& a) const;
1159 ComplexMatrix stack (const ComplexDiagMatrix& a) const;
1160
1161 ComplexMatrix hermitian (void) const; // complex conjugate transpose
1162 ComplexMatrix transpose (void) const;
1163
1164 friend Matrix real (const ComplexMatrix& a);
1165 friend Matrix imag (const ComplexMatrix& a);
1166 friend ComplexMatrix conj (const ComplexMatrix& a);
1167
1168 // resize is the destructive equivalent for this one
1169
1170 ComplexMatrix extract (int r1, int c1, int r2, int c2) const;
1171
1172 // extract row or column i.
1173
1174 ComplexRowVector row (int i) const;
1175 ComplexRowVector row (char *s) const;
1176
1177 ComplexColumnVector column (int i) const;
1178 ComplexColumnVector column (char *s) const;
1179
1180 ComplexMatrix inverse (int& info, double& rcond) const;
1181 ComplexMatrix inverse (int& info) const;
1182 ComplexMatrix inverse (void) const;
1183
1184 ComplexMatrix fourier (void) const;
1185 ComplexMatrix ifourier (void) const;
1186
1187 ComplexDET determinant (void) const;
1188 ComplexDET determinant (int& info) const;
1189 ComplexDET determinant (int& info, double& rcond) const;
1190
1191 ComplexMatrix solve (const Matrix& b) const;
1192 ComplexMatrix solve (const Matrix& b, int& info) const;
1193 ComplexMatrix solve (const Matrix& b, int& info, double& rcond) const;
1194
1195 ComplexMatrix solve (const ComplexMatrix& b) const;
1196 ComplexMatrix solve (const ComplexMatrix& b, int& info) const;
1197 ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const;
1198
1199 ComplexColumnVector solve (const ColumnVector& b) const;
1200 ComplexColumnVector solve (const ColumnVector& b, int& info) const;
1201 ComplexColumnVector solve (const ColumnVector& b, int& info,
1202 double& rcond) const;
1203
1204 ComplexColumnVector solve (const ComplexColumnVector& b) const;
1205 ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const;
1206 ComplexColumnVector solve (const ComplexColumnVector& b, int& info,
1207 double& rcond) const;
1208
1209 ComplexMatrix lssolve (const Matrix& b) const;
1210 ComplexMatrix lssolve (const Matrix& b, int& info) const;
1211 ComplexMatrix lssolve (const Matrix& b, int& info, int& rank) const;
1212
1213 ComplexMatrix lssolve (const ComplexMatrix& b) const;
1214 ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const;
1215 ComplexMatrix lssolve (const ComplexMatrix& b, int& info,
1216 int& rank) const;
1217
1218 ComplexColumnVector lssolve (const ColumnVector& b) const;
1219 ComplexColumnVector lssolve (const ColumnVector& b, int& info) const;
1220 ComplexColumnVector lssolve (const ColumnVector& b, int& info,
1221 int& rank) const;
1222
1223 ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
1224 ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const;
1225 ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
1226 int& rank) const;
1227
1228 // matrix by scalar -> matrix operations
1229
1230 ComplexMatrix operator + (double s) const;
1231 ComplexMatrix operator - (double s) const;
1232 ComplexMatrix operator * (double s) const;
1233 ComplexMatrix operator / (double s) const;
1234
1235 ComplexMatrix operator + (Complex s) const;
1236 ComplexMatrix operator - (Complex s) const;
1237 ComplexMatrix operator * (Complex s) const;
1238 ComplexMatrix operator / (Complex s) const;
1239
1240 // scalar by matrix -> matrix operations
1241
1242 friend ComplexMatrix operator + (double s, const ComplexMatrix& a);
1243 friend ComplexMatrix operator - (double s, const ComplexMatrix& a);
1244 friend ComplexMatrix operator * (double s, const ComplexMatrix& a);
1245 friend ComplexMatrix operator / (double s, const ComplexMatrix& a);
1246
1247 friend ComplexMatrix operator + (Complex s, const ComplexMatrix& a);
1248 friend ComplexMatrix operator - (Complex s, const ComplexMatrix& a);
1249 friend ComplexMatrix operator * (Complex s, const ComplexMatrix& a);
1250 friend ComplexMatrix operator / (Complex s, const ComplexMatrix& a);
1251
1252 // matrix by column vector -> column vector operations
1253
1254 ComplexColumnVector operator * (const ColumnVector& a) const;
1255
1256 ComplexColumnVector operator * (const ComplexColumnVector& a) const;
1257
1258 // matrix by diagonal matrix -> matrix operations
1259
1260 ComplexMatrix operator + (const DiagMatrix& a) const;
1261 ComplexMatrix operator - (const DiagMatrix& a) const;
1262 ComplexMatrix operator * (const DiagMatrix& a) const;
1263
1264 ComplexMatrix operator + (const ComplexDiagMatrix& a) const;
1265 ComplexMatrix operator - (const ComplexDiagMatrix& a) const;
1266 ComplexMatrix operator * (const ComplexDiagMatrix& a) const;
1267
1268 ComplexMatrix& operator += (const DiagMatrix& a);
1269 ComplexMatrix& operator -= (const DiagMatrix& a);
1270
1271 ComplexMatrix& operator += (const ComplexDiagMatrix& a);
1272 ComplexMatrix& operator -= (const ComplexDiagMatrix& a);
1273
1274 // matrix by matrix -> matrix operations
1275
1276 ComplexMatrix operator + (const Matrix& a) const;
1277 ComplexMatrix operator - (const Matrix& a) const;
1278 ComplexMatrix operator * (const Matrix& a) const;
1279
1280 ComplexMatrix operator + (const ComplexMatrix& a) const;
1281 ComplexMatrix operator - (const ComplexMatrix& a) const;
1282 ComplexMatrix operator * (const ComplexMatrix& a) const;
1283
1284 ComplexMatrix product (const Matrix& a) const; // element by element
1285 ComplexMatrix quotient (const Matrix& a) const; // element by element
1286
1287 ComplexMatrix product (const ComplexMatrix& a) const; // element by element
1288 ComplexMatrix quotient (const ComplexMatrix& a) const; // element by element
1289
1290 ComplexMatrix& operator += (const Matrix& a);
1291 ComplexMatrix& operator -= (const Matrix& a);
1292
1293 ComplexMatrix& operator += (const ComplexMatrix& a);
1294 ComplexMatrix& operator -= (const ComplexMatrix& a);
1295
1296 // unary operations
1297
1298 ComplexMatrix operator - (void) const;
1299 Matrix operator ! (void) const;
1300
1301 // other operations
1302
1303 friend ComplexMatrix map (c_c_Mapper f, const ComplexMatrix& a);
1304 friend Matrix map (d_c_Mapper f, const ComplexMatrix& a);
1305 void map (c_c_Mapper f);
1306
1307 Matrix all (void) const;
1308 Matrix any (void) const;
1309
1310 ComplexMatrix cumprod (void) const;
1311 ComplexMatrix cumsum (void) const;
1312 ComplexMatrix prod (void) const;
1313 ComplexMatrix sum (void) const;
1314 ComplexMatrix sumsq (void) const;
1315
1316 ComplexColumnVector diag (void) const;
1317 ComplexColumnVector diag (int k) const;
1318
1319 ComplexColumnVector row_min (void) const;
1320 ComplexColumnVector row_max (void) const;
1321
1322 ComplexRowVector column_min (void) const;
1323 ComplexRowVector column_max (void) const;
1324
1325 // i/o
1326
1327 friend ostream& operator << (ostream& os, const ComplexMatrix& a);
1328 friend istream& operator >> (istream& is, ComplexMatrix& a);
1329
1330 // conversions
1331
1332 Complex *fortran_vec (void);
1333
1334 private:
1335 int nr;
1336 int nc;
1337 int len;
1338 Complex *data;
1339
1340 ComplexMatrix (Complex *d, int r, int c);
1341 };
1342
1343 inline ComplexMatrix::ComplexMatrix (void)
1344 { nr = 0; nc = 0; len = 0; data = 0; }
1345
1346 inline ComplexMatrix::ComplexMatrix (Complex *d, int r, int c)
1347 { nr = r; nc = c; len = nr*nc; data = d; }
1348
1349 inline ComplexMatrix::~ComplexMatrix (void) { delete [] data; data = 0; }
1350
1351 inline int ComplexMatrix::rows (void) const { return nr; }
1352 inline int ComplexMatrix::cols (void) const { return nc; }
1353 inline int ComplexMatrix::columns (void) const { return nc; }
1354
1355 inline Complex& ComplexMatrix::elem (int r, int c) { return data[nr*c+r]; }
1356
1357 inline Complex& ComplexMatrix::checkelem (int r, int c)
1358 {
1359 #ifndef NO_RANGE_CHECK
1360 if (r < 0 || r >= nr || c < 0 || c >= nc)
1361 FAIL;
1362 #endif
1363
1364 return elem (r, c);
1365 }
1366
1367 inline Complex& ComplexMatrix::operator () (int r, int c)
1368 { return checkelem (r, c); }
1369
1370 inline Complex ComplexMatrix::elem (int r, int c) const
1371 { return data[nr*c+r]; }
1372
1373 inline Complex ComplexMatrix::checkelem (int r, int c) const
1374 {
1375 #ifndef NO_RANGE_CHECK
1376 if (r < 0 || r >= nr || c < 0 || c >= nc)
1377 FAIL;
1378 #endif
1379
1380 return elem (r, c);
1381 }
1382
1383 inline Complex ComplexMatrix::operator () (int r, int c) const
1384 { return checkelem (r, c); }
1385
1386 inline Complex *ComplexMatrix::fortran_vec (void) { return data; }
1387
1388 /*
1389 * Complex Column Vector class
1390 */
1391
1392 class ComplexColumnVector
1393 {
1394 friend class DiagMatrix;
1395 friend class ComplexMatrix;
1396 friend class ColumnVector;
1397 friend class ComplexDiagMatrix;
1398
1399 public:
1400 ComplexColumnVector (void);
1401 ComplexColumnVector (int n);
1402 ComplexColumnVector (int n, double val);
1403 ComplexColumnVector (int n, Complex val);
1404 ComplexColumnVector (const ColumnVector& a);
1405 ComplexColumnVector (const ComplexColumnVector& a);
1406 ComplexColumnVector (double a);
1407 ComplexColumnVector (Complex a);
1408 ~ComplexColumnVector (void);
1409
1410 ComplexColumnVector& operator = (const ColumnVector& a);
1411 ComplexColumnVector& operator = (const ComplexColumnVector& a);
1412
1413 int capacity (void) const;
1414 int length (void) const;
1415
1416 Complex& elem (int n);
1417 Complex& checkelem (int n);
1418 Complex& operator () (int n);
1419
1420 Complex elem (int n) const; // const access
1421 Complex checkelem (int n) const;
1422 Complex operator () (int n) const;
1423
1424 ComplexColumnVector& resize (int n);
1425 ComplexColumnVector& resize (int n, double val);
1426 ComplexColumnVector& resize (int n, Complex val);
1427
1428 int operator == (const ComplexColumnVector& a) const;
1429 int operator != (const ComplexColumnVector& a) const;
1430
1431 // destructive insert/delete/reorder operations
1432
1433 ComplexColumnVector& insert (const ColumnVector& a, int r);
1434 ComplexColumnVector& insert (const ComplexColumnVector& a, int r);
1435
1436 ComplexColumnVector& fill (double val);
1437 ComplexColumnVector& fill (Complex val);
1438 ComplexColumnVector& fill (double val, int r1, int r2);
1439 ComplexColumnVector& fill (Complex val, int r1, int r2);
1440
1441 ComplexColumnVector stack (const ColumnVector& a) const;
1442 ComplexColumnVector stack (const ComplexColumnVector& a) const;
1443
1444 ComplexRowVector hermitian (void) const; // complex conjugate transpose.
1445 ComplexRowVector transpose (void) const;
1446
1447 friend ColumnVector real (const ComplexColumnVector& a);
1448 friend ColumnVector imag (const ComplexColumnVector& a);
1449 friend ComplexColumnVector conj (const ComplexColumnVector& a);
1450
1451 // resize is the destructive equivalent for this one
1452
1453 ComplexColumnVector extract (int r1, int r2) const;
1454
1455 // column vector by scalar -> column vector operations
1456
1457 ComplexColumnVector operator + (double s) const;
1458 ComplexColumnVector operator - (double s) const;
1459 ComplexColumnVector operator * (double s) const;
1460 ComplexColumnVector operator / (double s) const;
1461
1462 ComplexColumnVector operator + (Complex s) const;
1463 ComplexColumnVector operator - (Complex s) const;
1464 ComplexColumnVector operator * (Complex s) const;
1465 ComplexColumnVector operator / (Complex s) const;
1466
1467 // scalar by column vector -> column vector operations
1468
1469 friend ComplexColumnVector operator + (double s,
1470 const ComplexColumnVector& a);
1471 friend ComplexColumnVector operator - (double s,
1472 const ComplexColumnVector& a);
1473 friend ComplexColumnVector operator * (double s,
1474 const ComplexColumnVector& a);
1475 friend ComplexColumnVector operator / (double s,
1476 const ComplexColumnVector& a);
1477
1478 friend ComplexColumnVector operator + (Complex s,
1479 const ComplexColumnVector& a);
1480 friend ComplexColumnVector operator - (Complex s,
1481 const ComplexColumnVector& a);
1482 friend ComplexColumnVector operator * (Complex s,
1483 const ComplexColumnVector& a);
1484 friend ComplexColumnVector operator / (Complex s,
1485 const ComplexColumnVector& a);
1486
1487 // column vector by row vector -> matrix operations
1488
1489 ComplexMatrix operator * (const RowVector& a) const;
1490
1491 ComplexMatrix operator * (const ComplexRowVector& a) const;
1492
1493 // column vector by column vector -> column vector operations
1494
1495 ComplexColumnVector operator + (const ColumnVector& a) const;
1496 ComplexColumnVector operator - (const ColumnVector& a) const;
1497
1498 ComplexColumnVector operator + (const ComplexColumnVector& a) const;
1499 ComplexColumnVector operator - (const ComplexColumnVector& a) const;
1500
1501 ComplexColumnVector product (const ColumnVector& a) const; // el by el
1502 ComplexColumnVector quotient (const ColumnVector& a) const; // el by el
1503
1504 ComplexColumnVector product (const ComplexColumnVector& a) const;
1505 ComplexColumnVector quotient (const ComplexColumnVector& a) const;
1506
1507 ComplexColumnVector& operator += (const ColumnVector& a);
1508 ComplexColumnVector& operator -= (const ColumnVector& a);
1509
1510 ComplexColumnVector& operator += (const ComplexColumnVector& a);
1511 ComplexColumnVector& operator -= (const ComplexColumnVector& a);
1512
1513 // unary operations
1514
1515 ComplexColumnVector operator - (void) const;
1516
1517 friend ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a);
1518 friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a);
1519 void map (c_c_Mapper f);
1520
1521 Complex min (void) const;
1522 Complex max (void) const;
1523
1524 // i/o
1525
1526 friend ostream& operator << (ostream& os, const ComplexColumnVector& a);
1527
1528 // conversions
1529
1530 Complex *fortran_vec (void);
1531
1532 private:
1533 int len;
1534 Complex *data;
1535
1536 ComplexColumnVector (Complex *d, int l);
1537 };
1538
1539 inline ComplexColumnVector::ComplexColumnVector (void) { len = 0; data = 0; }
1540 inline ComplexColumnVector::ComplexColumnVector (Complex *d, int l)
1541 { len = l; data = d; }
1542 inline ComplexColumnVector::~ComplexColumnVector (void)
1543 { delete [] data; data = 0; }
1544
1545 inline int ComplexColumnVector::capacity (void) const { return len; }
1546 inline int ComplexColumnVector::length (void) const { return len; }
1547
1548 inline Complex& ComplexColumnVector::elem (int n) { return data[n]; }
1549
1550 inline Complex&
1551 ComplexColumnVector::checkelem (int n)
1552 {
1553 #ifndef NO_RANGE_CHECK
1554 if (n < 0 || n >= len)
1555 FAIL;
1556 #endif
1557
1558 return elem (n);
1559 }
1560
1561 inline Complex& ComplexColumnVector::operator () (int n)
1562 { return checkelem (n); }
1563
1564 inline Complex ComplexColumnVector::elem (int n) const { return data[n]; }
1565
1566 inline Complex
1567 ComplexColumnVector::checkelem (int n) const
1568 {
1569 #ifndef NO_RANGE_CHECK
1570 if (n < 0 || n >= len)
1571 FAIL;
1572 #endif
1573
1574 return elem (n);
1575 }
1576
1577 inline Complex ComplexColumnVector::operator () (int n) const
1578 { return checkelem (n); }
1579
1580 inline Complex *ComplexColumnVector::fortran_vec (void) { return data; }
1581
1582 /*
1583 * Complex Row Vector class
1584 */
1585
1586 class ComplexRowVector
1587 {
1588 friend class RowVector;
1589 friend class ComplexMatrix;
1590 friend class ComplexColumnVector;
1591 friend class ComplexDiagMatrix;
1592
1593 public:
1594 ComplexRowVector (void);
1595 ComplexRowVector (int n);
1596 ComplexRowVector (int n, double val);
1597 ComplexRowVector (int n, Complex val);
1598 ComplexRowVector (const RowVector& a);
1599 ComplexRowVector (const ComplexRowVector& a);
1600 ComplexRowVector (double a);
1601 ComplexRowVector (Complex a);
1602 ~ComplexRowVector (void);
1603
1604 ComplexRowVector& operator = (const RowVector& a);
1605 ComplexRowVector& operator = (const ComplexRowVector& a);
1606
1607 int capacity (void) const;
1608 int length (void) const;
1609
1610 Complex& checkelem (int n);
1611 Complex& elem (int n);
1612 Complex& operator () (int n);
1613
1614 Complex checkelem (int n) const; // const access
1615 Complex elem (int n) const;
1616 Complex operator () (int n) const;
1617
1618 ComplexRowVector& resize (int n);
1619 ComplexRowVector& resize (int n, double val);
1620 ComplexRowVector& resize (int n, Complex val);
1621
1622 int operator == (const ComplexRowVector& a) const;
1623 int operator != (const ComplexRowVector& a) const;
1624
1625 // destructive insert/delete/reorder operations
1626
1627 ComplexRowVector& insert (const RowVector& a, int c);
1628 ComplexRowVector& insert (const ComplexRowVector& a, int c);
1629
1630 ComplexRowVector& fill (double val);
1631 ComplexRowVector& fill (Complex val);
1632 ComplexRowVector& fill (double val, int c1, int c2);
1633 ComplexRowVector& fill (Complex val, int c1, int c2);
1634
1635 ComplexRowVector append (const RowVector& a) const;
1636 ComplexRowVector append (const ComplexRowVector& a) const;
1637
1638 ComplexColumnVector hermitian (void) const; // complex conjugate transpose.
1639 ComplexColumnVector transpose (void) const;
1640
1641 friend RowVector real (const ComplexRowVector& a);
1642 friend RowVector imag (const ComplexRowVector& a);
1643 friend ComplexRowVector conj (const ComplexRowVector& a);
1644
1645 // resize is the destructive equivalent for this one
1646
1647 ComplexRowVector extract (int c1, int c2) const;
1648
1649 // row vector by scalar -> row vector operations
1650
1651 ComplexRowVector operator + (double s) const;
1652 ComplexRowVector operator - (double s) const;
1653 ComplexRowVector operator * (double s) const;
1654 ComplexRowVector operator / (double s) const;
1655
1656 ComplexRowVector operator + (Complex s) const;
1657 ComplexRowVector operator - (Complex s) const;
1658 ComplexRowVector operator * (Complex s) const;
1659 ComplexRowVector operator / (Complex s) const;
1660
1661 // scalar by row vector -> row vector operations
1662
1663 friend ComplexRowVector operator + (double s, const ComplexRowVector& a);
1664 friend ComplexRowVector operator - (double s, const ComplexRowVector& a);
1665 friend ComplexRowVector operator * (double s, const ComplexRowVector& a);
1666 friend ComplexRowVector operator / (double s, const ComplexRowVector& a);
1667
1668 friend ComplexRowVector operator + (Complex s, const ComplexRowVector& a);
1669 friend ComplexRowVector operator - (Complex s, const ComplexRowVector& a);
1670 friend ComplexRowVector operator * (Complex s, const ComplexRowVector& a);
1671 friend ComplexRowVector operator / (Complex s, const ComplexRowVector& a);
1672
1673 // row vector by column vector -> scalar
1674
1675 Complex operator * (const ColumnVector& a) const;
1676
1677 Complex operator * (const ComplexColumnVector& a) const;
1678
1679 // row vector by matrix -> row vector
1680
1681 ComplexRowVector operator * (const Matrix& a) const;
1682
1683 ComplexRowVector operator * (const ComplexMatrix& a) const;
1684
1685 // row vector by row vector -> row vector operations
1686
1687 ComplexRowVector operator + (const RowVector& a) const;
1688 ComplexRowVector operator - (const RowVector& a) const;
1689
1690 ComplexRowVector operator + (const ComplexRowVector& a) const;
1691 ComplexRowVector operator - (const ComplexRowVector& a) const;
1692
1693 ComplexRowVector product (const RowVector& a) const; // element by element
1694 ComplexRowVector quotient (const RowVector& a) const; // element by element
1695
1696 ComplexRowVector product (const ComplexRowVector& a) const; // el by el
1697 ComplexRowVector quotient (const ComplexRowVector& a) const; // el by el
1698
1699 ComplexRowVector& operator += (const RowVector& a);
1700 ComplexRowVector& operator -= (const RowVector& a);
1701
1702 ComplexRowVector& operator += (const ComplexRowVector& a);
1703 ComplexRowVector& operator -= (const ComplexRowVector& a);
1704
1705 // unary operations
1706
1707 ComplexRowVector operator - (void) const;
1708
1709 friend ComplexRowVector map (c_c_Mapper f, const ComplexRowVector& a);
1710 friend RowVector map (d_c_Mapper f, const ComplexRowVector& a);
1711 void map (c_c_Mapper f);
1712
1713 Complex min (void) const;
1714 Complex max (void) const;
1715
1716 // i/o
1717
1718 friend ostream& operator << (ostream& os, const ComplexRowVector& a);
1719
1720 // conversions
1721
1722 Complex *fortran_vec (void);
1723
1724 private:
1725 int len;
1726 Complex *data;
1727
1728 ComplexRowVector (Complex *d, int l);
1729 };
1730
1731 inline ComplexRowVector::ComplexRowVector (void) { len = 0; data = 0; }
1732 inline ComplexRowVector::ComplexRowVector (Complex *d, int l)
1733 { len = l; data = d; }
1734 inline ComplexRowVector::~ComplexRowVector (void) { delete [] data; data = 0; }
1735
1736 inline int ComplexRowVector::capacity (void) const { return len; }
1737 inline int ComplexRowVector::length (void) const { return len; }
1738
1739 inline Complex& ComplexRowVector::elem (int n) { return data[n]; }
1740
1741 inline Complex&
1742 ComplexRowVector::checkelem (int n)
1743 {
1744 #ifndef NO_RANGE_CHECK
1745 if (n < 0 || n >= len)
1746 FAIL;
1747 #endif
1748
1749 return elem (n);
1750 }
1751
1752 inline Complex& ComplexRowVector::operator () (int n) { return checkelem (n); }
1753
1754 inline Complex ComplexRowVector::elem (int n) const { return data[n]; }
1755
1756 inline Complex
1757 ComplexRowVector::checkelem (int n) const
1758 {
1759 #ifndef NO_RANGE_CHECK
1760 if (n < 0 || n >= len)
1761 FAIL;
1762 #endif
1763
1764 return elem (n);
1765 }
1766
1767 inline Complex ComplexRowVector::operator () (int n) const
1768 { return checkelem (n); }
1769
1770 inline Complex *ComplexRowVector::fortran_vec (void) { return data; }
1771
1772 /*
1773 * Complex Diagonal Matrix class
1774 */
1775
1776 class ComplexDiagMatrix
1777 {
1778 friend class Matrix;
1779 friend class DiagMatrix;
1780 friend class ComplexMatrix;
1781
1782 public:
1783 ComplexDiagMatrix (void);
1784 ComplexDiagMatrix (int n);
1785 ComplexDiagMatrix (int n, double val);
1786 ComplexDiagMatrix (int n, Complex val);
1787 ComplexDiagMatrix (int r, int c);
1788 ComplexDiagMatrix (int r, int c, double val);
1789 ComplexDiagMatrix (int r, int c, Complex val);
1790 ComplexDiagMatrix (const RowVector& a);
1791 ComplexDiagMatrix (const ComplexRowVector& a);
1792 ComplexDiagMatrix (const ColumnVector& a);
1793 ComplexDiagMatrix (const ComplexColumnVector& a);
1794 ComplexDiagMatrix (const DiagMatrix& a);
1795 ComplexDiagMatrix (const ComplexDiagMatrix& a);
1796 ComplexDiagMatrix (double a);
1797 ComplexDiagMatrix (Complex a);
1798 ~ComplexDiagMatrix (void);
1799
1800 ComplexDiagMatrix& operator = (const DiagMatrix& a);
1801 ComplexDiagMatrix& operator = (const ComplexDiagMatrix& a);
1802
1803 int rows (void) const;
1804 int cols (void) const;
1805 int columns (void) const;
1806
1807 Complex& checkelem (int r, int c);
1808 Complex& elem (int r, int c);
1809 Complex& operator () (int r, int c);
1810
1811 Complex checkelem (int r, int c) const; // const access
1812 Complex elem (int r, int c) const;
1813 Complex operator () (int r, int c) const;
1814
1815 ComplexDiagMatrix& resize (int r, int c);
1816 ComplexDiagMatrix& resize (int r, int c, double val);
1817 ComplexDiagMatrix& resize (int r, int c, Complex val);
1818
1819 int operator == (const ComplexDiagMatrix& a) const;
1820 int operator != (const ComplexDiagMatrix& a) const;
1821
1822 ComplexDiagMatrix& fill (double val);
1823 ComplexDiagMatrix& fill (Complex val);
1824 ComplexDiagMatrix& fill (double val, int beg, int end);
1825 ComplexDiagMatrix& fill (Complex val, int beg, int end);
1826 ComplexDiagMatrix& fill (const ColumnVector& a);
1827 ComplexDiagMatrix& fill (const ComplexColumnVector& a);
1828 ComplexDiagMatrix& fill (const RowVector& a);
1829 ComplexDiagMatrix& fill (const ComplexRowVector& a);
1830 ComplexDiagMatrix& fill (const ColumnVector& a, int beg);
1831 ComplexDiagMatrix& fill (const ComplexColumnVector& a, int beg);
1832 ComplexDiagMatrix& fill (const RowVector& a, int beg);
1833 ComplexDiagMatrix& fill (const ComplexRowVector& a, int beg);
1834
1835 ComplexDiagMatrix hermitian (void) const; // complex conjugate transpose
1836 ComplexDiagMatrix transpose (void) const;
1837
1838 friend DiagMatrix real (const ComplexDiagMatrix& a);
1839 friend DiagMatrix imag (const ComplexDiagMatrix& a);
1840 friend ComplexDiagMatrix conj (const ComplexDiagMatrix& a);
1841
1842 // resize is the destructive analog for this one
1843
1844 ComplexMatrix extract (int r1, int c1, int r2, int c2) const;
1845
1846 // extract row or column i.
1847
1848 ComplexRowVector row (int i) const;
1849 ComplexRowVector row (char *s) const;
1850
1851 ComplexColumnVector column (int i) const;
1852 ComplexColumnVector column (char *s) const;
1853
1854 ComplexDiagMatrix inverse (int& info) const;
1855 ComplexDiagMatrix inverse (void) const;
1856
1857 // diagonal matrix by scalar -> matrix operations
1858
1859 ComplexMatrix operator + (double s) const;
1860 ComplexMatrix operator - (double s) const;
1861
1862 ComplexMatrix operator + (Complex s) const;
1863 ComplexMatrix operator - (Complex s) const;
1864
1865 // diagonal matrix by scalar -> diagonal matrix operations
1866
1867 ComplexDiagMatrix operator * (double s) const;
1868 ComplexDiagMatrix operator / (double s) const;
1869
1870 ComplexDiagMatrix operator * (Complex s) const;
1871 ComplexDiagMatrix operator / (Complex s) const;
1872
1873 // scalar by diagonal matrix -> matrix operations
1874
1875 friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a);
1876 friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a);
1877
1878 friend ComplexMatrix operator + (Complex s, const ComplexDiagMatrix& a);
1879 friend ComplexMatrix operator - (Complex s, const ComplexDiagMatrix& a);
1880
1881 // scalar by diagonal matrix -> diagonal matrix operations
1882
1883 friend ComplexDiagMatrix operator * (double s, const ComplexDiagMatrix& a);
1884 friend ComplexDiagMatrix operator / (double s, const ComplexDiagMatrix& a);
1885
1886 friend ComplexDiagMatrix operator * (Complex s, const ComplexDiagMatrix& a);
1887 friend ComplexDiagMatrix operator / (Complex s, const ComplexDiagMatrix& a);
1888
1889 // diagonal matrix by column vector -> column vector operations
1890
1891 ComplexColumnVector operator * (const ColumnVector& a) const;
1892
1893 ComplexColumnVector operator * (const ComplexColumnVector& a) const;
1894
1895 // diagonal matrix by diagonal matrix -> diagonal matrix operations
1896
1897 ComplexDiagMatrix operator + (const DiagMatrix& a) const;
1898 ComplexDiagMatrix operator - (const DiagMatrix& a) const;
1899 ComplexDiagMatrix operator * (const DiagMatrix& a) const;
1900
1901 ComplexDiagMatrix operator + (const ComplexDiagMatrix& a) const;
1902 ComplexDiagMatrix operator - (const ComplexDiagMatrix& a) const;
1903 ComplexDiagMatrix operator * (const ComplexDiagMatrix& a) const;
1904
1905 ComplexDiagMatrix product (const DiagMatrix& a) const; // element by element
1906 ComplexDiagMatrix quotient (const DiagMatrix& a) const; // element by element
1907
1908 ComplexDiagMatrix product (const ComplexDiagMatrix& a) const; // el by el
1909 ComplexDiagMatrix quotient (const ComplexDiagMatrix& a) const; // el by el
1910
1911 ComplexDiagMatrix& operator += (const DiagMatrix& a);
1912 ComplexDiagMatrix& operator -= (const DiagMatrix& a);
1913
1914 ComplexDiagMatrix& operator += (const ComplexDiagMatrix& a);
1915 ComplexDiagMatrix& operator -= (const ComplexDiagMatrix& a);
1916
1917 // diagonal matrix by matrix -> matrix operations
1918
1919 ComplexMatrix operator + (const Matrix& a) const;
1920 ComplexMatrix operator - (const Matrix& a) const;
1921 ComplexMatrix operator * (const Matrix& a) const;
1922
1923 ComplexMatrix operator + (const ComplexMatrix& a) const;
1924 ComplexMatrix operator - (const ComplexMatrix& a) const;
1925 ComplexMatrix operator * (const ComplexMatrix& a) const;
1926
1927 // unary operations
1928
1929 ComplexDiagMatrix operator - (void) const;
1930
1931 ComplexColumnVector diag (void) const;
1932 ComplexColumnVector diag (int k) const;
1933
1934 // i/o
1935
1936 friend ostream& operator << (ostream& os, const ComplexDiagMatrix& a);
1937
1938 private:
1939 int nr;
1940 int nc;
1941 int len;
1942 Complex *data;
1943
1944 ComplexDiagMatrix (Complex *d, int nr, int nc);
1945 };
1946
1947 inline ComplexDiagMatrix::ComplexDiagMatrix (void)
1948 { nr = 0; nc = 0; len = 0; data = 0; }
1949
1950 inline ComplexDiagMatrix::ComplexDiagMatrix (Complex *d, int r, int c)
1951 { nr = r; nc = c; len = nr < nc ? nr : nc; data = d; }
1952
1953 inline ComplexDiagMatrix::~ComplexDiagMatrix (void)
1954 { delete [] data; data = 0; }
1955
1956 inline int ComplexDiagMatrix::rows (void) const { return len; }
1957 inline int ComplexDiagMatrix::cols (void) const { return len; }
1958 inline int ComplexDiagMatrix::columns (void) const { return len; }
1959
1960 // Would be nice to be able to avoid compiler warning and make this
1961 // fail on assignment.
1962 inline Complex& ComplexDiagMatrix::elem (int r, int c)
1963 { Complex czero (0.0, 0.0); return (r == c) ? data[r] : czero; }
1964
1965 inline Complex& ComplexDiagMatrix::checkelem (int r, int c)
1966 {
1967 #ifndef NO_RANGE_CHECK
1968 if (r < 0 || r >= nr || c < 0 || c >= nc)
1969 FAIL;
1970 #endif
1971
1972 return elem (r, c);
1973 }
1974
1975 inline Complex& ComplexDiagMatrix::operator () (int r, int c)
1976 { return checkelem (r, c); }
1977
1978 inline Complex ComplexDiagMatrix::elem (int r, int c) const
1979 { Complex czero (0.0, 0.0); return (r == c) ? data[r] : czero; }
1980
1981 inline Complex ComplexDiagMatrix::checkelem (int r, int c) const
1982 {
1983 #ifndef NO_RANGE_CHECK
1984 if (r < 0 || r >= nr || c < 0 || c >= nc)
1985 FAIL;
1986 #endif
1987
1988 return elem (r, c);
1989 }
1990
1991 inline Complex ComplexDiagMatrix::operator () (int r, int c) const
1992 { return checkelem (r, c); }
1993
1994 /*
1995 * Result of a Determinant calculation.
1996 */
1997
1998 class DET
1999 {
2000 public:
2001 DET (void) {}
2002
2003 DET (const DET& a);
2004
2005 DET& operator = (const DET& a);
2006
2007 int value_will_overflow (void) const;
2008 int value_will_underflow (void) const;
2009 double coefficient (void) const;
2010 int exponent (void) const;
2011 double value (void) const;
2012
2013 friend ostream& operator << (ostream& os, const DET& a);
2014
2015 private:
2016 DET (const double *d);
2017
2018 double det [2];
2019 };
2020
2021 inline DET::DET (const DET& a) { det[0] = a.det[0]; det[1] = a.det[1]; }
2022
2023 inline DET& DET::operator = (const DET& a)
2024 { det[0] = a.det[0]; det[1] = a.det[1]; return *this; }
2025
2026 inline int DET::value_will_overflow (void) const
2027 { return det[2] + 1 > log10 (MAXDOUBLE) ? 1 : 0; }
2028
2029 inline int DET::value_will_underflow (void) const
2030 { return det[2] - 1 < log10 (MINDOUBLE) ? 1 : 0; }
2031
2032 inline double DET::coefficient (void) const { return det[0]; }
2033 inline int DET::exponent (void) const { return (int) det[1]; }
2034 inline double DET::value (void) const { return det[0] * pow (10.0, det[1]); }
2035
2036 inline DET::DET (const double *d) { det[0] = d[0]; det[1] = d[1]; }
2037
2038 /*
2039 * Result of a Determinant calculation.
2040 */
2041
2042 class ComplexDET
2043 {
2044 public:
2045 ComplexDET (void) {}
2046
2047 ComplexDET (const ComplexDET& a);
2048
2049 ComplexDET& operator = (const ComplexDET& a);
2050
2051 int value_will_overflow (void) const;
2052 int value_will_underflow (void) const;
2053 Complex coefficient (void) const;
2054 int exponent (void) const;
2055 Complex value (void) const;
2056
2057 friend ostream& operator << (ostream& os, const ComplexDET& a);
2058
2059 private:
2060 ComplexDET (const Complex *d);
2061
2062 Complex det [2];
2063 };
2064
2065 inline ComplexDET::ComplexDET (const ComplexDET& a)
2066 { det[0] = a.det[0]; det[1] = a.det[1]; }
2067
2068 inline ComplexDET& ComplexDET::operator = (const ComplexDET& a)
2069 { det[0] = a.det[0]; det[1] = a.det[1]; return *this; }
2070
2071 inline int ComplexDET::value_will_overflow (void) const
2072 { return real (det[2]) + 1 > log10 (MAXDOUBLE) ? 1 : 0; }
2073
2074 inline int ComplexDET::value_will_underflow (void) const
2075 { return real (det[2]) - 1 < log10 (MINDOUBLE) ? 1 : 0; }
2076
2077 inline Complex ComplexDET::coefficient (void) const { return det[0]; }
2078
2079 inline int ComplexDET::exponent (void) const { return (int) real (det[1]); }
2080
2081 inline Complex ComplexDET::value (void) const
2082 { return det[0] * pow (10.0, real (det[1])); }
2083
2084 inline ComplexDET::ComplexDET (const Complex *d)
2085 { det[0] = d[0]; det[1] = d[1]; }
2086
2087 /*
2088 * Result of a Hessenbug Decomposition
2089 */
2090
2091 class HESS
2092 {
2093 friend class Matrix;
2094
2095 public:
2096 HESS (void) {}
2097
2098 HESS (const Matrix& a);
2099 HESS (const Matrix&a, int& info);
2100
2101 HESS (const HESS& a);
2102
2103 HESS& operator = (const HESS& a);
2104 Matrix hess_matrix (void) const;
2105 Matrix unitary_hess_matrix (void) const;
2106 friend ostream& operator << (ostream& os, const HESS& a);
2107
2108 private:
2109 int init (const Matrix& a);
2110
2111 Matrix hess_mat;
2112 Matrix unitary_hess_mat;
2113 };
2114
2115 inline HESS::HESS (const Matrix& a) {init (a); }
2116 inline HESS::HESS (const Matrix& a, int& info) { info = init(a); }
2117 inline HESS::HESS (const HESS& a)
2118 {
2119 hess_mat = a.hess_mat;
2120 unitary_hess_mat = a.unitary_hess_mat;
2121 }
2122 inline HESS&
2123 HESS::operator = (const HESS& a)
2124 {
2125 hess_mat = a.hess_mat;
2126 unitary_hess_mat = a.unitary_hess_mat;
2127
2128 return *this;
2129 }
2130 inline Matrix HESS::hess_matrix (void) const { return hess_mat; }
2131 inline Matrix HESS::unitary_hess_matrix (void) const {return unitary_hess_mat;}
2132
2133 /*
2134 * Result of a Hessenburg Decomposition
2135 */
2136
2137 class ComplexHESS
2138 {
2139 friend class ComplexMatrix;
2140
2141 public:
2142 ComplexHESS (void) {}
2143 ComplexHESS (const ComplexMatrix& a);
2144 ComplexHESS (const ComplexMatrix& a, int& info);
2145 ComplexHESS (const ComplexHESS& a);
2146 ComplexHESS& operator = (const ComplexHESS& a);
2147 ComplexMatrix hess_matrix (void) const;
2148 ComplexMatrix unitary_hess_matrix (void) const;
2149
2150 friend ostream& operator << (ostream& os, const ComplexHESS& a);
2151
2152 private:
2153 int init (const ComplexMatrix& a);
2154
2155 ComplexMatrix hess_mat;
2156 ComplexMatrix unitary_hess_mat;
2157 };
2158
2159 inline ComplexHESS::ComplexHESS (const ComplexMatrix& a) { init(a); }
2160 inline ComplexHESS::ComplexHESS (const ComplexMatrix& a, int& info)
2161 { info = init(a); }
2162
2163 inline ComplexHESS::ComplexHESS (const ComplexHESS& a)
2164 {
2165 hess_mat = a.hess_mat;
2166 unitary_hess_mat = a.unitary_hess_mat;
2167 }
2168
2169 inline ComplexHESS&
2170 ComplexHESS::operator = (const ComplexHESS& a)
2171 {
2172 hess_mat = a.hess_mat;
2173 unitary_hess_mat = a.unitary_hess_mat;
2174
2175 return *this;
2176 }
2177
2178 inline ComplexMatrix ComplexHESS::hess_matrix (void) const
2179 { return hess_mat; }
2180
2181 inline ComplexMatrix ComplexHESS::unitary_hess_matrix (void) const
2182 { return unitary_hess_mat; }
2183
2184 /*
2185 * Result of a Schur Decomposition
2186 */
2187
2188 class SCHUR
2189 {
2190 friend class Matrix;
2191
2192 public:
2193 SCHUR (void) {}
2194
2195 SCHUR (const Matrix& a, const char *ord);
2196 SCHUR (const Matrix& a, const char *ord, int& info);
2197
2198 SCHUR (const SCHUR& a, const char *ord);
2199
2200 SCHUR& operator = (const SCHUR& a, const char *ord);
2201
2202 Matrix schur_matrix (void) const;
2203 Matrix unitary_matrix (void) const;
2204
2205 friend ostream& operator << (ostream& os, const SCHUR& a);
2206
2207 private:
2208 int init (const Matrix& a, const char *ord);
2209
2210 Matrix schur_mat;
2211 Matrix unitary_mat;
2212 };
2213
2214 inline SCHUR::SCHUR (const Matrix& a, const char *ord) { init (a, ord); }
2215
2216 inline SCHUR::SCHUR (const Matrix& a, const char *ord, int& info)
2217 { info = init (a, ord); }
2218
2219 inline SCHUR::SCHUR (const SCHUR& a, const char *ord)
2220 {
2221 schur_mat = a.schur_mat;
2222 unitary_mat = a.unitary_mat;
2223 }
2224
2225 inline SCHUR&
2226 SCHUR::operator = (const SCHUR& a, const char *ord)
2227 {
2228 schur_mat = a.schur_mat;
2229 unitary_mat = a.unitary_mat;
2230
2231 return *this;
2232 }
2233
2234 inline Matrix SCHUR::schur_matrix (void) const { return schur_mat; }
2235 inline Matrix SCHUR::unitary_matrix (void) const { return unitary_mat; }
2236
2237 /*
2238 * Result of a Schur Decomposition
2239 */
2240
2241 class ComplexSCHUR
2242 {
2243 friend class ComplexMatrix;
2244
2245 public:
2246 ComplexSCHUR (void) {}
2247
2248 ComplexSCHUR (const ComplexMatrix& a, const char *ord);
2249 ComplexSCHUR (const ComplexMatrix& a, const char *ord, int& info);
2250
2251 ComplexSCHUR (const ComplexSCHUR& a, const char *ord);
2252
2253 ComplexSCHUR& operator = (const ComplexSCHUR& a, const char *ord);
2254
2255 ComplexMatrix schur_matrix (void) const;
2256 ComplexMatrix unitary_matrix (void) const;
2257
2258 friend ostream& operator << (ostream& os, const ComplexSCHUR& a);
2259
2260 private:
2261 int init (const ComplexMatrix& a, const char *ord);
2262
2263 ComplexMatrix schur_mat;
2264 ComplexMatrix unitary_mat;
2265 };
2266
2267 inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord)
2268 { init (a,ord); }
2269
2270 inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord,
2271 int& info)
2272 { info = init (a,ord); }
2273
2274 inline ComplexSCHUR::ComplexSCHUR (const ComplexSCHUR& a, const char *ord)
2275 {
2276 schur_mat = a.schur_mat;
2277 unitary_mat = a.unitary_mat;
2278 }
2279
2280 inline ComplexSCHUR&
2281 ComplexSCHUR::operator = (const ComplexSCHUR& a, const char *ord)
2282 {
2283 schur_mat = a.schur_mat;
2284 unitary_mat = a.unitary_mat;
2285
2286 return *this;
2287 }
2288 inline ComplexMatrix ComplexSCHUR::schur_matrix (void) const
2289 { return schur_mat; }
2290
2291 inline ComplexMatrix ComplexSCHUR::unitary_matrix (void) const
2292 { return unitary_mat; }
2293
2294
2295 /*
2296 * Result of a Singular Value Decomposition.
2297 */
2298
2299 class SVD
2300 {
2301 friend class Matrix;
2302
2303 public:
2304 SVD (void) {}
2305
2306 SVD (const Matrix& a);
2307 SVD (const Matrix& a, int& info);
2308
2309 SVD (const SVD& a);
2310
2311 SVD& operator = (const SVD& a);
2312
2313 DiagMatrix singular_values (void) const;
2314 Matrix left_singular_matrix (void) const;
2315 Matrix right_singular_matrix (void) const;
2316
2317 friend ostream& operator << (ostream& os, const SVD& a);
2318
2319 private:
2320 int init (const Matrix& a);
2321
2322 DiagMatrix sigma;
2323 Matrix left_sm;
2324 Matrix right_sm;
2325 };
2326
2327 inline SVD::SVD (const Matrix& a) { init (a); }
2328
2329 inline SVD::SVD (const Matrix& a, int& info) { info = init (a); }
2330
2331 inline SVD::SVD (const SVD& a)
2332 {
2333 sigma = a.sigma;
2334 left_sm = a.left_sm;
2335 right_sm = a.right_sm;
2336 }
2337
2338 inline SVD&
2339 SVD::operator = (const SVD& a)
2340 {
2341 sigma = a.sigma;
2342 left_sm = a.left_sm;
2343 right_sm = a.right_sm;
2344
2345 return *this;
2346 }
2347
2348 inline DiagMatrix SVD::singular_values (void) const { return sigma; }
2349 inline Matrix SVD::left_singular_matrix (void) const { return left_sm; }
2350 inline Matrix SVD::right_singular_matrix (void) const { return right_sm; }
2351
2352 /*
2353 * Result of a Singular Value Decomposition.
2354 */
2355
2356 class ComplexSVD
2357 {
2358 friend class ComplexMatrix;
2359
2360 public:
2361 ComplexSVD (void) {}
2362
2363 ComplexSVD (const ComplexMatrix& a);
2364 ComplexSVD (const ComplexMatrix& a, int& info);
2365
2366 ComplexSVD (const ComplexSVD& a);
2367
2368 ComplexSVD& operator = (const ComplexSVD& a);
2369
2370 DiagMatrix singular_values (void) const;
2371 ComplexMatrix left_singular_matrix (void) const;
2372 ComplexMatrix right_singular_matrix (void) const;
2373
2374 friend ostream& operator << (ostream& os, const ComplexSVD& a);
2375
2376 private:
2377 int init (const ComplexMatrix& a);
2378
2379 DiagMatrix sigma;
2380 ComplexMatrix left_sm;
2381 ComplexMatrix right_sm;
2382 };
2383
2384 inline ComplexSVD::ComplexSVD (const ComplexMatrix& a) { init (a); }
2385 inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, int& info)
2386 { info = init (a); }
2387
2388 inline ComplexSVD::ComplexSVD (const ComplexSVD& a)
2389 {
2390 sigma = a.sigma;
2391 left_sm = a.left_sm;
2392 right_sm = a.right_sm;
2393 }
2394
2395 inline ComplexSVD&
2396 ComplexSVD::operator = (const ComplexSVD& a)
2397 {
2398 sigma = a.sigma;
2399 left_sm = a.left_sm;
2400 right_sm = a.right_sm;
2401
2402 return *this;
2403 }
2404
2405 inline DiagMatrix ComplexSVD::singular_values (void) const
2406 { return sigma; }
2407
2408 inline ComplexMatrix ComplexSVD::left_singular_matrix (void) const
2409 { return left_sm; }
2410
2411 inline ComplexMatrix ComplexSVD::right_singular_matrix (void) const
2412 { return right_sm; }
2413
2414 /*
2415 * Result of an Eigenvalue computation.
2416 */
2417
2418 class EIG
2419 {
2420 friend class Matrix;
2421 friend class ComplexMatrix;
2422
2423 public:
2424 EIG (void) {}
2425
2426 EIG (const Matrix& a);
2427 EIG (const Matrix& a, int& info);
2428
2429 EIG (const ComplexMatrix& a);
2430 EIG (const ComplexMatrix& a, int& info);
2431
2432 EIG (const EIG& a);
2433
2434 EIG& operator = (const EIG& a);
2435
2436 ComplexColumnVector eigenvalues (void) const;
2437 ComplexMatrix eigenvectors (void) const;
2438
2439 friend ostream& operator << (ostream& os, const EIG& a);
2440
2441 private:
2442 int init (const Matrix& a);
2443 int init (const ComplexMatrix& a);
2444
2445 ComplexColumnVector lambda;
2446 ComplexMatrix v;
2447 };
2448
2449 inline EIG::EIG (const Matrix& a) { init (a); }
2450 inline EIG::EIG (const Matrix& a, int& info) { info = init (a); }
2451
2452 inline EIG::EIG (const ComplexMatrix& a) { init (a); }
2453 inline EIG::EIG (const ComplexMatrix& a, int& info) { info = init (a); }
2454
2455 inline EIG::EIG (const EIG& a) { lambda = a.lambda; v = a.v; }
2456
2457 inline EIG& EIG::operator = (const EIG& a)
2458 { lambda = a.lambda; v = a.v; return *this; }
2459
2460 inline ComplexColumnVector EIG::eigenvalues (void) const { return lambda; }
2461
2462 inline ComplexMatrix EIG::eigenvectors (void) const { return v; }
2463
2464 /*
2465 * Result of an LU decomposition.
2466 */
2467
2468 class LU
2469 {
2470 friend class Matrix;
2471
2472 public:
2473 LU (void) {}
2474
2475 LU (const Matrix& a);
2476
2477 LU (const LU& a);
2478
2479 LU& operator = (const LU& a);
2480
2481 Matrix L (void) const;
2482 Matrix U (void) const;
2483 Matrix P (void) const;
2484
2485 friend ostream& operator << (ostream& os, const LU& a);
2486
2487 private:
2488
2489 Matrix l;
2490 Matrix u;
2491 Matrix p;
2492 };
2493
2494 inline LU::LU (const LU& a) { l = a.l; u = a.u; p = a.p; }
2495
2496 inline LU& LU::operator = (const LU& a)
2497 { l = a.l; u = a.u; p = a.p; return *this; }
2498
2499 inline Matrix LU::L (void) const { return l; }
2500 inline Matrix LU::U (void) const { return u; }
2501 inline Matrix LU::P (void) const { return p; }
2502
2503 class ComplexLU
2504 {
2505 friend class ComplexMatrix;
2506
2507 public:
2508 ComplexLU (void) {}
2509
2510 ComplexLU (const ComplexMatrix& a);
2511
2512 ComplexLU (const ComplexLU& a);
2513
2514 ComplexLU& operator = (const ComplexLU& a);
2515
2516 ComplexMatrix L (void) const;
2517 ComplexMatrix U (void) const;
2518 Matrix P (void) const;
2519
2520 friend ostream& operator << (ostream& os, const ComplexLU& a);
2521
2522 private:
2523
2524 ComplexMatrix l;
2525 ComplexMatrix u;
2526 Matrix p;
2527 };
2528
2529 inline ComplexLU::ComplexLU (const ComplexLU& a) { l = a.l; u = a.u; p = a.p; }
2530
2531 inline ComplexLU& ComplexLU::operator = (const ComplexLU& a)
2532 { l = a.l; u = a.u; p = a.p; return *this; }
2533
2534 inline ComplexMatrix ComplexLU::L (void) const { return l; }
2535 inline ComplexMatrix ComplexLU::U (void) const { return u; }
2536 inline Matrix ComplexLU::P (void) const { return p; }
2537
2538 /*
2539 * Result of a QR decomposition.
2540 */
2541
2542 class QR
2543 {
2544 public:
2545 QR (void) {}
2546
2547 QR (const Matrix& A);
2548
2549 QR (const QR& a);
2550
2551 QR& operator = (const QR& a);
2552
2553 Matrix Q (void) const;
2554 Matrix R (void) const;
2555
2556 friend ostream& operator << (ostream& os, const QR& a);
2557
2558 private:
2559 Matrix q;
2560 Matrix r;
2561 };
2562
2563 inline QR::QR (const QR& a) { q = a.q; r = a.r; }
2564
2565 inline QR& QR::operator = (const QR& a) { q = a.q; r = a.r; return *this; }
2566
2567 inline Matrix QR::Q (void) const { return q; }
2568 inline Matrix QR::R (void) const { return r; }
2569
2570 class ComplexQR
2571 {
2572 public:
2573 ComplexQR (void) {}
2574
2575 ComplexQR (const ComplexMatrix& A);
2576
2577 ComplexQR (const ComplexQR& a);
2578
2579 ComplexQR& operator = (const ComplexQR& a);
2580
2581 ComplexMatrix Q (void) const;
2582 ComplexMatrix R (void) const;
2583
2584 friend ostream& operator << (ostream& os, const ComplexQR& a);
2585
2586 private:
2587 ComplexMatrix q;
2588 ComplexMatrix r;
2589 };
2590
2591 inline ComplexQR::ComplexQR (const ComplexQR& a) { q = a.q; r = a.r; }
2592
2593 inline ComplexQR& ComplexQR::operator = (const ComplexQR& a)
2594 { q = a.q; r = a.r; return *this; }
2595
2596 inline ComplexMatrix ComplexQR::Q (void) const { return q; }
2597 inline ComplexMatrix ComplexQR::R (void) const { return r; }
2598
2599 #endif
2600
2601 /*
2602 ;;; Local Variables: ***
2603 ;;; mode: C++ ***
2604 ;;; page-delimiter: "^/\\*" ***
2605 ;;; End: ***
2606 */