Mercurial > hg > octave-jordi

comparison libinterp/corefcn/qz.cc @ 20830:35241c4b696c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
eliminate return statements after calls to error * Cell.cc, __ichol__.cc, __lin_interpn__.cc, __pchip_deriv__.cc, besselj.cc, cellfun.cc, colloc.cc, debug.cc, dlmread.cc, dynamic-ld.cc, filter.cc, find.cc, gl2ps-renderer.cc, load-path.cc, load-save.cc, ls-mat4.cc, ls-mat5.cc, ls-oct-text.cc, luinc.cc, max.cc, nproc.cc, oct-hist.cc, oct-map.cc, oct-obj.cc, oct-stream.cc, ordschur.cc, pinv.cc, pr-output.cc, profiler.cc, psi.cc, quadcc.cc, qz.cc, rand.cc, strfind.cc, strfns.cc, sysdep.cc, toplev.cc, tril.cc, typecast.cc, urlwrite.cc, utils.cc, variables.cc: Eliminate return statements after calls to error.
author John W. Eaton <jwe@octave.org>
date Wed, 09 Dec 2015 14:00:43 -0500
parents 8bb38ba1bad6
children 1142cf6abc0d
comparison
equal deleted inserted replaced
20829:ae0bd73671f3 20830:35241c4b696c
407 407
408 if (! (ord_job == 'N' || ord_job == 'n' 408 if (! (ord_job == 'N' || ord_job == 'n'
409 || ord_job == 'S' || ord_job == 's' 409 || ord_job == 'S' || ord_job == 's'
410 || ord_job == 'B' || ord_job == 'b' 410 || ord_job == 'B' || ord_job == 'b'
411 || ord_job == '+' || ord_job == '-')) 411 || ord_job == '+' || ord_job == '-'))
412 { 412 error ("qz: invalid order option");
413 error ("qz: invalid order option");
414 return retval;
415 }
416 413
417 // overflow constant required by dlag2 414 // overflow constant required by dlag2
418 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1), 415 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
419 safmin 416 safmin
420 F77_CHAR_ARG_LEN (1)); 417 F77_CHAR_ARG_LEN (1));
510 507
511 volatile int complex_case 508 volatile int complex_case
512 = (args(0).is_complex_type () || args(1).is_complex_type ()); 509 = (args(0).is_complex_type () || args(1).is_complex_type ());
513 510
514 if (nargin == 3 && complex_case) 511 if (nargin == 3 && complex_case)
515 { 512 error ("qz: cannot re-order complex qz decomposition");
516 error ("qz: cannot re-order complex qz decomposition");
517 return retval;
518 }
519 513
520 // First, declare variables used in both the real and complex case. 514 // First, declare variables used in both the real and complex case.
521 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); 515 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
522 RowVector alphar(nn), alphai(nn), betar(nn); 516 RowVector alphar(nn), alphai(nn), betar(nn);
523 ComplexRowVector xalpha(nn), xbeta(nn); 517 ComplexRowVector xalpha(nn), xbeta(nn);
792 786
793 // Order the QZ decomposition? 787 // Order the QZ decomposition?
794 if (! (ord_job == 'N' || ord_job == 'n')) 788 if (! (ord_job == 'N' || ord_job == 'n'))
795 { 789 {
796 if (complex_case) 790 if (complex_case)
797 { 791 // Probably not needed, but better be safe.
798 // Probably not needed, but better be safe. 792 error ("qz: cannot re-order complex qz decomposition");
799 error ("qz: cannot re-order complex qz decomposition"); 793
800 return retval;
801 }
802 else
803 {
804 #ifdef DEBUG_SORT 794 #ifdef DEBUG_SORT
805 std::cout << "qz: ordering eigenvalues: ord_job = " 795 std::cout << "qz: ordering eigenvalues: ord_job = "
806 << ord_job << std::endl; 796 << ord_job << std::endl;
807 #endif 797 #endif
808 798
809 // Declared static to avoid vfork/long jump compiler complaints. 799 // Declared static to avoid vfork/long jump compiler complaints.
810 static sort_function sort_test; 800 static sort_function sort_test;
811 sort_test = 0; 801 sort_test = 0;
812 802
813 switch (ord_job) 803 switch (ord_job)
804 {
805 case 'S':
806 case 's':
807 sort_test = &fin;
808 break;
809
810 case 'B':
811 case 'b':
812 sort_test = &fout;
813 break;
814
815 case '+':
816 sort_test = &fcrhp;
817 break;
818
819 case '-':
820 sort_test = &folhp;
821 break;
822
823 default:
824 // Invalid order option (should never happen, since we
825 // checked the options at the top).
826 panic_impossible ();
827 break;
828 }
829
830 octave_idx_type ndim, fail;
831 double inf_norm;
832
833 F77_XFCN (xdlange, XDLANGE,
834 (F77_CONST_CHAR_ARG2 ("I", 1),
835 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
836 F77_CHAR_ARG_LEN (1)));
837
838 double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;
839
840 #ifdef DEBUG_SORT
841 std::cout << "qz: calling dsubsp: aa=" << std::endl;
842 octave_print_internal (std::cout, aa, 0);
843 std::cout << std::endl << "bb=" << std::endl;
844 octave_print_internal (std::cout, bb, 0);
845 if (compz == 'V')
846 {
847 std::cout << std::endl << "ZZ=" << std::endl;
848 octave_print_internal (std::cout, ZZ, 0);
849 }
850 std::cout << std::endl;
851 std::cout << "alphar = " << std::endl;
852 octave_print_internal (std::cout, (Matrix) alphar, 0);
853 std::cout << std::endl << "alphai = " << std::endl;
854 octave_print_internal (std::cout, (Matrix) alphai, 0);
855 std::cout << std::endl << "beta = " << std::endl;
856 octave_print_internal (std::cout, (Matrix) betar, 0);
857 std::cout << std::endl;
858 #endif
859
860 Array<octave_idx_type> ind (dim_vector (nn, 1));
861
862 F77_XFCN (dsubsp, DSUBSP,
863 (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
864 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
865 ind.fortran_vec ()));
866
867 #ifdef DEBUG
868 std::cout << "qz: back from dsubsp: aa=" << std::endl;
869 octave_print_internal (std::cout, aa, 0);
870 std::cout << std::endl << "bb=" << std::endl;
871 octave_print_internal (std::cout, bb, 0);
872 if (compz == 'V')
873 {
874 std::cout << std::endl << "ZZ=" << std::endl;
875 octave_print_internal (std::cout, ZZ, 0);
876 }
877 std::cout << std::endl;
878 #endif
879
880 // Manually update alphar, alphai, betar.
881 static int jj;
882
883 jj = 0;
884 while (jj < nn)
885 {
886 #ifdef DEBUG_EIG
887 std::cout << "computing gen eig #" << jj << std::endl;
888 #endif
889
890 // Number of zeros in this block.
891 static int zcnt;
892
893 if (jj == (nn-1))
894 zcnt = 1;
895 else if (aa(jj+1,jj) == 0)
896 zcnt = 1;
897 else zcnt = 2;
898
899 if (zcnt == 1)
814 { 900 {
815 case 'S': 901 // Real zero.
816 case 's': 902 #ifdef DEBUG_EIG
817 sort_test = &fin; 903 std::cout << " single gen eig:" << std::endl;
818 break; 904 std::cout << " alphar(" << jj << ") = " << aa(jj,jj)
819 905 << std::endl;
820 case 'B': 906 std::cout << " betar(" << jj << ") = " << bb(jj,jj)
821 case 'b': 907 << std::endl;
822 sort_test = &fout; 908 std::cout << " alphai(" << jj << ") = 0" << std::endl;
823 break; 909 #endif
824 910
825 case '+': 911 alphar(jj) = aa(jj,jj);
826 sort_test = &fcrhp; 912 alphai(jj) = 0;
827 break; 913 betar(jj) = bb(jj,jj);
828
829 case '-':
830 sort_test = &folhp;
831 break;
832
833 default:
834 // Invalid order option (should never happen, since we
835 // checked the options at the top).
836 panic_impossible ();
837 break;
838 } 914 }
839 915 else
840 octave_idx_type ndim, fail;
841 double inf_norm;
842
843 F77_XFCN (xdlange, XDLANGE,
844 (F77_CONST_CHAR_ARG2 ("I", 1),
845 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
846 F77_CHAR_ARG_LEN (1)));
847
848 double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn;
849
850 #ifdef DEBUG_SORT
851 std::cout << "qz: calling dsubsp: aa=" << std::endl;
852 octave_print_internal (std::cout, aa, 0);
853 std::cout << std::endl << "bb=" << std::endl;
854 octave_print_internal (std::cout, bb, 0);
855 if (compz == 'V')
856 { 916 {
857 std::cout << std::endl << "ZZ=" << std::endl; 917 // Complex conjugate pair.
858 octave_print_internal (std::cout, ZZ, 0);
859 }
860 std::cout << std::endl;
861 std::cout << "alphar = " << std::endl;
862 octave_print_internal (std::cout, (Matrix) alphar, 0);
863 std::cout << std::endl << "alphai = " << std::endl;
864 octave_print_internal (std::cout, (Matrix) alphai, 0);
865 std::cout << std::endl << "beta = " << std::endl;
866 octave_print_internal (std::cout, (Matrix) betar, 0);
867 std::cout << std::endl;
868 #endif
869
870 Array<octave_idx_type> ind (dim_vector (nn, 1));
871
872 F77_XFCN (dsubsp, DSUBSP,
873 (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
874 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
875 ind.fortran_vec ()));
876
877 #ifdef DEBUG
878 std::cout << "qz: back from dsubsp: aa=" << std::endl;
879 octave_print_internal (std::cout, aa, 0);
880 std::cout << std::endl << "bb=" << std::endl;
881 octave_print_internal (std::cout, bb, 0);
882 if (compz == 'V')
883 {
884 std::cout << std::endl << "ZZ=" << std::endl;
885 octave_print_internal (std::cout, ZZ, 0);
886 }
887 std::cout << std::endl;
888 #endif
889
890 // Manually update alphar, alphai, betar.
891 static int jj;
892
893 jj = 0;
894 while (jj < nn)
895 {
896 #ifdef DEBUG_EIG 918 #ifdef DEBUG_EIG
897 std::cout << "computing gen eig #" << jj << std::endl; 919 std::cout << "qz: calling dlag2:" << std::endl;
898 #endif 920 std::cout << "safmin="
899 921 << setiosflags (std::ios::scientific)
900 // Number of zeros in this block. 922 << safmin << std::endl;
901 static int zcnt; 923
902 924 for (int idr = jj; idr <= jj+1; idr++)
903 if (jj == (nn-1))
904 zcnt = 1;
905 else if (aa(jj+1,jj) == 0)
906 zcnt = 1;
907 else zcnt = 2;
908
909 if (zcnt == 1)
910 { 925 {
911 // Real zero. 926 for (int idc = jj; idc <= jj+1; idc++)
927 {
928 std::cout << "aa(" << idr << "," << idc << ")="
929 << aa(idr,idc) << std::endl;
930 std::cout << "bb(" << idr << "," << idc << ")="
931 << bb(idr,idc) << std::endl;
932 }
933 }
934 #endif
935
936 // FIXME: probably should be using
937 // fortran_vec instead of &aa(jj,jj) here.
938
939 double scale1, scale2, wr1, wr2, wi;
940 const double *aa_ptr = aa.data () + jj * nn + jj;
941 const double *bb_ptr = bb.data () + jj * nn + jj;
942 F77_XFCN (dlag2, DLAG2,
943 (aa_ptr, nn, bb_ptr, nn, safmin,
944 scale1, scale2, wr1, wr2, wi));
945
912 #ifdef DEBUG_EIG 946 #ifdef DEBUG_EIG
913 std::cout << " single gen eig:" << std::endl; 947 std::cout << "dlag2 returns: scale1=" << scale1
914 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) 948 << "\tscale2=" << scale2 << std::endl
915 << std::endl; 949 << "\twr1=" << wr1 << "\twr2=" << wr2
916 std::cout << " betar(" << jj << ") = " << bb(jj,jj) 950 << "\twi=" << wi << std::endl;
917 << std::endl; 951 #endif
918 std::cout << " alphai(" << jj << ") = 0" << std::endl; 952
919 #endif 953 // Just to be safe, check if it's a real pair.
920 954 if (wi == 0)
921 alphar(jj) = aa(jj,jj); 955 {
956 alphar(jj) = wr1;
922 alphai(jj) = 0; 957 alphai(jj) = 0;
923 betar(jj) = bb(jj,jj); 958 betar(jj) = scale1;
959 alphar(jj+1) = wr2;
960 alphai(jj+1) = 0;
961 betar(jj+1) = scale2;
924 } 962 }
925 else 963 else
926 { 964 {
927 // Complex conjugate pair. 965 alphar(jj) = alphar(jj+1) = wr1;
928 #ifdef DEBUG_EIG 966 alphai(jj) = -(alphai(jj+1) = wi);
929 std::cout << "qz: calling dlag2:" << std::endl; 967 betar(jj) = betar(jj+1) = scale1;
930 std::cout << "safmin="
931 << setiosflags (std::ios::scientific)
932 << safmin << std::endl;
933
934 for (int idr = jj; idr <= jj+1; idr++)
935 {
936 for (int idc = jj; idc <= jj+1; idc++)
937 {
938 std::cout << "aa(" << idr << "," << idc << ")="
939 << aa(idr,idc) << std::endl;
940 std::cout << "bb(" << idr << "," << idc << ")="
941 << bb(idr,idc) << std::endl;
942 }
943 }
944 #endif
945
946 // FIXME: probably should be using
947 // fortran_vec instead of &aa(jj,jj) here.
948
949 double scale1, scale2, wr1, wr2, wi;
950 const double *aa_ptr = aa.data () + jj * nn + jj;
951 const double *bb_ptr = bb.data () + jj * nn + jj;
952 F77_XFCN (dlag2, DLAG2,
953 (aa_ptr, nn, bb_ptr, nn, safmin,
954 scale1, scale2, wr1, wr2, wi));
955
956 #ifdef DEBUG_EIG
957 std::cout << "dlag2 returns: scale1=" << scale1
958 << "\tscale2=" << scale2 << std::endl
959 << "\twr1=" << wr1 << "\twr2=" << wr2
960 << "\twi=" << wi << std::endl;
961 #endif
962
963 // Just to be safe, check if it's a real pair.
964 if (wi == 0)
965 {
966 alphar(jj) = wr1;
967 alphai(jj) = 0;
968 betar(jj) = scale1;
969 alphar(jj+1) = wr2;
970 alphai(jj+1) = 0;
971 betar(jj+1) = scale2;
972 }
973 else
974 {
975 alphar(jj) = alphar(jj+1) = wr1;
976 alphai(jj) = -(alphai(jj+1) = wi);
977 betar(jj) = betar(jj+1) = scale1;
978 }
979 } 968 }
980
981 // Advance past this block.
982 jj += zcnt;
983 } 969 }
984 970
971 // Advance past this block.
972 jj += zcnt;
973 }
974
985 #ifdef DEBUG_SORT 975 #ifdef DEBUG_SORT
986 std::cout << "qz: back from dsubsp: aa=" << std::endl; 976 std::cout << "qz: back from dsubsp: aa=" << std::endl;
987 octave_print_internal (std::cout, aa, 0); 977 octave_print_internal (std::cout, aa, 0);
988 std::cout << std::endl << "bb=" << std::endl; 978 std::cout << std::endl << "bb=" << std::endl;
989 octave_print_internal (std::cout, bb, 0); 979 octave_print_internal (std::cout, bb, 0);
990 980
991 if (compz == 'V') 981 if (compz == 'V')
992 { 982 {
993 std::cout << std::endl << "ZZ=" << std::endl; 983 std::cout << std::endl << "ZZ=" << std::endl;
994 octave_print_internal (std::cout, ZZ, 0); 984 octave_print_internal (std::cout, ZZ, 0);
995 } 985 }
996 std::cout << std::endl << "qz: ndim=" << ndim << std::endl 986 std::cout << std::endl << "qz: ndim=" << ndim << std::endl
997 << "fail=" << fail << std::endl; 987 << "fail=" << fail << std::endl;
998 std::cout << "alphar = " << std::endl; 988 std::cout << "alphar = " << std::endl;
999 octave_print_internal (std::cout, (Matrix) alphar, 0); 989 octave_print_internal (std::cout, (Matrix) alphar, 0);
1000 std::cout << std::endl << "alphai = " << std::endl; 990 std::cout << std::endl << "alphai = " << std::endl;
1001 octave_print_internal (std::cout, (Matrix) alphai, 0); 991 octave_print_internal (std::cout, (Matrix) alphai, 0);
1002 std::cout << std::endl << "beta = " << std::endl; 992 std::cout << std::endl << "beta = " << std::endl;
1003 octave_print_internal (std::cout, (Matrix) betar, 0); 993 octave_print_internal (std::cout, (Matrix) betar, 0);
1004 std::cout << std::endl; 994 std::cout << std::endl;
1005 #endif 995 #endif
1006 }
1007 } 996 }
1008 997
1009 // Compute generalized eigenvalues? 998 // Compute generalized eigenvalues?
1010 ComplexColumnVector gev; 999 ComplexColumnVector gev;
1011 1000