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 "CmplxHESS.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034
00035 extern "C"
00036 {
00037 F77_RET_T
00038 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
00039 const int&, Complex*, const int&,
00040 int&, int&, double*, int&
00041 F77_CHAR_ARG_LEN_DECL);
00042
00043 F77_RET_T
00044 F77_FUNC (zgehrd, ZGEHRD) (const int&, const int&, const int&,
00045 Complex*, const int&, Complex*,
00046 Complex*, const int&, int&);
00047
00048 F77_RET_T
00049 F77_FUNC (zunghr, ZUNGHR) (const int&, const int&, const int&,
00050 Complex*, const int&, Complex*,
00051 Complex*, const int&, int&);
00052
00053 F77_RET_T
00054 F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL,
00055 F77_CONST_CHAR_ARG_DECL,
00056 const int&, const int&, const int&, double*,
00057 const int&, Complex*, const int&, int&
00058 F77_CHAR_ARG_LEN_DECL
00059 F77_CHAR_ARG_LEN_DECL);
00060 }
00061
00062 int
00063 ComplexHESS::init (const ComplexMatrix& a)
00064 {
00065 int a_nr = a.rows ();
00066 int a_nc = a.cols ();
00067
00068 if (a_nr != a_nc)
00069 {
00070 (*current_liboctave_error_handler)
00071 ("ComplexHESS requires square matrix");
00072 return -1;
00073 }
00074
00075 char job = 'N';
00076 char side = 'R';
00077
00078 int n = a_nc;
00079 int lwork = 32 * n;
00080 int info;
00081 int ilo;
00082 int ihi;
00083
00084 hess_mat = a;
00085 Complex *h = hess_mat.fortran_vec ();
00086
00087 Array<double> scale (n);
00088 double *pscale = scale.fortran_vec ();
00089
00090 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
00091 n, h, n, ilo, ihi, pscale, info
00092 F77_CHAR_ARG_LEN (1)));
00093
00094 if (f77_exception_encountered)
00095 (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
00096 else
00097 {
00098 Array<Complex> tau (n-1);
00099 Complex *ptau = tau.fortran_vec ();
00100
00101 Array<Complex> work (lwork);
00102 Complex *pwork = work.fortran_vec ();
00103
00104 F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
00105
00106 if (f77_exception_encountered)
00107 (*current_liboctave_error_handler) ("unrecoverable error in zgehrd");
00108 else
00109 {
00110 unitary_hess_mat = hess_mat;
00111 Complex *z = unitary_hess_mat.fortran_vec ();
00112
00113 F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
00114 lwork, info));
00115
00116 if (f77_exception_encountered)
00117 (*current_liboctave_error_handler)
00118 ("unrecoverable error in zunghr");
00119 else
00120 {
00121 F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
00122 F77_CONST_CHAR_ARG2 (&side, 1),
00123 n, ilo, ihi, pscale, n, z, n, info
00124 F77_CHAR_ARG_LEN (1)
00125 F77_CHAR_ARG_LEN (1)));
00126
00127 if (f77_exception_encountered)
00128 (*current_liboctave_error_handler)
00129 ("unrecoverable error in zgebak");
00130 else
00131 {
00132
00133
00134
00135
00136 if (n > 2)
00137 for (int j = 0; j < a_nc; j++)
00138 for (int i = j+2; i < a_nr; i++)
00139 hess_mat.elem (i, j) = 0;
00140 }
00141 }
00142 }
00143 }
00144
00145 return info;
00146 }
00147
00148
00149
00150
00151
00152