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 "EIG.h"
00032 #include "dColVector.h"
00033 #include "f77-fcn.h"
00034 #include "lo-error.h"
00035
00036 extern "C"
00037 {
00038 F77_RET_T
00039 F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL,
00040 F77_CONST_CHAR_ARG_DECL,
00041 const int&, double*, const int&, double*,
00042 double*, double*, const int&, double*,
00043 const int&, double*, const int&, int&
00044 F77_CHAR_ARG_LEN_DECL
00045 F77_CHAR_ARG_LEN_DECL);
00046
00047 F77_RET_T
00048 F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL,
00049 F77_CONST_CHAR_ARG_DECL,
00050 const int&, Complex*, const int&, Complex*,
00051 Complex*, const int&, Complex*, const int&,
00052 Complex*, const int&, double*, int&
00053 F77_CHAR_ARG_LEN_DECL
00054 F77_CHAR_ARG_LEN_DECL);
00055
00056 F77_RET_T
00057 F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL,
00058 F77_CONST_CHAR_ARG_DECL,
00059 const int&, double*, const int&, double*,
00060 double*, const int&, int&
00061 F77_CHAR_ARG_LEN_DECL
00062 F77_CHAR_ARG_LEN_DECL);
00063
00064 F77_RET_T
00065 F77_FUNC (zheev, ZHEEV) (F77_CONST_CHAR_ARG_DECL,
00066 F77_CONST_CHAR_ARG_DECL,
00067 const int&, Complex*, const int&, double*,
00068 Complex*, const int&, double*, int&
00069 F77_CHAR_ARG_LEN_DECL
00070 F77_CHAR_ARG_LEN_DECL);
00071 }
00072
00073 int
00074 EIG::init (const Matrix& a, bool calc_ev)
00075 {
00076 if (a.is_symmetric ())
00077 return symmetric_init (a, calc_ev);
00078
00079 int n = a.rows ();
00080
00081 if (n != a.cols ())
00082 {
00083 (*current_liboctave_error_handler) ("EIG requires square matrix");
00084 return -1;
00085 }
00086
00087 int info = 0;
00088
00089 Matrix atmp = a;
00090 double *tmp_data = atmp.fortran_vec ();
00091
00092 Array<double> wr (n);
00093 double *pwr = wr.fortran_vec ();
00094
00095 Array<double> wi (n);
00096 double *pwi = wi.fortran_vec ();
00097
00098 volatile int nvr = calc_ev ? n : 0;
00099 Matrix vr (nvr, nvr);
00100 double *pvr = vr.fortran_vec ();
00101
00102 int lwork = -1;
00103 double dummy_work;
00104
00105 double *dummy = 0;
00106 int idummy = 1;
00107
00108 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
00109 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00110 n, tmp_data, n, pwr, pwi, dummy,
00111 idummy, pvr, n, &dummy_work, lwork, info
00112 F77_CHAR_ARG_LEN (1)
00113 F77_CHAR_ARG_LEN (1)));
00114
00115 if (! f77_exception_encountered && info == 0)
00116 {
00117 lwork = static_cast<int> (dummy_work);
00118 Array<double> work (lwork);
00119 double *pwork = work.fortran_vec ();
00120
00121 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
00122 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00123 n, tmp_data, n, pwr, pwi, dummy,
00124 idummy, pvr, n, pwork, lwork, info
00125 F77_CHAR_ARG_LEN (1)
00126 F77_CHAR_ARG_LEN (1)));
00127
00128 if (f77_exception_encountered || info < 0)
00129 {
00130 (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
00131 return info;
00132 }
00133
00134 if (info > 0)
00135 {
00136 (*current_liboctave_error_handler) ("dgeev failed to converge");
00137 return info;
00138 }
00139
00140 lambda.resize (n);
00141 v.resize (nvr, nvr);
00142
00143 for (int j = 0; j < n; j++)
00144 {
00145 if (wi.elem (j) == 0.0)
00146 {
00147 lambda.elem (j) = Complex (wr.elem (j));
00148 for (int i = 0; i < nvr; i++)
00149 v.elem (i, j) = vr.elem (i, j);
00150 }
00151 else
00152 {
00153 if (j+1 >= n)
00154 {
00155 (*current_liboctave_error_handler) ("EIG: internal error");
00156 return -1;
00157 }
00158
00159 lambda.elem(j) = Complex (wr.elem(j), wi.elem(j));
00160 lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1));
00161
00162 for (int i = 0; i < nvr; i++)
00163 {
00164 double real_part = vr.elem (i, j);
00165 double imag_part = vr.elem (i, j+1);
00166 v.elem (i, j) = Complex (real_part, imag_part);
00167 v.elem (i, j+1) = Complex (real_part, -imag_part);
00168 }
00169 j++;
00170 }
00171 }
00172 }
00173 else
00174 (*current_liboctave_error_handler) ("dgeev workspace query failed");
00175
00176 return info;
00177 }
00178
00179 int
00180 EIG::symmetric_init (const Matrix& a, bool calc_ev)
00181 {
00182 int n = a.rows ();
00183
00184 if (n != a.cols ())
00185 {
00186 (*current_liboctave_error_handler) ("EIG requires square matrix");
00187 return -1;
00188 }
00189
00190 int info = 0;
00191
00192 Matrix atmp = a;
00193 double *tmp_data = atmp.fortran_vec ();
00194
00195 ColumnVector wr (n);
00196 double *pwr = wr.fortran_vec ();
00197
00198 int lwork = -1;
00199 double dummy_work;
00200
00201 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00202 F77_CONST_CHAR_ARG2 ("U", 1),
00203 n, tmp_data, n, pwr, &dummy_work, lwork, info
00204 F77_CHAR_ARG_LEN (1)
00205 F77_CHAR_ARG_LEN (1)));
00206
00207 if (! f77_exception_encountered && info == 0)
00208 {
00209 lwork = static_cast<int> (dummy_work);
00210 Array<double> work (lwork);
00211 double *pwork = work.fortran_vec ();
00212
00213 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00214 F77_CONST_CHAR_ARG2 ("U", 1),
00215 n, tmp_data, n, pwr, pwork, lwork, info
00216 F77_CHAR_ARG_LEN (1)
00217 F77_CHAR_ARG_LEN (1)));
00218
00219 if (f77_exception_encountered || info < 0)
00220 {
00221 (*current_liboctave_error_handler) ("unrecoverable error in dsyev");
00222 return info;
00223 }
00224
00225 if (info > 0)
00226 {
00227 (*current_liboctave_error_handler) ("dsyev failed to converge");
00228 return info;
00229 }
00230
00231 lambda = ComplexColumnVector (wr);
00232 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
00233 }
00234 else
00235 (*current_liboctave_error_handler) ("dsyev workspace query failed");
00236
00237 return info;
00238 }
00239
00240 int
00241 EIG::init (const ComplexMatrix& a, bool calc_ev)
00242 {
00243 if (a.is_hermitian ())
00244 return hermitian_init (a, calc_ev);
00245
00246 int n = a.rows ();
00247
00248 if (n != a.cols ())
00249 {
00250 (*current_liboctave_error_handler) ("EIG requires square matrix");
00251 return -1;
00252 }
00253
00254 int info = 0;
00255
00256 ComplexMatrix atmp = a;
00257 Complex *tmp_data = atmp.fortran_vec ();
00258
00259 ComplexColumnVector w (n);
00260 Complex *pw = w.fortran_vec ();
00261
00262 int nvr = calc_ev ? n : 0;
00263 ComplexMatrix vtmp (nvr, nvr);
00264 Complex *pv = vtmp.fortran_vec ();
00265
00266 int lwork = -1;
00267 Complex dummy_work;
00268
00269 int lrwork = 2*n;
00270 Array<double> rwork (lrwork);
00271 double *prwork = rwork.fortran_vec ();
00272
00273 Complex *dummy = 0;
00274 int idummy = 1;
00275
00276 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
00277 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00278 n, tmp_data, n, pw, dummy, idummy,
00279 pv, n, &dummy_work, lwork, prwork, info
00280 F77_CHAR_ARG_LEN (1)
00281 F77_CHAR_ARG_LEN (1)));
00282
00283 if (! f77_exception_encountered && info == 0)
00284 {
00285 lwork = static_cast<int> (dummy_work.real ());
00286 Array<Complex> work (lwork);
00287 Complex *pwork = work.fortran_vec ();
00288
00289 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
00290 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00291 n, tmp_data, n, pw, dummy, idummy,
00292 pv, n, pwork, lwork, prwork, info
00293 F77_CHAR_ARG_LEN (1)
00294 F77_CHAR_ARG_LEN (1)));
00295
00296 if (f77_exception_encountered || info < 0)
00297 {
00298 (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
00299 return info;
00300 }
00301
00302 if (info > 0)
00303 {
00304 (*current_liboctave_error_handler) ("zgeev failed to converge");
00305 return info;
00306 }
00307
00308 lambda = w;
00309 v = vtmp;
00310 }
00311 else
00312 (*current_liboctave_error_handler) ("zgeev workspace query failed");
00313
00314 return info;
00315 }
00316
00317 int
00318 EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev)
00319 {
00320 int n = a.rows ();
00321
00322 if (n != a.cols ())
00323 {
00324 (*current_liboctave_error_handler) ("EIG requires square matrix");
00325 return -1;
00326 }
00327
00328 int info = 0;
00329
00330 ComplexMatrix atmp = a;
00331 Complex *tmp_data = atmp.fortran_vec ();
00332
00333 ColumnVector wr (n);
00334 double *pwr = wr.fortran_vec ();
00335
00336 int lwork = -1;
00337 Complex dummy_work;
00338
00339 int lrwork = 3*n;
00340 Array<double> rwork (lrwork);
00341 double *prwork = rwork.fortran_vec ();
00342
00343 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00344 F77_CONST_CHAR_ARG2 ("U", 1),
00345 n, tmp_data, n, pwr, &dummy_work, lwork,
00346 prwork, info
00347 F77_CHAR_ARG_LEN (1)
00348 F77_CHAR_ARG_LEN (1)));
00349
00350 if (! f77_exception_encountered && info == 0)
00351 {
00352 lwork = static_cast<int> (dummy_work.real ());
00353 Array<Complex> work (lwork);
00354 Complex *pwork = work.fortran_vec ();
00355
00356 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
00357 F77_CONST_CHAR_ARG2 ("U", 1),
00358 n, tmp_data, n, pwr, pwork, lwork, prwork, info
00359 F77_CHAR_ARG_LEN (1)
00360 F77_CHAR_ARG_LEN (1)));
00361
00362 if (f77_exception_encountered || info < 0)
00363 {
00364 (*current_liboctave_error_handler) ("unrecoverable error in zheev");
00365 return info;
00366 }
00367
00368 if (info > 0)
00369 {
00370 (*current_liboctave_error_handler) ("zheev failed to converge");
00371 return info;
00372 }
00373
00374 lambda = ComplexColumnVector (wr);
00375 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
00376 }
00377 else
00378 (*current_liboctave_error_handler) ("zheev workspace query failed");
00379
00380 return info;
00381 }
00382
00383
00384
00385
00386
00387