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