00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include <iostream>
00032
00033 #include "dbleSCHUR.h"
00034 #include "f77-fcn.h"
00035 #include "lo-error.h"
00036
00037 extern "C"
00038 {
00039 F77_RET_T
00040 F77_FUNC (dgeesx, DGEESX) (F77_CONST_CHAR_ARG_DECL,
00041 F77_CONST_CHAR_ARG_DECL,
00042 SCHUR::select_function,
00043 F77_CONST_CHAR_ARG_DECL,
00044 const int&, double*, const int&, int&,
00045 double*, double*, double*, const int&,
00046 double&, double&, double*, const int&,
00047 int*, const int&, int*, int&
00048 F77_CHAR_ARG_LEN_DECL
00049 F77_CHAR_ARG_LEN_DECL
00050 F77_CHAR_ARG_LEN_DECL);
00051 }
00052
00053 static int
00054 select_ana (const double& a, const double&)
00055 {
00056 return (a < 0.0);
00057 }
00058
00059 static int
00060 select_dig (const double& a, const double& b)
00061 {
00062 return (hypot (a, b) < 1.0);
00063 }
00064
00065 int
00066 SCHUR::init (const Matrix& a, const std::string& ord, bool calc_unitary)
00067 {
00068 int a_nr = a.rows ();
00069 int a_nc = a.cols ();
00070
00071 if (a_nr != a_nc)
00072 {
00073 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
00074 return -1;
00075 }
00076
00077
00078
00079
00080 char jobvs;
00081 char sense = 'N';
00082 char sort = 'N';
00083
00084 if (calc_unitary)
00085 jobvs = 'V';
00086 else
00087 jobvs = 'N';
00088
00089 char ord_char = ord.empty () ? 'U' : ord[0];
00090
00091 if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
00092 sort = 'S';
00093
00094 if (ord_char == 'A' || ord_char == 'a')
00095 selector = select_ana;
00096 else if (ord_char == 'D' || ord_char == 'd')
00097 selector = select_dig;
00098 else
00099 selector = 0;
00100
00101 int n = a_nc;
00102 int lwork = 8 * n;
00103 int liwork = 1;
00104 int info;
00105 int sdim;
00106 double rconde;
00107 double rcondv;
00108
00109 schur_mat = a;
00110
00111 if (calc_unitary)
00112 unitary_mat.resize (n, n);
00113
00114 double *s = schur_mat.fortran_vec ();
00115 double *q = unitary_mat.fortran_vec ();
00116
00117 Array<double> wr (n);
00118 double *pwr = wr.fortran_vec ();
00119
00120 Array<double> wi (n);
00121 double *pwi = wi.fortran_vec ();
00122
00123 Array<double> work (lwork);
00124 double *pwork = work.fortran_vec ();
00125
00126
00127 Array<int> bwork ((ord_char == 'N' || ord_char == 'n') ? 0 : n);
00128 int *pbwork = bwork.fortran_vec ();
00129
00130 Array<int> iwork (liwork);
00131 int *piwork = iwork.fortran_vec ();
00132
00133 F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
00134 F77_CONST_CHAR_ARG2 (&sort, 1),
00135 selector,
00136 F77_CONST_CHAR_ARG2 (&sense, 1),
00137 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
00138 pwork, lwork, piwork, liwork, pbwork, info
00139 F77_CHAR_ARG_LEN (1)
00140 F77_CHAR_ARG_LEN (1)
00141 F77_CHAR_ARG_LEN (1)));
00142
00143 if (f77_exception_encountered)
00144 (*current_liboctave_error_handler) ("unrecoverable error in dgeesx");
00145
00146 return info;
00147 }
00148
00149 std::ostream&
00150 operator << (std::ostream& os, const SCHUR& a)
00151 {
00152 os << a.schur_matrix () << "\n";
00153 os << a.unitary_matrix () << "\n";
00154
00155 return os;
00156 }
00157
00158
00159
00160
00161
00162