Mercurial > hg > octave-lyh

comparison liboctave/chMatrix.cc @ 1573:403c60daa8c7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-10-19 04:34:20 by jwe] Initial revision
author jwe
date Thu, 19 Oct 1995 04:34:20 +0000
parents
children 42b4f904f1af
comparison
equal deleted inserted replaced
1572:0d9e10d10bd7 1573:403c60daa8c7
1 // Matrix manipulations. -*- C++ -*-
2 /*
3
4 Copyright (C) 1992, 1993, 1994, 1995 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, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21
22 */
23
24 #if defined (__GNUG__)
25 #pragma implementation
26 #endif
27
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31
32 #include <cstdio>
33 #include <cstring>
34
35 #include <iostream.h>
36
37 // #include <sys/types.h> // XXX FIXME XXX
38
39 #include "lo-error.h"
40 #include "mx-base.h"
41 #include "mx-inlines.cc"
42
43 // charMatrix class.
44
45 charMatrix::charMatrix (const char *s)
46 : MArray2<char> ((s ? 1 : 0), (s ? strlen (s) : 0))
47 {
48 int nc = cols ();
49 for (int i = 0; i < nc; i++)
50 elem (0, i) = s[i];
51 }
52
53 int
54 charMatrix::operator == (const charMatrix& a) const
55 {
56 if (rows () != a.rows () || cols () != a.cols ())
57 return 0;
58
59 return equal (data (), a.data (), length ());
60 }
61
62 int
63 charMatrix::operator != (const charMatrix& a) const
64 {
65 return !(*this == a);
66 }
67
68 charMatrix&
69 charMatrix::insert (const char *s, int r, int c)
70 {
71 if (s)
72 {
73 int s_len = strlen (s);
74
75 if (r < 0 || r >= rows () || c < 0 || c + s_len - 1 > cols ())
76 {
77 (*current_liboctave_error_handler) ("range error for insert");
78 return *this;
79 }
80
81 for (int i = 0; i < s_len; i++)
82 elem (r, c+i) = s[i];
83 }
84 return *this;
85 }
86
87 charMatrix&
88 charMatrix::insert (const charMatrix& a, int r, int c)
89 {
90 Array2<char>::insert (a, r, c);
91 return *this;
92 }
93
94 char *
95 charMatrix::row_as_string (int r) const
96 {
97 if (r < 0 || r >= rows ())
98 {
99 (*current_liboctave_error_handler) ("range error for row_as_string");
100 return 0;
101 }
102
103 int nc = cols ();
104
105 char *retval = new char [nc + 1];
106
107 retval[nc] = '\0';
108
109 for (int i = 0; i < nc; i++)
110 retval[i] = elem (r, i);
111
112 return retval;
113 }
114
115 #if 0
116 Matrix&
117 Matrix::insert (const RowVector& a, int r, int c)
118 {
119 int a_len = a.length ();
120 if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ())
121 {
122 (*current_liboctave_error_handler) ("range error for insert");
123 return *this;
124 }
125
126 for (int i = 0; i < a_len; i++)
127 elem (r, c+i) = a.elem (i);
128
129 return *this;
130 }
131
132 Matrix&
133 Matrix::insert (const ColumnVector& a, int r, int c)
134 {
135 int a_len = a.length ();
136 if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ())
137 {
138 (*current_liboctave_error_handler) ("range error for insert");
139 return *this;
140 }
141
142 for (int i = 0; i < a_len; i++)
143 elem (r+i, c) = a.elem (i);
144
145 return *this;
146 }
147
148 Matrix&
149 Matrix::insert (const DiagMatrix& a, int r, int c)
150 {
151 if (r < 0 || r + a.rows () - 1 > rows ()
152 || c < 0 || c + a.cols () - 1 > cols ())
153 {
154 (*current_liboctave_error_handler) ("range error for insert");
155 return *this;
156 }
157
158 for (int i = 0; i < a.length (); i++)
159 elem (r+i, c+i) = a.elem (i, i);
160
161 return *this;
162 }
163
164 Matrix&
165 Matrix::fill (double val)
166 {
167 int nr = rows ();
168 int nc = cols ();
169 if (nr > 0 && nc > 0)
170 for (int j = 0; j < nc; j++)
171 for (int i = 0; i < nr; i++)
172 elem (i, j) = val;
173
174 return *this;
175 }
176
177 Matrix&
178 Matrix::fill (double val, int r1, int c1, int r2, int c2)
179 {
180 int nr = rows ();
181 int nc = cols ();
182 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
183 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
184 {
185 (*current_liboctave_error_handler) ("range error for fill");
186 return *this;
187 }
188
189 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
190 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
191
192 for (int j = c1; j <= c2; j++)
193 for (int i = r1; i <= r2; i++)
194 elem (i, j) = val;
195
196 return *this;
197 }
198
199 Matrix
200 Matrix::append (const Matrix& a) const
201 {
202 int nr = rows ();
203 int nc = cols ();
204 if (nr != a.rows ())
205 {
206 (*current_liboctave_error_handler) ("row dimension mismatch for append");
207 return Matrix ();
208 }
209
210 int nc_insert = nc;
211 Matrix retval (nr, nc + a.cols ());
212 retval.insert (*this, 0, 0);
213 retval.insert (a, 0, nc_insert);
214 return retval;
215 }
216
217 Matrix
218 Matrix::append (const RowVector& a) const
219 {
220 int nr = rows ();
221 int nc = cols ();
222 if (nr != 1)
223 {
224 (*current_liboctave_error_handler) ("row dimension mismatch for append");
225 return Matrix ();
226 }
227
228 int nc_insert = nc;
229 Matrix retval (nr, nc + a.length ());
230 retval.insert (*this, 0, 0);
231 retval.insert (a, 0, nc_insert);
232 return retval;
233 }
234
235 Matrix
236 Matrix::append (const ColumnVector& a) const
237 {
238 int nr = rows ();
239 int nc = cols ();
240 if (nr != a.length ())
241 {
242 (*current_liboctave_error_handler) ("row dimension mismatch for append");
243 return Matrix ();
244 }
245
246 int nc_insert = nc;
247 Matrix retval (nr, nc + 1);
248 retval.insert (*this, 0, 0);
249 retval.insert (a, 0, nc_insert);
250 return retval;
251 }
252
253 Matrix
254 Matrix::append (const DiagMatrix& a) const
255 {
256 int nr = rows ();
257 int nc = cols ();
258 if (nr != a.rows ())
259 {
260 (*current_liboctave_error_handler) ("row dimension mismatch for append");
261 return *this;
262 }
263
264 int nc_insert = nc;
265 Matrix retval (nr, nc + a.cols ());
266 retval.insert (*this, 0, 0);
267 retval.insert (a, 0, nc_insert);
268 return retval;
269 }
270
271 Matrix
272 Matrix::stack (const Matrix& a) const
273 {
274 int nr = rows ();
275 int nc = cols ();
276 if (nc != a.cols ())
277 {
278 (*current_liboctave_error_handler)
279 ("column dimension mismatch for stack");
280 return Matrix ();
281 }
282
283 int nr_insert = nr;
284 Matrix retval (nr + a.rows (), nc);
285 retval.insert (*this, 0, 0);
286 retval.insert (a, nr_insert, 0);
287 return retval;
288 }
289
290 Matrix
291 Matrix::stack (const RowVector& a) const
292 {
293 int nr = rows ();
294 int nc = cols ();
295 if (nc != a.length ())
296 {
297 (*current_liboctave_error_handler)
298 ("column dimension mismatch for stack");
299 return Matrix ();
300 }
301
302 int nr_insert = nr;
303 Matrix retval (nr + 1, nc);
304 retval.insert (*this, 0, 0);
305 retval.insert (a, nr_insert, 0);
306 return retval;
307 }
308
309 Matrix
310 Matrix::stack (const ColumnVector& a) const
311 {
312 int nr = rows ();
313 int nc = cols ();
314 if (nc != 1)
315 {
316 (*current_liboctave_error_handler)
317 ("column dimension mismatch for stack");
318 return Matrix ();
319 }
320
321 int nr_insert = nr;
322 Matrix retval (nr + a.length (), nc);
323 retval.insert (*this, 0, 0);
324 retval.insert (a, nr_insert, 0);
325 return retval;
326 }
327
328 Matrix
329 Matrix::stack (const DiagMatrix& a) const
330 {
331 int nr = rows ();
332 int nc = cols ();
333 if (nc != a.cols ())
334 {
335 (*current_liboctave_error_handler)
336 ("column dimension mismatch for stack");
337 return Matrix ();
338 }
339
340 int nr_insert = nr;
341 Matrix retval (nr + a.rows (), nc);
342 retval.insert (*this, 0, 0);
343 retval.insert (a, nr_insert, 0);
344 return retval;
345 }
346
347 Matrix
348 Matrix::transpose (void) const
349 {
350 int nr = rows ();
351 int nc = cols ();
352 Matrix result (nc, nr);
353 if (length () > 0)
354 {
355 for (int j = 0; j < nc; j++)
356 for (int i = 0; i < nr; i++)
357 result.elem (j, i) = elem (i, j);
358 }
359 return result;
360 }
361
362 Matrix
363 real (const ComplexMatrix& a)
364 {
365 int a_len = a.length ();
366 Matrix retval;
367 if (a_len > 0)
368 retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ());
369 return retval;
370 }
371
372 Matrix
373 imag (const ComplexMatrix& a)
374 {
375 int a_len = a.length ();
376 Matrix retval;
377 if (a_len > 0)
378 retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ());
379 return retval;
380 }
381
382 Matrix
383 Matrix::extract (int r1, int c1, int r2, int c2) const
384 {
385 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
386 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
387
388 int new_r = r2 - r1 + 1;
389 int new_c = c2 - c1 + 1;
390
391 Matrix result (new_r, new_c);
392
393 for (int j = 0; j < new_c; j++)
394 for (int i = 0; i < new_r; i++)
395 result.elem (i, j) = elem (r1+i, c1+j);
396
397 return result;
398 }
399
400 // extract row or column i.
401
402 RowVector
403 Matrix::row (int i) const
404 {
405 int nc = cols ();
406 if (i < 0 || i >= rows ())
407 {
408 (*current_liboctave_error_handler) ("invalid row selection");
409 return RowVector ();
410 }
411
412 RowVector retval (nc);
413 for (int j = 0; j < nc; j++)
414 retval.elem (j) = elem (i, j);
415
416 return retval;
417 }
418
419 RowVector
420 Matrix::row (char *s) const
421 {
422 if (! s)
423 {
424 (*current_liboctave_error_handler) ("invalid row selection");
425 return RowVector ();
426 }
427
428 char c = *s;
429 if (c == 'f' || c == 'F')
430 return row (0);
431 else if (c == 'l' || c == 'L')
432 return row (rows () - 1);
433 else
434 {
435 (*current_liboctave_error_handler) ("invalid row selection");
436 return RowVector ();
437 }
438 }
439
440 ColumnVector
441 Matrix::column (int i) const
442 {
443 int nr = rows ();
444 if (i < 0 || i >= cols ())
445 {
446 (*current_liboctave_error_handler) ("invalid column selection");
447 return ColumnVector ();
448 }
449
450 ColumnVector retval (nr);
451 for (int j = 0; j < nr; j++)
452 retval.elem (j) = elem (j, i);
453
454 return retval;
455 }
456
457 ColumnVector
458 Matrix::column (char *s) const
459 {
460 if (! s)
461 {
462 (*current_liboctave_error_handler) ("invalid column selection");
463 return ColumnVector ();
464 }
465
466 char c = *s;
467 if (c == 'f' || c == 'F')
468 return column (0);
469 else if (c == 'l' || c == 'L')
470 return column (cols () - 1);
471 else
472 {
473 (*current_liboctave_error_handler) ("invalid column selection");
474 return ColumnVector ();
475 }
476 }
477
478 Matrix
479 Matrix::inverse (void) const
480 {
481 int info;
482 double rcond;
483 return inverse (info, rcond);
484 }
485
486 Matrix
487 Matrix::inverse (int& info) const
488 {
489 double rcond;
490 return inverse (info, rcond);
491 }
492
493 Matrix
494 Matrix::inverse (int& info, double& rcond) const
495 {
496 int nr = rows ();
497 int nc = cols ();
498 int len = length ();
499 if (nr != nc || nr == 0 || nc == 0)
500 {
501 (*current_liboctave_error_handler) ("inverse requires square matrix");
502 return Matrix ();
503 }
504
505 info = 0;
506
507 int *ipvt = new int [nr];
508 double *z = new double [nr];
509 double *tmp_data = dup (data (), len);
510
511 F77_FCN (dgeco, DGECO) (tmp_data, nr, nc, ipvt, rcond, z);
512
513 volatile double rcond_plus_one = rcond + 1.0;
514 if (rcond_plus_one == 1.0)
515 {
516 info = -1;
517 copy (tmp_data, data (), len); // Restore matrix contents.
518 }
519 else
520 {
521 double *dummy = 0;
522
523 F77_FCN (dgedi, DGEDI) (tmp_data, nr, nc, ipvt, dummy, z, 1);
524 }
525
526 delete [] ipvt;
527 delete [] z;
528
529 return Matrix (tmp_data, nr, nc);
530 }
531
532 Matrix
533 Matrix::pseudo_inverse (double tol)
534 {
535 SVD result (*this);
536
537 DiagMatrix S = result.singular_values ();
538 Matrix U = result.left_singular_matrix ();
539 Matrix V = result.right_singular_matrix ();
540
541 ColumnVector sigma = S.diag ();
542
543 int r = sigma.length () - 1;
544 int nr = rows ();
545 int nc = cols ();
546
547 if (tol <= 0.0)
548 {
549 if (nr > nc)
550 tol = nr * sigma.elem (0) * DBL_EPSILON;
551 else
552 tol = nc * sigma.elem (0) * DBL_EPSILON;
553 }
554
555 while (r >= 0 && sigma.elem (r) < tol)
556 r--;
557
558 if (r < 0)
559 return Matrix (nc, nr, 0.0);
560 else
561 {
562 Matrix Ur = U.extract (0, 0, nr-1, r);
563 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
564 Matrix Vr = V.extract (0, 0, nc-1, r);
565 return Vr * D * Ur.transpose ();
566 }
567 }
568
569 ComplexMatrix
570 Matrix::fourier (void) const
571 {
572 int nr = rows ();
573 int nc = cols ();
574 int npts, nsamples;
575 if (nr == 1 || nc == 1)
576 {
577 npts = nr > nc ? nr : nc;
578 nsamples = 1;
579 }
580 else
581 {
582 npts = nr;
583 nsamples = nc;
584 }
585
586 int nn = 4*npts+15;
587 Complex *wsave = new Complex [nn];
588 Complex *tmp_data = make_complex (data (), length ());
589
590 F77_FCN (cffti, CFFTI) (npts, wsave);
591
592 for (int j = 0; j < nsamples; j++)
593 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], wsave);
594
595 delete [] wsave;
596
597 return ComplexMatrix (tmp_data, nr, nc);
598 }
599
600 ComplexMatrix
601 Matrix::ifourier (void) const
602 {
603 int nr = rows ();
604 int nc = cols ();
605 int npts, nsamples;
606 if (nr == 1 || nc == 1)
607 {
608 npts = nr > nc ? nr : nc;
609 nsamples = 1;
610 }
611 else
612 {
613 npts = nr;
614 nsamples = nc;
615 }
616
617 int nn = 4*npts+15;
618 Complex *wsave = new Complex [nn];
619 Complex *tmp_data = make_complex (data (), length ());
620
621 F77_FCN (cffti, CFFTI) (npts, wsave);
622
623 for (int j = 0; j < nsamples; j++)
624 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], wsave);
625
626 for (int j = 0; j < npts*nsamples; j++)
627 tmp_data[j] = tmp_data[j] / (double) npts;
628
629 delete [] wsave;
630
631 return ComplexMatrix (tmp_data, nr, nc);
632 }
633
634 ComplexMatrix
635 Matrix::fourier2d (void) const
636 {
637 int nr = rows ();
638 int nc = cols ();
639 int npts, nsamples;
640 if (nr == 1 || nc == 1)
641 {
642 npts = nr > nc ? nr : nc;
643 nsamples = 1;
644 }
645 else
646 {
647 npts = nr;
648 nsamples = nc;
649 }
650
651 int nn = 4*npts+15;
652 Complex *wsave = new Complex [nn];
653 Complex *tmp_data = make_complex (data (), length ());
654
655 F77_FCN (cffti, CFFTI) (npts, wsave);
656
657 for (int j = 0; j < nsamples; j++)
658 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], wsave);
659
660 delete [] wsave;
661
662 npts = nc;
663 nsamples = nr;
664 nn = 4*npts+15;
665 wsave = new Complex [nn];
666 Complex *row = new Complex[npts];
667
668 F77_FCN (cffti, CFFTI) (npts, wsave);
669
670 for (int j = 0; j < nsamples; j++)
671 {
672 for (int i = 0; i < npts; i++)
673 row[i] = tmp_data[i*nr + j];
674
675 F77_FCN (cfftf, CFFTF) (npts, row, wsave);
676
677 for (int i = 0; i < npts; i++)
678 tmp_data[i*nr + j] = row[i];
679 }
680
681 delete [] wsave;
682 delete [] row;
683
684 return ComplexMatrix (tmp_data, nr, nc);
685 }
686
687 ComplexMatrix
688 Matrix::ifourier2d (void) const
689 {
690 int nr = rows ();
691 int nc = cols ();
692 int npts, nsamples;
693 if (nr == 1 || nc == 1)
694 {
695 npts = nr > nc ? nr : nc;
696 nsamples = 1;
697 }
698 else
699 {
700 npts = nr;
701 nsamples = nc;
702 }
703
704 int nn = 4*npts+15;
705 Complex *wsave = new Complex [nn];
706 Complex *tmp_data = make_complex (data (), length ());
707
708 F77_FCN (cffti, CFFTI) (npts, wsave);
709
710 for (int j = 0; j < nsamples; j++)
711 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], wsave);
712
713 delete [] wsave;
714
715 for (int j = 0; j < npts*nsamples; j++)
716 tmp_data[j] = tmp_data[j] / (double) npts;
717
718 npts = nc;
719 nsamples = nr;
720 nn = 4*npts+15;
721 wsave = new Complex [nn];
722 Complex *row = new Complex[npts];
723
724 F77_FCN (cffti, CFFTI) (npts, wsave);
725
726 for (int j = 0; j < nsamples; j++)
727 {
728 for (int i = 0; i < npts; i++)
729 row[i] = tmp_data[i*nr + j];
730
731 F77_FCN (cfftb, CFFTB) (npts, row, wsave);
732
733 for (int i = 0; i < npts; i++)
734 tmp_data[i*nr + j] = row[i] / (double) npts;
735 }
736
737 delete [] wsave;
738 delete [] row;
739
740 return ComplexMatrix (tmp_data, nr, nc);
741 }
742
743 DET
744 Matrix::determinant (void) const
745 {
746 int info;
747 double rcond;
748 return determinant (info, rcond);
749 }
750
751 DET
752 Matrix::determinant (int& info) const
753 {
754 double rcond;
755 return determinant (info, rcond);
756 }
757
758 DET
759 Matrix::determinant (int& info, double& rcond) const
760 {
761 DET retval;
762
763 int nr = rows ();
764 int nc = cols ();
765
766 if (nr == 0 || nc == 0)
767 {
768 double d[2];
769 d[0] = 1.0;
770 d[1] = 0.0;
771 retval = DET (d);
772 }
773 else
774 {
775 info = 0;
776 int *ipvt = new int [nr];
777
778 double *z = new double [nr];
779 double *tmp_data = dup (data (), length ());
780
781 F77_FCN (dgeco, DGECO) (tmp_data, nr, nr, ipvt, rcond, z);
782
783 volatile double rcond_plus_one = rcond + 1.0;
784 if (rcond_plus_one == 1.0)
785 {
786 info = -1;
787 retval = DET ();
788 }
789 else
790 {
791 double d[2];
792 F77_FCN (dgedi, DGEDI) (tmp_data, nr, nr, ipvt, d, z, 10);
793 retval = DET (d);
794 }
795
796 delete [] tmp_data;
797 delete [] ipvt;
798 delete [] z;
799 }
800
801 return retval;
802 }
803
804 Matrix
805 Matrix::solve (const Matrix& b) const
806 {
807 int info;
808 double rcond;
809 return solve (b, info, rcond);
810 }
811
812 Matrix
813 Matrix::solve (const Matrix& b, int& info) const
814 {
815 double rcond;
816 return solve (b, info, rcond);
817 }
818
819 Matrix
820 Matrix::solve (const Matrix& b, int& info, double& rcond) const
821 {
822 Matrix retval;
823
824 int nr = rows ();
825 int nc = cols ();
826 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
827 {
828 (*current_liboctave_error_handler)
829 ("matrix dimension mismatch solution of linear equations");
830 return Matrix ();
831 }
832
833 info = 0;
834 int *ipvt = new int [nr];
835
836 double *z = new double [nr];
837 double *tmp_data = dup (data (), length ());
838
839 F77_FCN (dgeco, DGECO) (tmp_data, nr, nr, ipvt, rcond, z);
840
841 volatile double rcond_plus_one = rcond + 1.0;
842 if (rcond_plus_one == 1.0)
843 {
844 info = -2;
845 }
846 else
847 {
848 double *result = dup (b.data (), b.length ());
849
850 int b_nc = b.cols ();
851 for (int j = 0; j < b_nc; j++)
852 F77_FCN (dgesl, DGESL) (tmp_data, nr, nr, ipvt, &result[nr*j], 0);
853
854 retval = Matrix (result, b.rows (), b_nc);
855 }
856
857 delete [] tmp_data;
858 delete [] ipvt;
859 delete [] z;
860
861 return retval;
862 }
863
864 ComplexMatrix
865 Matrix::solve (const ComplexMatrix& b) const
866 {
867 ComplexMatrix tmp (*this);
868 return tmp.solve (b);
869 }
870
871 ComplexMatrix
872 Matrix::solve (const ComplexMatrix& b, int& info) const
873 {
874 ComplexMatrix tmp (*this);
875 return tmp.solve (b, info);
876 }
877
878 ComplexMatrix
879 Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
880 {
881 ComplexMatrix tmp (*this);
882 return tmp.solve (b, info, rcond);
883 }
884
885 ColumnVector
886 Matrix::solve (const ColumnVector& b) const
887 {
888 int info; double rcond;
889 return solve (b, info, rcond);
890 }
891
892 ColumnVector
893 Matrix::solve (const ColumnVector& b, int& info) const
894 {
895 double rcond;
896 return solve (b, info, rcond);
897 }
898
899 ColumnVector
900 Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
901 {
902 ColumnVector retval;
903
904 int nr = rows ();
905 int nc = cols ();
906 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
907 {
908 (*current_liboctave_error_handler)
909 ("matrix dimension mismatch solution of linear equations");
910 return ColumnVector ();
911 }
912
913 info = 0;
914 int *ipvt = new int [nr];
915
916 double *z = new double [nr];
917 double *tmp_data = dup (data (), length ());
918
919 F77_FCN (dgeco, DGECO) (tmp_data, nr, nr, ipvt, rcond, z);
920
921 volatile double rcond_plus_one = rcond + 1.0;
922 if (rcond_plus_one == 1.0)
923 {
924 info = -2;
925 }
926 else
927 {
928 int b_len = b.length ();
929
930 double *result = dup (b.data (), b_len);
931
932 F77_FCN (dgesl, DGESL) (tmp_data, nr, nr, ipvt, result, 0);
933
934 retval = ColumnVector (result, b_len);
935 }
936
937 delete [] tmp_data;
938 delete [] ipvt;
939 delete [] z;
940
941 return retval;
942 }
943
944 ComplexColumnVector
945 Matrix::solve (const ComplexColumnVector& b) const
946 {
947 ComplexMatrix tmp (*this);
948 return tmp.solve (b);
949 }
950
951 ComplexColumnVector
952 Matrix::solve (const ComplexColumnVector& b, int& info) const
953 {
954 ComplexMatrix tmp (*this);
955 return tmp.solve (b, info);
956 }
957
958 ComplexColumnVector
959 Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
960 {
961 ComplexMatrix tmp (*this);
962 return tmp.solve (b, info, rcond);
963 }
964
965 Matrix
966 Matrix::lssolve (const Matrix& b) const
967 {
968 int info;
969 int rank;
970 return lssolve (b, info, rank);
971 }
972
973 Matrix
974 Matrix::lssolve (const Matrix& b, int& info) const
975 {
976 int rank;
977 return lssolve (b, info, rank);
978 }
979
980 Matrix
981 Matrix::lssolve (const Matrix& b, int& info, int& rank) const
982 {
983 int nrhs = b.cols ();
984
985 int m = rows ();
986 int n = cols ();
987
988 if (m == 0 || n == 0 || m != b.rows ())
989 {
990 (*current_liboctave_error_handler)
991 ("matrix dimension mismatch in solution of least squares problem");
992 return Matrix ();
993 }
994
995 double *tmp_data = dup (data (), length ());
996
997 int nrr = m > n ? m : n;
998 Matrix result (nrr, nrhs);
999
1000 for (int j = 0; j < nrhs; j++)
1001 for (int i = 0; i < m; i++)
1002 result.elem (i, j) = b.elem (i, j);
1003
1004 double *presult = result.fortran_vec ();
1005
1006 int len_s = m < n ? m : n;
1007 double *s = new double [len_s];
1008 double rcond = -1.0;
1009 int lwork;
1010 if (m < n)
1011 lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
1012 else
1013 lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
1014
1015 double *work = new double [lwork];
1016
1017 F77_FCN (dgelss, DGELSS) (m, n, nrhs, tmp_data, m, presult, nrr, s,
1018 rcond, rank, work, lwork, info);
1019
1020 Matrix retval (n, nrhs);
1021 for (int j = 0; j < nrhs; j++)
1022 for (int i = 0; i < n; i++)
1023 retval.elem (i, j) = result.elem (i, j);
1024
1025 delete [] tmp_data;
1026 delete [] s;
1027 delete [] work;
1028
1029 return retval;
1030 }
1031
1032 ComplexMatrix
1033 Matrix::lssolve (const ComplexMatrix& b) const
1034 {
1035 ComplexMatrix tmp (*this);
1036 int info;
1037 int rank;
1038 return tmp.lssolve (b, info, rank);
1039 }
1040
1041 ComplexMatrix
1042 Matrix::lssolve (const ComplexMatrix& b, int& info) const
1043 {
1044 ComplexMatrix tmp (*this);
1045 int rank;
1046 return tmp.lssolve (b, info, rank);
1047 }
1048
1049 ComplexMatrix
1050 Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
1051 {
1052 ComplexMatrix tmp (*this);
1053 return tmp.lssolve (b, info, rank);
1054 }
1055
1056 ColumnVector
1057 Matrix::lssolve (const ColumnVector& b) const
1058 {
1059 int info;
1060 int rank;
1061 return lssolve (b, info, rank);
1062 }
1063
1064 ColumnVector
1065 Matrix::lssolve (const ColumnVector& b, int& info) const
1066 {
1067 int rank;
1068 return lssolve (b, info, rank);
1069 }
1070
1071 ColumnVector
1072 Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const
1073 {
1074 int nrhs = 1;
1075
1076 int m = rows ();
1077 int n = cols ();
1078
1079 if (m == 0 || n == 0 || m != b.length ())
1080 {
1081 (*current_liboctave_error_handler)
1082 ("matrix dimension mismatch in solution of least squares problem");
1083 return ColumnVector ();
1084 }
1085
1086 double *tmp_data = dup (data (), length ());
1087
1088 int nrr = m > n ? m : n;
1089 ColumnVector result (nrr);
1090
1091 for (int i = 0; i < m; i++)
1092 result.elem (i) = b.elem (i);
1093
1094 double *presult = result.fortran_vec ();
1095
1096 int len_s = m < n ? m : n;
1097 double *s = new double [len_s];
1098 double rcond = -1.0;
1099 int lwork;
1100 if (m < n)
1101 lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
1102 else
1103 lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
1104
1105 double *work = new double [lwork];
1106
1107 F77_FCN (dgelss, DGELSS) (m, n, nrhs, tmp_data, m, presult, nrr, s,
1108 rcond, rank, work, lwork, info);
1109
1110 ColumnVector retval (n);
1111 for (int i = 0; i < n; i++)
1112 retval.elem (i) = result.elem (i);
1113
1114 delete [] tmp_data;
1115 delete [] s;
1116 delete [] work;
1117
1118 return retval;
1119 }
1120
1121 ComplexColumnVector
1122 Matrix::lssolve (const ComplexColumnVector& b) const
1123 {
1124 ComplexMatrix tmp (*this);
1125 return tmp.lssolve (b);
1126 }
1127
1128 ComplexColumnVector
1129 Matrix::lssolve (const ComplexColumnVector& b, int& info) const
1130 {
1131 ComplexMatrix tmp (*this);
1132 return tmp.lssolve (b, info);
1133 }
1134
1135 ComplexColumnVector
1136 Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
1137 {
1138 ComplexMatrix tmp (*this);
1139 return tmp.lssolve (b, info, rank);
1140 }
1141
1142 Matrix&
1143 Matrix::operator += (const Matrix& a)
1144 {
1145 int nr = rows ();
1146 int nc = cols ();
1147 if (nr != a.rows () || nc != a.cols ())
1148 {
1149 (*current_liboctave_error_handler)
1150 ("nonconformant matrix += operation attempted");
1151 return *this;
1152 }
1153
1154 if (nr == 0 || nc == 0)
1155 return *this;
1156
1157 double *d = fortran_vec (); // Ensures only one reference to my privates!
1158
1159 add2 (d, a.data (), length ());
1160
1161 return *this;
1162 }
1163
1164 Matrix&
1165 Matrix::operator -= (const Matrix& a)
1166 {
1167 int nr = rows ();
1168 int nc = cols ();
1169 if (nr != a.rows () || nc != a.cols ())
1170 {
1171 (*current_liboctave_error_handler)
1172 ("nonconformant matrix -= operation attempted");
1173 return *this;
1174 }
1175
1176 if (nr == 0 || nc == 0)
1177 return *this;
1178
1179 double *d = fortran_vec (); // Ensures only one reference to my privates!
1180
1181 subtract2 (d, a.data (), length ());
1182
1183 return *this;
1184 }
1185
1186 Matrix&
1187 Matrix::operator += (const DiagMatrix& a)
1188 {
1189 if (rows () != a.rows () || cols () != a.cols ())
1190 {
1191 (*current_liboctave_error_handler)
1192 ("nonconformant matrix += operation attempted");
1193 return *this;
1194 }
1195
1196 for (int i = 0; i < a.length (); i++)
1197 elem (i, i) += a.elem (i, i);
1198
1199 return *this;
1200 }
1201
1202 Matrix&
1203 Matrix::operator -= (const DiagMatrix& a)
1204 {
1205 if (rows () != a.rows () || cols () != a.cols ())
1206 {
1207 (*current_liboctave_error_handler)
1208 ("nonconformant matrix += operation attempted");
1209 return *this;
1210 }
1211
1212 for (int i = 0; i < a.length (); i++)
1213 elem (i, i) -= a.elem (i, i);
1214
1215 return *this;
1216 }
1217
1218 // unary operations
1219
1220 Matrix
1221 Matrix::operator ! (void) const
1222 {
1223 int nr = rows ();
1224 int nc = cols ();
1225
1226 Matrix b (nr, nc);
1227
1228 for (int j = 0; j < nc; j++)
1229 for (int i = 0; i < nr; i++)
1230 b.elem (i, j) = ! elem (i, j);
1231
1232 return b;
1233 }
1234
1235 // column vector by row vector -> matrix operations
1236
1237 Matrix
1238 operator * (const ColumnVector& v, const RowVector& a)
1239 {
1240 int len = v.length ();
1241 int a_len = a.length ();
1242 if (len != a_len)
1243 {
1244 (*current_liboctave_error_handler)
1245 ("nonconformant vector multiplication attempted");
1246 return Matrix ();
1247 }
1248
1249 if (len == 0)
1250 return Matrix (len, len, 0.0);
1251
1252 double *c = new double [len * a_len];
1253
1254 F77_FCN (dgemm, DGEMM) ("N", "N", len, a_len, 1, 1.0, v.data (),
1255 len, a.data (), 1, 0.0, c, len, 1L, 1L);
1256
1257 return Matrix (c, len, a_len);
1258 }
1259
1260 // diagonal matrix by scalar -> matrix operations
1261
1262 Matrix
1263 operator + (const DiagMatrix& a, double s)
1264 {
1265 Matrix tmp (a.rows (), a.cols (), s);
1266 return a + tmp;
1267 }
1268
1269 Matrix
1270 operator - (const DiagMatrix& a, double s)
1271 {
1272 Matrix tmp (a.rows (), a.cols (), -s);
1273 return a + tmp;
1274 }
1275
1276 // scalar by diagonal matrix -> matrix operations
1277
1278 Matrix
1279 operator + (double s, const DiagMatrix& a)
1280 {
1281 Matrix tmp (a.rows (), a.cols (), s);
1282 return tmp + a;
1283 }
1284
1285 Matrix
1286 operator - (double s, const DiagMatrix& a)
1287 {
1288 Matrix tmp (a.rows (), a.cols (), s);
1289 return tmp - a;
1290 }
1291
1292 // matrix by diagonal matrix -> matrix operations
1293
1294 Matrix
1295 operator + (const Matrix& m, const DiagMatrix& a)
1296 {
1297 int nr = m.rows ();
1298 int nc = m.cols ();
1299 if (nr != a.rows () || nc != a.cols ())
1300 {
1301 (*current_liboctave_error_handler)
1302 ("nonconformant matrix addition attempted");
1303 return Matrix ();
1304 }
1305
1306 if (nr == 0 || nc == 0)
1307 return Matrix (nr, nc);
1308
1309 Matrix result (m);
1310 int a_len = a.length ();
1311 for (int i = 0; i < a_len; i++)
1312 result.elem (i, i) += a.elem (i, i);
1313
1314 return result;
1315 }
1316
1317 Matrix
1318 operator - (const Matrix& m, const DiagMatrix& a)
1319 {
1320 int nr = m.rows ();
1321 int nc = m.cols ();
1322 if (nr != a.rows () || nc != a.cols ())
1323 {
1324 (*current_liboctave_error_handler)
1325 ("nonconformant matrix subtraction attempted");
1326 return Matrix ();
1327 }
1328
1329 if (nr == 0 || nc == 0)
1330 return Matrix (nr, nc);
1331
1332 Matrix result (m);
1333 int a_len = a.length ();
1334 for (int i = 0; i < a_len; i++)
1335 result.elem (i, i) -= a.elem (i, i);
1336
1337 return result;
1338 }
1339
1340 Matrix
1341 operator * (const Matrix& m, const DiagMatrix& a)
1342 {
1343 int nr = m.rows ();
1344 int nc = m.cols ();
1345 int a_nr = a.rows ();
1346 int a_nc = a.cols ();
1347 if (nc != a_nr)
1348 {
1349 (*current_liboctave_error_handler)
1350 ("nonconformant matrix multiplication attempted");
1351 return Matrix ();
1352 }
1353
1354 if (nr == 0 || nc == 0 || a_nc == 0)
1355 return Matrix (nr, a_nc, 0.0);
1356
1357 double *c = new double [nr*a_nc];
1358 double *ctmp = 0;
1359
1360 int a_len = a.length ();
1361 for (int j = 0; j < a_len; j++)
1362 {
1363 int idx = j * nr;
1364 ctmp = c + idx;
1365 if (a.elem (j, j) == 1.0)
1366 {
1367 for (int i = 0; i < nr; i++)
1368 ctmp[i] = m.elem (i, j);
1369 }
1370 else if (a.elem (j, j) == 0.0)
1371 {
1372 for (int i = 0; i < nr; i++)
1373 ctmp[i] = 0.0;
1374 }
1375 else
1376 {
1377 for (int i = 0; i < nr; i++)
1378 ctmp[i] = a.elem (j, j) * m.elem (i, j);
1379 }
1380 }
1381
1382 if (a_nr < a_nc)
1383 {
1384 for (int i = nr * nc; i < nr * a_nc; i++)
1385 ctmp[i] = 0.0;
1386 }
1387
1388 return Matrix (c, nr, a_nc);
1389 }
1390
1391 // diagonal matrix by matrix -> matrix operations
1392
1393 Matrix
1394 operator + (const DiagMatrix& m, const Matrix& a)
1395 {
1396 int nr = m.rows ();
1397 int nc = m.cols ();
1398 if (nr != a.rows () || nc != a.cols ())
1399 {
1400 (*current_liboctave_error_handler)
1401 ("nonconformant matrix addition attempted");
1402 return Matrix ();
1403 }
1404
1405 if (nr == 0 || nc == 0)
1406 return Matrix (nr, nc);
1407
1408 Matrix result (a);
1409 for (int i = 0; i < m.length (); i++)
1410 result.elem (i, i) += m.elem (i, i);
1411
1412 return result;
1413 }
1414
1415 Matrix
1416 operator - (const DiagMatrix& m, const Matrix& a)
1417 {
1418 int nr = m.rows ();
1419 int nc = m.cols ();
1420 if (nr != a.rows () || nc != a.cols ())
1421 {
1422 (*current_liboctave_error_handler)
1423 ("nonconformant matrix subtraction attempted");
1424 return Matrix ();
1425 }
1426
1427 if (nr == 0 || nc == 0)
1428 return Matrix (nr, nc);
1429
1430 Matrix result (-a);
1431 for (int i = 0; i < m.length (); i++)
1432 result.elem (i, i) += m.elem (i, i);
1433
1434 return result;
1435 }
1436
1437 Matrix
1438 operator * (const DiagMatrix& m, const Matrix& a)
1439 {
1440 int nr = m.rows ();
1441 int nc = m.cols ();
1442 int a_nr = a.rows ();
1443 int a_nc = a.cols ();
1444 if (nc != a_nr)
1445 {
1446 (*current_liboctave_error_handler)
1447 ("nonconformant matrix multiplication attempted");
1448 return Matrix ();
1449 }
1450
1451 if (nr == 0 || nc == 0 || a_nc == 0)
1452 return Matrix (nr, a_nc, 0.0);
1453
1454 Matrix c (nr, a_nc);
1455
1456 for (int i = 0; i < m.length (); i++)
1457 {
1458 if (m.elem (i, i) == 1.0)
1459 {
1460 for (int j = 0; j < a_nc; j++)
1461 c.elem (i, j) = a.elem (i, j);
1462 }
1463 else if (m.elem (i, i) == 0.0)
1464 {
1465 for (int j = 0; j < a_nc; j++)
1466 c.elem (i, j) = 0.0;
1467 }
1468 else
1469 {
1470 for (int j = 0; j < a_nc; j++)
1471 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
1472 }
1473 }
1474
1475 if (nr > nc)
1476 {
1477 for (int j = 0; j < a_nc; j++)
1478 for (int i = a_nr; i < nr; i++)
1479 c.elem (i, j) = 0.0;
1480 }
1481
1482 return c;
1483 }
1484
1485 // matrix by matrix -> matrix operations
1486
1487 Matrix
1488 operator * (const Matrix& m, const Matrix& a)
1489 {
1490 int nr = m.rows ();
1491 int nc = m.cols ();
1492 int a_nr = a.rows ();
1493 int a_nc = a.cols ();
1494 if (nc != a_nr)
1495 {
1496 (*current_liboctave_error_handler)
1497 ("nonconformant matrix multiplication attempted");
1498 return Matrix ();
1499 }
1500
1501 if (nr == 0 || nc == 0 || a_nc == 0)
1502 return Matrix (nr, a_nc, 0.0);
1503
1504 int ld = nr;
1505 int lda = a_nr;
1506
1507 double *c = new double [nr*a_nc];
1508
1509 F77_FCN (dgemm, DGEMM) ("N", "N", nr, a_nc, nc, 1.0, m.data (),
1510 ld, a.data (), lda, 0.0, c, nr, 1L, 1L);
1511
1512 return Matrix (c, nr, a_nc);
1513 }
1514
1515 // other operations.
1516
1517 Matrix
1518 map (d_d_Mapper f, const Matrix& a)
1519 {
1520 Matrix b (a);
1521 b.map (f);
1522 return b;
1523 }
1524
1525 Matrix
1526 map (d_c_Mapper f, const ComplexMatrix& a)
1527 {
1528 int a_nc = a.cols ();
1529 int a_nr = a.rows ();
1530 Matrix b (a_nr, a_nc);
1531 for (int j = 0; j < a_nc; j++)
1532 for (int i = 0; i < a_nr; i++)
1533 b.elem (i, j) = f (a.elem (i, j));
1534 return b;
1535 }
1536
1537 void
1538 Matrix::map (d_d_Mapper f)
1539 {
1540 double *d = fortran_vec (); // Ensures only one reference to my privates!
1541
1542 for (int i = 0; i < length (); i++)
1543 d[i] = f (d[i]);
1544 }
1545
1546 // XXX FIXME XXX Do these really belong here? They should maybe be
1547 // cleaned up a bit, no? What about corresponding functions for the
1548 // Vectors?
1549
1550 Matrix
1551 Matrix::all (void) const
1552 {
1553 int nr = rows ();
1554 int nc = cols ();
1555 Matrix retval;
1556 if (nr > 0 && nc > 0)
1557 {
1558 if (nr == 1)
1559 {
1560 retval.resize (1, 1);
1561 retval.elem (0, 0) = 1.0;
1562 for (int j = 0; j < nc; j++)
1563 {
1564 if (elem (0, j) == 0.0)
1565 {
1566 retval.elem (0, 0) = 0.0;
1567 break;
1568 }
1569 }
1570 }
1571 else if (nc == 1)
1572 {
1573 retval.resize (1, 1);
1574 retval.elem (0, 0) = 1.0;
1575 for (int i = 0; i < nr; i++)
1576 {
1577 if (elem (i, 0) == 0.0)
1578 {
1579 retval.elem (0, 0) = 0.0;
1580 break;
1581 }
1582 }
1583 }
1584 else
1585 {
1586 retval.resize (1, nc);
1587 for (int j = 0; j < nc; j++)
1588 {
1589 retval.elem (0, j) = 1.0;
1590 for (int i = 0; i < nr; i++)
1591 {
1592 if (elem (i, j) == 0.0)
1593 {
1594 retval.elem (0, j) = 0.0;
1595 break;
1596 }
1597 }
1598 }
1599 }
1600 }
1601 return retval;
1602 }
1603
1604 Matrix
1605 Matrix::any (void) const
1606 {
1607 int nr = rows ();
1608 int nc = cols ();
1609 Matrix retval;
1610 if (nr > 0 && nc > 0)
1611 {
1612 if (nr == 1)
1613 {
1614 retval.resize (1, 1);
1615 retval.elem (0, 0) = 0.0;
1616 for (int j = 0; j < nc; j++)
1617 {
1618 if (elem (0, j) != 0.0)
1619 {
1620 retval.elem (0, 0) = 1.0;
1621 break;
1622 }
1623 }
1624 }
1625 else if (nc == 1)
1626 {
1627 retval.resize (1, 1);
1628 retval.elem (0, 0) = 0.0;
1629 for (int i = 0; i < nr; i++)
1630 {
1631 if (elem (i, 0) != 0.0)
1632 {
1633 retval.elem (0, 0) = 1.0;
1634 break;
1635 }
1636 }
1637 }
1638 else
1639 {
1640 retval.resize (1, nc);
1641 for (int j = 0; j < nc; j++)
1642 {
1643 retval.elem (0, j) = 0.0;
1644 for (int i = 0; i < nr; i++)
1645 {
1646 if (elem (i, j) != 0.0)
1647 {
1648 retval.elem (0, j) = 1.0;
1649 break;
1650 }
1651 }
1652 }
1653 }
1654 }
1655 return retval;
1656 }
1657
1658 Matrix
1659 Matrix::cumprod (void) const
1660 {
1661 Matrix retval;
1662
1663 int nr = rows ();
1664 int nc = cols ();
1665
1666 if (nr == 1)
1667 {
1668 retval.resize (1, nc);
1669 if (nc > 0)
1670 {
1671 double prod = elem (0, 0);
1672 for (int j = 0; j < nc; j++)
1673 {
1674 retval.elem (0, j) = prod;
1675 if (j < nc - 1)
1676 prod *= elem (0, j+1);
1677 }
1678 }
1679 }
1680 else if (nc == 1)
1681 {
1682 retval.resize (nr, 1);
1683 if (nr > 0)
1684 {
1685 double prod = elem (0, 0);
1686 for (int i = 0; i < nr; i++)
1687 {
1688 retval.elem (i, 0) = prod;
1689 if (i < nr - 1)
1690 prod *= elem (i+1, 0);
1691 }
1692 }
1693 }
1694 else
1695 {
1696 retval.resize (nr, nc);
1697 if (nr > 0 && nc > 0)
1698 {
1699 for (int j = 0; j < nc; j++)
1700 {
1701 double prod = elem (0, j);
1702 for (int i = 0; i < nr; i++)
1703 {
1704 retval.elem (i, j) = prod;
1705 if (i < nr - 1)
1706 prod *= elem (i+1, j);
1707 }
1708 }
1709 }
1710 }
1711 return retval;
1712 }
1713
1714 Matrix
1715 Matrix::cumsum (void) const
1716 {
1717 Matrix retval;
1718
1719 int nr = rows ();
1720 int nc = cols ();
1721
1722 if (nr == 1)
1723 {
1724 retval.resize (1, nc);
1725 if (nc > 0)
1726 {
1727 double sum = elem (0, 0);
1728 for (int j = 0; j < nc; j++)
1729 {
1730 retval.elem (0, j) = sum;
1731 if (j < nc - 1)
1732 sum += elem (0, j+1);
1733 }
1734 }
1735 }
1736 else if (nc == 1)
1737 {
1738 retval.resize (nr, 1);
1739 if (nr > 0)
1740 {
1741 double sum = elem (0, 0);
1742 for (int i = 0; i < nr; i++)
1743 {
1744 retval.elem (i, 0) = sum;
1745 if (i < nr - 1)
1746 sum += elem (i+1, 0);
1747 }
1748 }
1749 }
1750 else
1751 {
1752 retval.resize (nr, nc);
1753 if (nr > 0 && nc > 0)
1754 {
1755 for (int j = 0; j < nc; j++)
1756 {
1757 double sum = elem (0, j);
1758 for (int i = 0; i < nr; i++)
1759 {
1760 retval.elem (i, j) = sum;
1761 if (i < nr - 1)
1762 sum += elem (i+1, j);
1763 }
1764 }
1765 }
1766 }
1767 return retval;
1768 }
1769
1770 Matrix
1771 Matrix::prod (void) const
1772 {
1773 Matrix retval;
1774
1775 int nr = rows ();
1776 int nc = cols ();
1777
1778 if (nr == 1)
1779 {
1780 retval.resize (1, 1);
1781 retval.elem (0, 0) = 1.0;
1782 for (int j = 0; j < nc; j++)
1783 retval.elem (0, 0) *= elem (0, j);
1784 }
1785 else if (nc == 1)
1786 {
1787 retval.resize (1, 1);
1788 retval.elem (0, 0) = 1.0;
1789 for (int i = 0; i < nr; i++)
1790 retval.elem (0, 0) *= elem (i, 0);
1791 }
1792 else
1793 {
1794 if (nc == 0)
1795 {
1796 retval.resize (1, 1);
1797 retval.elem (0, 0) = 1.0;
1798 }
1799 else
1800 retval.resize (1, nc);
1801
1802 for (int j = 0; j < nc; j++)
1803 {
1804 retval.elem (0, j) = 1.0;
1805 for (int i = 0; i < nr; i++)
1806 retval.elem (0, j) *= elem (i, j);
1807 }
1808 }
1809 return retval;
1810 }
1811
1812 Matrix
1813 Matrix::sum (void) const
1814 {
1815 Matrix retval;
1816
1817 int nr = rows ();
1818 int nc = cols ();
1819
1820 if (nr == 1)
1821 {
1822 retval.resize (1, 1);
1823 retval.elem (0, 0) = 0.0;
1824 for (int j = 0; j < nc; j++)
1825 retval.elem (0, 0) += elem (0, j);
1826 }
1827 else if (nc == 1)
1828 {
1829 retval.resize (1, 1);
1830 retval.elem (0, 0) = 0.0;
1831 for (int i = 0; i < nr; i++)
1832 retval.elem (0, 0) += elem (i, 0);
1833 }
1834 else
1835 {
1836 if (nc == 0)
1837 {
1838 retval.resize (1, 1);
1839 retval.elem (0, 0) = 0.0;
1840 }
1841 else
1842 retval.resize (1, nc);
1843
1844 for (int j = 0; j < nc; j++)
1845 {
1846 retval.elem (0, j) = 0.0;
1847 for (int i = 0; i < nr; i++)
1848 retval.elem (0, j) += elem (i, j);
1849 }
1850 }
1851 return retval;
1852 }
1853
1854 Matrix
1855 Matrix::sumsq (void) const
1856 {
1857 Matrix retval;
1858
1859 int nr = rows ();
1860 int nc = cols ();
1861
1862 if (nr == 1)
1863 {
1864 retval.resize (1, 1);
1865 retval.elem (0, 0) = 0.0;
1866 for (int j = 0; j < nc; j++)
1867 {
1868 double d = elem (0, j);
1869 retval.elem (0, 0) += d * d;
1870 }
1871 }
1872 else if (nc == 1)
1873 {
1874 retval.resize (1, 1);
1875 retval.elem (0, 0) = 0.0;
1876 for (int i = 0; i < nr; i++)
1877 {
1878 double d = elem (i, 0);
1879 retval.elem (0, 0) += d * d;
1880 }
1881 }
1882 else
1883 {
1884 retval.resize (1, nc);
1885 for (int j = 0; j < nc; j++)
1886 {
1887 retval.elem (0, j) = 0.0;
1888 for (int i = 0; i < nr; i++)
1889 {
1890 double d = elem (i, j);
1891 retval.elem (0, j) += d * d;
1892 }
1893 }
1894 }
1895 return retval;
1896 }
1897
1898 ColumnVector
1899 Matrix::diag (void) const
1900 {
1901 return diag (0);
1902 }
1903
1904 ColumnVector
1905 Matrix::diag (int k) const
1906 {
1907 int nnr = rows ();
1908 int nnc = cols ();
1909 if (k > 0)
1910 nnc -= k;
1911 else if (k < 0)
1912 nnr += k;
1913
1914 ColumnVector d;
1915
1916 if (nnr > 0 && nnc > 0)
1917 {
1918 int ndiag = (nnr < nnc) ? nnr : nnc;
1919
1920 d.resize (ndiag);
1921
1922 if (k > 0)
1923 {
1924 for (int i = 0; i < ndiag; i++)
1925 d.elem (i) = elem (i, i+k);
1926 }
1927 else if ( k < 0)
1928 {
1929 for (int i = 0; i < ndiag; i++)
1930 d.elem (i) = elem (i-k, i);
1931 }
1932 else
1933 {
1934 for (int i = 0; i < ndiag; i++)
1935 d.elem (i) = elem (i, i);
1936 }
1937 }
1938 else
1939 cerr << "diag: requested diagonal out of range\n";
1940
1941 return d;
1942 }
1943
1944 ColumnVector
1945 Matrix::row_min (void) const
1946 {
1947 ColumnVector result;
1948
1949 int nr = rows ();
1950 int nc = cols ();
1951
1952 if (nr > 0 && nc > 0)
1953 {
1954 result.resize (nr);
1955
1956 for (int i = 0; i < nr; i++)
1957 {
1958 double res = elem (i, 0);
1959 for (int j = 1; j < nc; j++)
1960 if (elem (i, j) < res)
1961 res = elem (i, j);
1962 result.elem (i) = res;
1963 }
1964 }
1965
1966 return result;
1967 }
1968
1969 ColumnVector
1970 Matrix::row_min_loc (void) const
1971 {
1972 ColumnVector result;
1973
1974 int nr = rows ();
1975 int nc = cols ();
1976
1977 if (nr > 0 && nc > 0)
1978 {
1979 result.resize (nr);
1980
1981 for (int i = 0; i < nr; i++)
1982 {
1983 int res = 0;
1984 for (int j = 0; j < nc; j++)
1985 if (elem (i, j) < elem (i, res))
1986 res = j;
1987 result.elem (i) = (double) (res + 1);
1988 }
1989 }
1990
1991 return result;
1992 }
1993
1994 ColumnVector
1995 Matrix::row_max (void) const
1996 {
1997 ColumnVector result;
1998
1999 int nr = rows ();
2000 int nc = cols ();
2001
2002 if (nr > 0 && nc > 0)
2003 {
2004 result.resize (nr);
2005
2006 for (int i = 0; i < nr; i++)
2007 {
2008 double res = elem (i, 0);
2009 for (int j = 1; j < nc; j++)
2010 if (elem (i, j) > res)
2011 res = elem (i, j);
2012 result.elem (i) = res;
2013 }
2014 }
2015
2016 return result;
2017 }
2018
2019 ColumnVector
2020 Matrix::row_max_loc (void) const
2021 {
2022 ColumnVector result;
2023
2024 int nr = rows ();
2025 int nc = cols ();
2026
2027 if (nr > 0 && nc > 0)
2028 {
2029 result.resize (nr);
2030
2031 for (int i = 0; i < nr; i++)
2032 {
2033 int res = 0;
2034 for (int j = 0; j < nc; j++)
2035 if (elem (i, j) > elem (i, res))
2036 res = j;
2037 result.elem (i) = (double) (res + 1);
2038 }
2039 }
2040
2041 return result;
2042 }
2043
2044 RowVector
2045 Matrix::column_min (void) const
2046 {
2047 RowVector result;
2048
2049 int nr = rows ();
2050 int nc = cols ();
2051
2052 if (nr > 0 && nc > 0)
2053 {
2054 result.resize (nc);
2055
2056 for (int j = 0; j < nc; j++)
2057 {
2058 double res = elem (0, j);
2059 for (int i = 1; i < nr; i++)
2060 if (elem (i, j) < res)
2061 res = elem (i, j);
2062 result.elem (j) = res;
2063 }
2064 }
2065
2066 return result;
2067 }
2068 RowVector
2069 Matrix::column_min_loc (void) const
2070 {
2071 RowVector result;
2072
2073 int nr = rows ();
2074 int nc = cols ();
2075
2076 if (nr > 0 && nc > 0)
2077 {
2078 result.resize (nc);
2079
2080 for (int j = 0; j < nc; j++)
2081 {
2082 int res = 0;
2083 for (int i = 0; i < nr; i++)
2084 if (elem (i, j) < elem (res, j))
2085 res = i;
2086 result.elem (j) = (double) (res + 1);
2087 }
2088 }
2089
2090 return result;
2091 }
2092
2093
2094 RowVector
2095 Matrix::column_max (void) const
2096 {
2097 RowVector result;
2098
2099 int nr = rows ();
2100 int nc = cols ();
2101
2102 if (nr > 0 && nc > 0)
2103 {
2104 result.resize (nc);
2105
2106 for (int j = 0; j < nc; j++)
2107 {
2108 double res = elem (0, j);
2109 for (int i = 1; i < nr; i++)
2110 if (elem (i, j) > res)
2111 res = elem (i, j);
2112 result.elem (j) = res;
2113 }
2114 }
2115
2116 return result;
2117 }
2118
2119 RowVector
2120 Matrix::column_max_loc (void) const
2121 {
2122 RowVector result;
2123
2124 int nr = rows ();
2125 int nc = cols ();
2126
2127 if (nr > 0 && nc > 0)
2128 {
2129 result.resize (nc);
2130
2131 for (int j = 0; j < nc; j++)
2132 {
2133 int res = 0;
2134 for (int i = 0; i < nr; i++)
2135 if (elem (i, j) > elem (res, j))
2136 res = i;
2137 result.elem (j) = (double) (res + 1);
2138 }
2139 }
2140
2141 return result;
2142 }
2143
2144 ostream&
2145 operator << (ostream& os, const Matrix& a)
2146 {
2147 // int field_width = os.precision () + 7;
2148
2149 for (int i = 0; i < a.rows (); i++)
2150 {
2151 for (int j = 0; j < a.cols (); j++)
2152 os << " " /* setw (field_width) */ << a.elem (i, j);
2153 os << "\n";
2154 }
2155 return os;
2156 }
2157
2158 istream&
2159 operator >> (istream& is, Matrix& a)
2160 {
2161 int nr = a.rows ();
2162 int nc = a.cols ();
2163
2164 if (nr < 1 || nc < 1)
2165 is.clear (ios::badbit);
2166 else
2167 {
2168 double tmp;
2169 for (int i = 0; i < nr; i++)
2170 for (int j = 0; j < nc; j++)
2171 {
2172 is >> tmp;
2173 if (is)
2174 a.elem (i, j) = tmp;
2175 else
2176 break;
2177 }
2178 }
2179
2180 return is;
2181 }
2182
2183 // Read an array of data from a file in binary format.
2184
2185 int
2186 Matrix::read (FILE *fptr, const char *type)
2187 {
2188 // Allocate buffer pointers.
2189
2190 union
2191 {
2192 void *vd;
2193 char *ch;
2194 u_char *uc;
2195 short *sh;
2196 u_short *us;
2197 int *in;
2198 u_int *ui;
2199 long *ln;
2200 u_long *ul;
2201 float *fl;
2202 double *db;
2203 }
2204 buf;
2205
2206 // Convert data to double.
2207
2208 if (! type)
2209 {
2210 (*current_liboctave_error_handler)
2211 ("fread: invalid NULL type parameter");
2212 return 0;
2213 }
2214
2215 int count;
2216 int nitems = length ();
2217
2218 double *d = fortran_vec (); // Ensures only one reference to my privates!
2219
2220 #define DO_FREAD(TYPE,ELEM) \
2221 do \
2222 { \
2223 size_t size = sizeof (TYPE); \
2224 buf.ch = new char [size * nitems]; \
2225 count = fread (buf.ch, size, nitems, fptr); \
2226 for (int k = 0; k < count; k++) \
2227 d[k] = buf.ELEM[k]; \
2228 delete [] buf.ch; \
2229 } \
2230 while (0)
2231
2232 if (strcasecmp (type, "double") == 0)
2233 DO_FREAD (double, db);
2234 else if (strcasecmp (type, "char") == 0)
2235 DO_FREAD (char, ch);
2236 else if (strcasecmp (type, "uchar") == 0)
2237 DO_FREAD (u_char, uc);
2238 else if (strcasecmp (type, "short") == 0)
2239 DO_FREAD (short, sh);
2240 else if (strcasecmp (type, "ushort") == 0)
2241 DO_FREAD (u_short, us);
2242 else if (strcasecmp (type, "int") == 0)
2243 DO_FREAD (int, in);
2244 else if (strcasecmp (type, "uint") == 0)
2245 DO_FREAD (u_int, ui);
2246 else if (strcasecmp (type, "long") == 0)
2247 DO_FREAD (long, ul);
2248 else if (strcasecmp (type, "float") == 0)
2249 DO_FREAD (float, fl);
2250 else
2251 {
2252 (*current_liboctave_error_handler)
2253 ("fread: invalid NULL type parameter");
2254 return 0;
2255 }
2256
2257 return count;
2258 }
2259
2260 // Write the data array to a file in binary format.
2261
2262 int
2263 Matrix::write (FILE *fptr, const char *type)
2264 {
2265 // Allocate buffer pointers.
2266
2267 union
2268 {
2269 void *vd;
2270 char *ch;
2271 u_char *uc;
2272 short *sh;
2273 u_short *us;
2274 int *in;
2275 u_int *ui;
2276 long *ln;
2277 u_long *ul;
2278 float *fl;
2279 double *db;
2280 }
2281 buf;
2282
2283 int nitems = length ();
2284
2285 double *d = fortran_vec ();
2286
2287 // Convert from double to correct size.
2288
2289 if (! type)
2290 {
2291 (*current_liboctave_error_handler)
2292 ("fwrite: invalid NULL type parameter");
2293 return 0;
2294 }
2295
2296 size_t size;
2297 int count;
2298
2299 #define DO_FWRITE(TYPE,ELEM) \
2300 do \
2301 { \
2302 size = sizeof (TYPE); \
2303 buf.ELEM = new TYPE [nitems]; \
2304 for (int k = 0; k < nitems; k++) \
2305 buf.ELEM[k] = (TYPE) d[k]; \
2306 count = fwrite (buf.ELEM, size, nitems, fptr); \
2307 delete [] buf.ELEM; \
2308 } \
2309 while (0)
2310
2311 if (strcasecmp (type, "double") == 0)
2312 DO_FWRITE (double, db);
2313 else if (strcasecmp (type, "char") == 0)
2314 DO_FWRITE (char, ch);
2315 else if (strcasecmp (type, "uchar") == 0)
2316 DO_FWRITE (u_char, uc);
2317 else if (strcasecmp (type, "short") == 0)
2318 DO_FWRITE (short, sh);
2319 else if (strcasecmp (type, "ushort") == 0)
2320 DO_FWRITE (u_short, us);
2321 else if (strcasecmp (type, "int") == 0)
2322 DO_FWRITE (int, in);
2323 else if (strcasecmp (type, "uint") == 0)
2324 DO_FWRITE (u_int, ui);
2325 else if (strcasecmp (type, "long") == 0)
2326 DO_FWRITE (long, ln);
2327 else if (strcasecmp (type, "ulong") == 0)
2328 DO_FWRITE (u_long, ul);
2329 else if (strcasecmp (type, "float") == 0)
2330 DO_FWRITE (float, fl);
2331 else
2332 {
2333 (*current_liboctave_error_handler)
2334 ("fwrite: unrecognized type parameter %s", type);
2335 return 0;
2336 }
2337
2338 return count;
2339 }
2340 #endif
2341
2342 /*
2343 ;;; Local Variables: ***
2344 ;;; mode: C++ ***
2345 ;;; page-delimiter: "^/\\*" ***
2346 ;;; End: ***
2347 */