Mercurial > hg > octave-nkf

comparison src/xdiv.cc @ 8366:8b1a2555c4e2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement diagonal matrix objects * * *
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 03 Dec 2008 13:32:57 +0100
parents 82be108cc558
children eb63fbe60fab
comparison
equal deleted inserted replaced
8365:65ca196fff28 8366:8b1a2555c4e2
1 /* 1 /*
2 2
3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003, 3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003,
4 2005, 2006, 2007 John W. Eaton 4 2005, 2006, 2007 John W. Eaton
5 Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
5 6
6 This file is part of Octave. 7 This file is part of Octave.
7 8
8 Octave is free software; you can redistribute it and/or modify it 9 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 under the terms of the GNU General Public License as published by the
35 #include "fCMatrix.h" 36 #include "fCMatrix.h"
36 #include "fMatrix.h" 37 #include "fMatrix.h"
37 #include "fCNDArray.h" 38 #include "fCNDArray.h"
38 #include "fNDArray.h" 39 #include "fNDArray.h"
39 #include "oct-cmplx.h" 40 #include "oct-cmplx.h"
41 #include "dDiagMatrix.h"
42 #include "fDiagMatrix.h"
43 #include "CDiagMatrix.h"
44 #include "fCDiagMatrix.h"
40 #include "quit.h" 45 #include "quit.h"
41 46
42 #include "error.h" 47 #include "error.h"
43 #include "xdiv.h" 48 #include "xdiv.h"
44 49
720 octave_idx_type info; 725 octave_idx_type info;
721 float rcond = 0.0; 726 float rcond = 0.0;
722 return a.solve (typ, b, info, rcond, solve_singularity_warning); 727 return a.solve (typ, b, info, rcond, solve_singularity_warning);
723 } 728 }
724 729
730 // Diagonal matrix division.
731
732 template <class MT, class DMT>
733 MT
734 mdm_div_impl (const MT& a, const DMT& d)
735 {
736 if (! mx_div_conform (a, d))
737 return MT ();
738
739 octave_idx_type m = a.rows (), n = d.rows (), k = d.cols ();
740 MT x (m, n);
741 const typename DMT::element_type zero = typename DMT::element_type ();
742
743 for (octave_idx_type j = 0; j < n; j++)
744 {
745 if (j < k && d(j, j) != zero)
746 {
747 for (octave_idx_type i = 0; i < m; i++)
748 x(i, j) = a(i, j) / d(j, j);
749 }
750 else
751 {
752 for (octave_idx_type i = 0; i < m; i++)
753 x(i, j) = zero;
754 }
755 }
756
757 return x;
758 }
759
760 // Right division functions.
761 //
762 // op2 / op1: dm cdm
763 // +-- +---+----+
764 // matrix | 1 | |
765 // +---+----+
766 // complex_matrix | 2 | 3 |
767 // +---+----+
768
769 // -*- 1 -*-
770 Matrix
771 xdiv (const Matrix& a, const DiagMatrix& b)
772 { return mdm_div_impl (a, b); }
773
774 // -*- 2 -*-
775 ComplexMatrix
776 xdiv (const ComplexMatrix& a, const DiagMatrix& b)
777 { return mdm_div_impl (a, b); }
778
779 // -*- 3 -*-
780 ComplexMatrix
781 xdiv (const ComplexMatrix& a, const ComplexDiagMatrix& b)
782 { return mdm_div_impl (a, b); }
783
784 // Right division functions, float type.
785 //
786 // op2 / op1: dm cdm
787 // +-- +---+----+
788 // matrix | 1 | |
789 // +---+----+
790 // complex_matrix | 2 | 3 |
791 // +---+----+
792
793 // -*- 1 -*-
794 FloatMatrix
795 xdiv (const FloatMatrix& a, const FloatDiagMatrix& b)
796 { return mdm_div_impl (a, b); }
797
798 // -*- 2 -*-
799 FloatComplexMatrix
800 xdiv (const FloatComplexMatrix& a, const FloatDiagMatrix& b)
801 { return mdm_div_impl (a, b); }
802
803 // -*- 3 -*-
804 FloatComplexMatrix
805 xdiv (const FloatComplexMatrix& a, const FloatComplexDiagMatrix& b)
806 { return mdm_div_impl (a, b); }
807
808 template <class MT, class DMT>
809 MT
810 dmm_leftdiv_impl (const DMT& d, const MT& a)
811 {
812 if (! mx_leftdiv_conform (d, a))
813 return MT ();
814
815 octave_idx_type m = d.cols (), n = a.cols (), k = d.rows ();
816 octave_idx_type mk = m < k ? m : k;
817 MT x (m, n);
818 const typename DMT::element_type zero = typename DMT::element_type ();
819
820 for (octave_idx_type j = 0; j < n; j++)
821 {
822 for (octave_idx_type i = 0; i < mk; i++)
823 x(i, j) = d(i, i) != zero ? a(i, j) / d(i, i) : 0;
824 for (octave_idx_type i = mk; i < m; i++)
825 x(i, j) = zero;
826 }
827
828 return x;
829 }
830
831 // Left division functions.
832 //
833 // op2 \ op1: m cm
834 // +---+----+
835 // diag_matrix | 1 | 2 |
836 // +---+----+
837 // complex_diag_matrix | | 3 |
838 // +---+----+
839
840 // -*- 1 -*-
841 Matrix
842 xleftdiv (const DiagMatrix& a, const Matrix& b)
843 { return dmm_leftdiv_impl (a, b); }
844
845 // -*- 2 -*-
846 ComplexMatrix
847 xleftdiv (const DiagMatrix& a, const ComplexMatrix& b)
848 { return dmm_leftdiv_impl (a, b); }
849
850 // -*- 3 -*-
851 ComplexMatrix
852 xleftdiv (const ComplexDiagMatrix& a, const ComplexMatrix& b)
853 { return dmm_leftdiv_impl (a, b); }
854
855 // Left division functions, float type.
856 //
857 // op2 \ op1: m cm
858 // +---+----+
859 // diag_matrix | 1 | 2 |
860 // +---+----+
861 // complex_diag_matrix | | 3 |
862 // +---+----+
863
864 // -*- 1 -*-
865 FloatMatrix
866 xleftdiv (const FloatDiagMatrix& a, const FloatMatrix& b)
867 { return dmm_leftdiv_impl (a, b); }
868
869 // -*- 2 -*-
870 FloatComplexMatrix
871 xleftdiv (const FloatDiagMatrix& a, const FloatComplexMatrix& b)
872 { return dmm_leftdiv_impl (a, b); }
873
874 // -*- 3 -*-
875 FloatComplexMatrix
876 xleftdiv (const FloatComplexDiagMatrix& a, const FloatComplexMatrix& b)
877 { return dmm_leftdiv_impl (a, b); }
878
879 // Diagonal by diagonal matrix division.
880
881 template <class MT, class DMT>
882 MT
883 dmdm_div_impl (const MT& a, const DMT& d)
884 {
885 if (! mx_div_conform (a, d))
886 return MT ();
887
888 octave_idx_type m = a.rows (), n = d.rows (), k = d.cols ();
889 octave_idx_type mn = m < n ? m : n;
890 MT x (m, n);
891 const typename DMT::element_type zero = typename DMT::element_type ();
892
893 for (octave_idx_type j = 0; j < mn; j++)
894 {
895 if (j < k && d(j, j) != zero)
896 x(j, j) = a(j, j) / d(j, j);
897 else
898 x(j, j) = zero;
899 }
900
901 return x;
902 }
903
904 // Right division functions.
905 //
906 // op2 / op1: dm cdm
907 // +-- +---+----+
908 // diag_matrix | 1 | |
909 // +---+----+
910 // complex_diag_matrix | 2 | 3 |
911 // +---+----+
912
913 // -*- 1 -*-
914 DiagMatrix
915 xdiv (const DiagMatrix& a, const DiagMatrix& b)
916 { return dmdm_div_impl (a, b); }
917
918 // -*- 2 -*-
919 ComplexDiagMatrix
920 xdiv (const ComplexDiagMatrix& a, const DiagMatrix& b)
921 { return dmdm_div_impl (a, b); }
922
923 // -*- 3 -*-
924 ComplexDiagMatrix
925 xdiv (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
926 { return dmdm_div_impl (a, b); }
927
928 // Right division functions, float type.
929 //
930 // op2 / op1: dm cdm
931 // +-- +---+----+
932 // diag_matrix | 1 | |
933 // +---+----+
934 // complex_diag_matrix | 2 | 3 |
935 // +---+----+
936
937 // -*- 1 -*-
938 FloatDiagMatrix
939 xdiv (const FloatDiagMatrix& a, const FloatDiagMatrix& b)
940 { return dmdm_div_impl (a, b); }
941
942 // -*- 2 -*-
943 FloatComplexDiagMatrix
944 xdiv (const FloatComplexDiagMatrix& a, const FloatDiagMatrix& b)
945 { return dmdm_div_impl (a, b); }
946
947 // -*- 3 -*-
948 FloatComplexDiagMatrix
949 xdiv (const FloatComplexDiagMatrix& a, const FloatComplexDiagMatrix& b)
950 { return dmdm_div_impl (a, b); }
951
952 template <class MT, class DMT>
953 MT
954 dmdm_leftdiv_impl (const DMT& d, const MT& a)
955 {
956 if (! mx_leftdiv_conform (d, a))
957 return MT ();
958
959 octave_idx_type m = d.cols (), n = a.cols (), k = d.rows ();
960 octave_idx_type mn = m < n ? m : n;
961 MT x (m, n);
962 const typename DMT::element_type zero = typename DMT::element_type ();
963
964 for (octave_idx_type j = 0; j < mn; j++)
965 {
966 if (j < k && d(j, j) != zero)
967 x(j, j) = a(j, j) / d(j, j);
968 else
969 x(j, j) = zero;
970 }
971
972 return x;
973 }
974
975 // Left division functions.
976 //
977 // op2 \ op1: dm cdm
978 // +---+----+
979 // diag_matrix | 1 | 2 |
980 // +---+----+
981 // complex_diag_matrix | | 3 |
982 // +---+----+
983
984 // -*- 1 -*-
985 DiagMatrix
986 xleftdiv (const DiagMatrix& a, const DiagMatrix& b)
987 { return dmdm_leftdiv_impl (a, b); }
988
989 // -*- 2 -*-
990 ComplexDiagMatrix
991 xleftdiv (const DiagMatrix& a, const ComplexDiagMatrix& b)
992 { return dmdm_leftdiv_impl (a, b); }
993
994 // -*- 3 -*-
995 ComplexDiagMatrix
996 xleftdiv (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
997 { return dmdm_leftdiv_impl (a, b); }
998
999 // Left division functions, float type.
1000 //
1001 // op2 \ op1: dm cdm
1002 // +---+----+
1003 // diag_matrix | 1 | 2 |
1004 // +---+----+
1005 // complex_diag_matrix | | 3 |
1006 // +---+----+
1007
1008 // -*- 1 -*-
1009 FloatDiagMatrix
1010 xleftdiv (const FloatDiagMatrix& a, const FloatDiagMatrix& b)
1011 { return dmdm_leftdiv_impl (a, b); }
1012
1013 // -*- 2 -*-
1014 FloatComplexDiagMatrix
1015 xleftdiv (const FloatDiagMatrix& a, const FloatComplexDiagMatrix& b)
1016 { return dmdm_leftdiv_impl (a, b); }
1017
1018 // -*- 3 -*-
1019 FloatComplexDiagMatrix
1020 xleftdiv (const FloatComplexDiagMatrix& a, const FloatComplexDiagMatrix& b)
1021 { return dmdm_leftdiv_impl (a, b); }
1022
725 /* 1023 /*
726 ;;; Local Variables: *** 1024 ;;; Local Variables: ***
727 ;;; mode: C++ *** 1025 ;;; mode: C++ ***
728 ;;; End: *** 1026 ;;; End: ***
729 */ 1027 */