Mercurial > hg > octave-nkf

comparison liboctave/dDiagMatrix.cc @ 1205:8302fab9fe24

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-04-04 02:05:01 by jwe]
author jwe
date Tue, 04 Apr 1995 02:05:01 +0000
parents b6360f2d4fa6
children 0bf4d2b7def4
comparison
equal deleted inserted replaced
1204:68d147abe7ca 1205:8302fab9fe24
43 #include "mx-kludge.cc" 43 #include "mx-kludge.cc"
44 #undef KLUDGE_DIAG_MATRICES 44 #undef KLUDGE_DIAG_MATRICES
45 #undef TYPE 45 #undef TYPE
46 #undef KL_DMAT_TYPE 46 #undef KL_DMAT_TYPE
47 47
48 #if 0
49 DiagMatrix&
50 DiagMatrix::resize (int r, int c)
51 {
52 if (r < 0 || c < 0)
53 {
54 (*current_liboctave_error_handler)
55 ("can't resize to negative dimensions");
56 return *this;
57 }
58
59 int new_len = r < c ? r : c;
60 double *new_data = 0;
61 if (new_len > 0)
62 {
63 new_data = new double [new_len];
64
65 int min_len = new_len < len ? new_len : len;
66
67 for (int i = 0; i < min_len; i++)
68 new_data[i] = data[i];
69 }
70
71 delete [] data;
72 nr = r;
73 nc = c;
74 len = new_len;
75 data = new_data;
76
77 return *this;
78 }
79
80 DiagMatrix&
81 DiagMatrix::resize (int r, int c, double val)
82 {
83 if (r < 0 || c < 0)
84 {
85 (*current_liboctave_error_handler)
86 ("can't resize to negative dimensions");
87 return *this;
88 }
89
90 int new_len = r < c ? r : c;
91 double *new_data = 0;
92 if (new_len > 0)
93 {
94 new_data = new double [new_len];
95
96 int min_len = new_len < len ? new_len : len;
97
98 for (int i = 0; i < min_len; i++)
99 new_data[i] = data[i];
100
101 for (i = min_len; i < new_len; i++)
102 new_data[i] = val;
103 }
104
105 delete [] data;
106 nr = r;
107 nc = c;
108 len = new_len;
109 data = new_data;
110
111 return *this;
112 }
113 #endif
114
115 int 48 int
116 DiagMatrix::operator == (const DiagMatrix& a) const 49 DiagMatrix::operator == (const DiagMatrix& a) const
117 { 50 {
118 if (rows () != a.rows () || cols () != a.cols ()) 51 if (rows () != a.rows () || cols () != a.cols ())
119 return 0; 52 return 0;
216 149
217 DiagMatrix 150 DiagMatrix
218 DiagMatrix::transpose (void) const 151 DiagMatrix::transpose (void) const
219 { 152 {
220 return DiagMatrix (dup (data (), length ()), cols (), rows ()); 153 return DiagMatrix (dup (data (), length ()), cols (), rows ());
154 }
155
156 DiagMatrix
157 real (const ComplexDiagMatrix& a)
158 {
159 DiagMatrix retval;
160 int a_len = a.length ();
161 if (a_len > 0)
162 retval = DiagMatrix (real_dup (a.data (), a_len), a.rows (),
163 a.cols ());
164 return retval;
165 }
166
167 DiagMatrix
168 imag (const ComplexDiagMatrix& a)
169 {
170 DiagMatrix retval;
171 int a_len = a.length ();
172 if (a_len > 0)
173 retval = DiagMatrix (imag_dup (a.data (), a_len), a.rows (),
174 a.cols ());
175 return retval;
221 } 176 }
222 177
223 Matrix 178 Matrix
224 DiagMatrix::extract (int r1, int c1, int r2, int c2) const 179 DiagMatrix::extract (int r1, int c1, int r2, int c2) const
225 { 180 {
396 351
397 double *d = fortran_vec (); // Ensures only one reference to my privates! 352 double *d = fortran_vec (); // Ensures only one reference to my privates!
398 353
399 subtract2 (d, a.data (), length ()); 354 subtract2 (d, a.data (), length ());
400 return *this; 355 return *this;
401 }
402
403 // diagonal matrix by scalar -> matrix operations
404
405 Matrix
406 operator + (const DiagMatrix& a, double s)
407 {
408 Matrix tmp (a.rows (), a.cols (), s);
409 return a + tmp;
410 }
411
412 Matrix
413 operator - (const DiagMatrix& a, double s)
414 {
415 Matrix tmp (a.rows (), a.cols (), -s);
416 return a + tmp;
417 }
418
419 ComplexMatrix
420 operator + (const DiagMatrix& a, const Complex& s)
421 {
422 ComplexMatrix tmp (a.rows (), a.cols (), s);
423 return a + tmp;
424 }
425
426 ComplexMatrix
427 operator - (const DiagMatrix& a, const Complex& s)
428 {
429 ComplexMatrix tmp (a.rows (), a.cols (), -s);
430 return a + tmp;
431 }
432
433 // diagonal matrix by scalar -> diagonal matrix operations
434
435 ComplexDiagMatrix
436 operator * (const DiagMatrix& a, const Complex& s)
437 {
438 return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
439 a.rows (), a.cols ());
440 }
441
442 ComplexDiagMatrix
443 operator / (const DiagMatrix& a, const Complex& s)
444 {
445 return ComplexDiagMatrix (divide (a.data (), a.length (), s),
446 a.rows (), a.cols ());
447 }
448
449 // scalar by diagonal matrix -> matrix operations
450
451 Matrix
452 operator + (double s, const DiagMatrix& a)
453 {
454 Matrix tmp (a.rows (), a.cols (), s);
455 return tmp + a;
456 }
457
458 Matrix
459 operator - (double s, const DiagMatrix& a)
460 {
461 Matrix tmp (a.rows (), a.cols (), s);
462 return tmp - a;
463 }
464
465 ComplexMatrix
466 operator + (const Complex& s, const DiagMatrix& a)
467 {
468 ComplexMatrix tmp (a.rows (), a.cols (), s);
469 return tmp + a;
470 }
471
472 ComplexMatrix
473 operator - (const Complex& s, const DiagMatrix& a)
474 {
475 ComplexMatrix tmp (a.rows (), a.cols (), s);
476 return tmp - a;
477 }
478
479 // scalar by diagonal matrix -> diagonal matrix operations
480
481 ComplexDiagMatrix
482 operator * (const Complex& s, const DiagMatrix& a)
483 {
484 return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
485 a.rows (), a.cols ());
486 }
487
488 // diagonal matrix by column vector -> column vector operations
489
490 ColumnVector
491 operator * (const DiagMatrix& m, const ColumnVector& a)
492 {
493 int nr = m.rows ();
494 int nc = m.cols ();
495 int a_len = a.length ();
496 if (nc != a_len)
497 {
498 (*current_liboctave_error_handler)
499 ("nonconformant matrix multiplication attempted");
500 return ColumnVector ();
501 }
502
503 if (nc == 0 || nr == 0)
504 return ColumnVector (0);
505
506 ColumnVector result (nr);
507
508 for (int i = 0; i < a_len; i++)
509 result.elem (i) = a.elem (i) * m.elem (i, i);
510
511 for (i = a_len; i < nr; i++)
512 result.elem (i) = 0.0;
513
514 return result;
515 }
516
517 ComplexColumnVector
518 operator * (const DiagMatrix& m, const ComplexColumnVector& a)
519 {
520 int nr = m.rows ();
521 int nc = m.cols ();
522 int a_len = a.length ();
523 if (nc != a_len)
524 {
525 (*current_liboctave_error_handler)
526 ("nonconformant matrix multiplication attempted");
527 return ColumnVector ();
528 }
529
530 if (nc == 0 || nr == 0)
531 return ComplexColumnVector (0);
532
533 ComplexColumnVector result (nr);
534
535 for (int i = 0; i < a_len; i++)
536 result.elem (i) = a.elem (i) * m.elem (i, i);
537
538 for (i = a_len; i < nr; i++)
539 result.elem (i) = 0.0;
540
541 return result;
542 } 356 }
543 357
544 // diagonal matrix by diagonal matrix -> diagonal matrix operations 358 // diagonal matrix by diagonal matrix -> diagonal matrix operations
545 359
546 DiagMatrix 360 DiagMatrix
580 } 394 }
581 395
582 return c; 396 return c;
583 } 397 }
584 398
585 ComplexDiagMatrix
586 operator + (const DiagMatrix& m, const ComplexDiagMatrix& a)
587 {
588 int nr = m.rows ();
589 int nc = m.cols ();
590 if (nr != a.rows () || nc != a.cols ())
591 {
592 (*current_liboctave_error_handler)
593 ("nonconformant matrix addition attempted");
594 return ComplexDiagMatrix ();
595 }
596
597 if (nc == 0 || nr == 0)
598 return ComplexDiagMatrix (nr, nc);
599
600 return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
601 }
602
603 ComplexDiagMatrix
604 operator - (const DiagMatrix& m, const ComplexDiagMatrix& a)
605 {
606 int nr = m.rows ();
607 int nc = m.cols ();
608 if (nr != a.rows () || nc != a.cols ())
609 {
610 (*current_liboctave_error_handler)
611 ("nonconformant matrix subtraction attempted");
612 return ComplexDiagMatrix ();
613 }
614
615 if (nc == 0 || nr == 0)
616 return ComplexDiagMatrix (nr, nc);
617
618 return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
619 nr, nc);
620 }
621
622 ComplexDiagMatrix
623 operator * (const DiagMatrix& a, const ComplexDiagMatrix& b)
624 {
625 int nr_a = a.rows ();
626 int nc_a = a.cols ();
627 int nr_b = b.rows ();
628 int nc_b = b.cols ();
629 if (nc_a != nr_b)
630 {
631 (*current_liboctave_error_handler)
632 ("nonconformant matrix multiplication attempted");
633 return ComplexDiagMatrix ();
634 }
635
636 if (nr_a == 0 || nc_a == 0 || nc_b == 0)
637 return ComplexDiagMatrix (nr_a, nc_a, 0.0);
638
639 ComplexDiagMatrix c (nr_a, nc_b);
640
641 int len = nr_a < nc_b ? nr_a : nc_b;
642
643 for (int i = 0; i < len; i++)
644 {
645 double a_element = a.elem (i, i);
646 Complex b_element = b.elem (i, i);
647
648 if (a_element == 0.0 || b_element == 0.0)
649 c.elem (i, i) = 0.0;
650 else if (a_element == 1.0)
651 c.elem (i, i) = b_element;
652 else if (b_element == 1.0)
653 c.elem (i, i) = a_element;
654 else
655 c.elem (i, i) = a_element * b_element;
656 }
657
658 return c;
659 }
660
661 ComplexDiagMatrix
662 product (const DiagMatrix& m, const ComplexDiagMatrix& a)
663 {
664 int nr = m.rows ();
665 int nc = m.cols ();
666 if (nr != a.rows () || nc != a.cols ())
667 {
668 (*current_liboctave_error_handler)
669 ("nonconformant matrix product attempted");
670 return ComplexDiagMatrix ();
671 }
672
673 if (nc == 0 || nr == 0)
674 return ComplexDiagMatrix (nr, nc);
675
676 return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
677 nr, nc);
678 }
679
680 // diagonal matrix by matrix -> matrix operations
681
682 Matrix
683 operator + (const DiagMatrix& m, const Matrix& a)
684 {
685 int nr = m.rows ();
686 int nc = m.cols ();
687 if (nr != a.rows () || nc != a.cols ())
688 {
689 (*current_liboctave_error_handler)
690 ("nonconformant matrix addition attempted");
691 return Matrix ();
692 }
693
694 if (nr == 0 || nc == 0)
695 return Matrix (nr, nc);
696
697 Matrix result (a);
698 for (int i = 0; i < m.length (); i++)
699 result.elem (i, i) += m.elem (i, i);
700
701 return result;
702 }
703
704 Matrix
705 operator - (const DiagMatrix& m, const Matrix& a)
706 {
707 int nr = m.rows ();
708 int nc = m.cols ();
709 if (nr != a.rows () || nc != a.cols ())
710 {
711 (*current_liboctave_error_handler)
712 ("nonconformant matrix subtraction attempted");
713 return Matrix ();
714 }
715
716 if (nr == 0 || nc == 0)
717 return Matrix (nr, nc);
718
719 Matrix result (-a);
720 for (int i = 0; i < m.length (); i++)
721 result.elem (i, i) += m.elem (i, i);
722
723 return result;
724 }
725
726 Matrix
727 operator * (const DiagMatrix& m, const Matrix& a)
728 {
729 int nr = m.rows ();
730 int nc = m.cols ();
731 int a_nr = a.rows ();
732 int a_nc = a.cols ();
733 if (nc != a_nr)
734 {
735 (*current_liboctave_error_handler)
736 ("nonconformant matrix multiplication attempted");
737 return Matrix ();
738 }
739
740 if (nr == 0 || nc == 0 || a_nc == 0)
741 return Matrix (nr, a_nc, 0.0);
742
743 Matrix c (nr, a_nc);
744
745 for (int i = 0; i < m.length (); i++)
746 {
747 if (m.elem (i, i) == 1.0)
748 {
749 for (int j = 0; j < a_nc; j++)
750 c.elem (i, j) = a.elem (i, j);
751 }
752 else if (m.elem (i, i) == 0.0)
753 {
754 for (int j = 0; j < a_nc; j++)
755 c.elem (i, j) = 0.0;
756 }
757 else
758 {
759 for (int j = 0; j < a_nc; j++)
760 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
761 }
762 }
763
764 if (nr > nc)
765 {
766 for (int j = 0; j < a_nc; j++)
767 for (int i = a_nr; i < nr; i++)
768 c.elem (i, j) = 0.0;
769 }
770
771 return c;
772 }
773
774 ComplexMatrix
775 operator + (const DiagMatrix& m, const ComplexMatrix& a)
776 {
777 int nr = m.rows ();
778 int nc = m.cols ();
779 if (nr != a.rows () || nc != a.cols ())
780 {
781 (*current_liboctave_error_handler)
782 ("nonconformant matrix addition attempted");
783 return ComplexMatrix ();
784 }
785
786 if (nr == 0 || nc == 0)
787 return ComplexMatrix (nr, nc);
788
789 ComplexMatrix result (a);
790 for (int i = 0; i < m.length (); i++)
791 result.elem (i, i) += m.elem (i, i);
792
793 return result;
794 }
795
796 ComplexMatrix
797 operator - (const DiagMatrix& m, const ComplexMatrix& a)
798 {
799 int nr = m.rows ();
800 int nc = m.cols ();
801 if (nr != a.rows () || nc != a.cols ())
802 {
803 (*current_liboctave_error_handler)
804 ("nonconformant matrix subtraction attempted");
805 return ComplexMatrix ();
806 }
807
808 if (nr == 0 || nc == 0)
809 return ComplexMatrix (nr, nc);
810
811 ComplexMatrix result (-a);
812 for (int i = 0; i < m.length (); i++)
813 result.elem (i, i) += m.elem (i, i);
814
815 return result;
816 }
817
818 ComplexMatrix
819 operator * (const DiagMatrix& m, const ComplexMatrix& a)
820 {
821 int nr = m.rows ();
822 int nc = m.cols ();
823 int a_nr = a.rows ();
824 int a_nc = a.cols ();
825 if (nc != a_nr)
826 {
827 (*current_liboctave_error_handler)
828 ("nonconformant matrix multiplication attempted");
829 return ComplexMatrix ();
830 }
831
832 if (nr == 0 || nc == 0 || a_nc == 0)
833 return ComplexMatrix (nr, nc, 0.0);
834
835 ComplexMatrix c (nr, a_nc);
836
837 for (int i = 0; i < m.length (); i++)
838 {
839 if (m.elem (i, i) == 1.0)
840 {
841 for (int j = 0; j < a_nc; j++)
842 c.elem (i, j) = a.elem (i, j);
843 }
844 else if (m.elem (i, i) == 0.0)
845 {
846 for (int j = 0; j < a_nc; j++)
847 c.elem (i, j) = 0.0;
848 }
849 else
850 {
851 for (int j = 0; j < a_nc; j++)
852 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
853 }
854 }
855
856 if (nr > nc)
857 {
858 for (int j = 0; j < a_nc; j++)
859 for (int i = a_nr; i < nr; i++)
860 c.elem (i, j) = 0.0;
861 }
862
863 return c;
864 }
865
866 // other operations 399 // other operations
867 400
868 ColumnVector 401 ColumnVector
869 DiagMatrix::diag (void) const 402 DiagMatrix::diag (void) const
870 { 403 {