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 "CmplxSCHUR.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034
00035 extern "C"
00036 {
00037 F77_RET_T
00038 F77_FUNC (zgeesx, ZGEESX) (F77_CONST_CHAR_ARG_DECL,
00039 F77_CONST_CHAR_ARG_DECL,
00040 ComplexSCHUR::select_function,
00041 F77_CONST_CHAR_ARG_DECL,
00042 const int&, Complex*, const int&, int&,
00043 Complex*, Complex*, const int&, double&,
00044 double&, Complex*, const int&, double*, int*,
00045 int&
00046 F77_CHAR_ARG_LEN_DECL
00047 F77_CHAR_ARG_LEN_DECL
00048 F77_CHAR_ARG_LEN_DECL);
00049 }
00050
00051 static int
00052 select_ana (const Complex& a)
00053 {
00054 return a.real () < 0.0;
00055 }
00056
00057 static int
00058 select_dig (const Complex& a)
00059 {
00060 return (abs (a) < 1.0);
00061 }
00062
00063 int
00064 ComplexSCHUR::init (const ComplexMatrix& a, const std::string& ord,
00065 bool calc_unitary)
00066 {
00067 int a_nr = a.rows ();
00068 int a_nc = a.cols ();
00069
00070 if (a_nr != a_nc)
00071 {
00072 (*current_liboctave_error_handler)
00073 ("ComplexSCHUR 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 info;
00104 int sdim;
00105 double rconde;
00106 double rcondv;
00107
00108 schur_mat = a;
00109 if (calc_unitary)
00110 unitary_mat.resize (n, n);
00111
00112 Complex *s = schur_mat.fortran_vec ();
00113 Complex *q = unitary_mat.fortran_vec ();
00114
00115 Array<double> rwork (n);
00116 double *prwork = rwork.fortran_vec ();
00117
00118 Array<Complex> w (n);
00119 Complex *pw = w.fortran_vec ();
00120
00121 Array<Complex> work (lwork);
00122 Complex *pwork = work.fortran_vec ();
00123
00124
00125 Array<int> bwork ((ord_char == 'N' || ord_char == 'n') ? 0 : n);
00126 int *pbwork = bwork.fortran_vec ();
00127
00128 F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
00129 F77_CONST_CHAR_ARG2 (&sort, 1),
00130 selector,
00131 F77_CONST_CHAR_ARG2 (&sense, 1),
00132 n, s, n, sdim, pw, q, n, rconde, rcondv,
00133 pwork, lwork, prwork, pbwork, info
00134 F77_CHAR_ARG_LEN (1)
00135 F77_CHAR_ARG_LEN (1)
00136 F77_CHAR_ARG_LEN (1)));
00137
00138 if (f77_exception_encountered)
00139 (*current_liboctave_error_handler) ("unrecoverable error in zgeesx");
00140
00141 return info;
00142 }
00143
00144
00145
00146
00147
00148