API tools faq
paste
Advertisement
BumiBarbi

pb.diff

BumiBarbi
Mar 29th, 2016
112
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Diff 28.36 KB | None | 0 0
raw download clone embed print report
  1. # HG changeset patch
  2. # User Barbara Locsi <locsi.barbara@gmail.com>
  3. # Date 1459272054 -7200
  4. #      Tue Mar 29 19:20:54 2016 +0200
  5. # Node ID c6383dae0bd5cabfd35647e26a67092cb19910fe
  6. # Parent  102b33b53ea466d7e87f4fed51f4b9b7d05b410e
  7. preliminary balancing for eig
  8.  
  9. diff -r 102b33b53ea4 -r c6383dae0bd5 libinterp/corefcn/eig.cc
  10. --- a/libinterp/corefcn/eig.cc  Mon Mar 28 20:00:37 2016 +1100
  11. +++ b/libinterp/corefcn/eig.cc  Tue Mar 29 19:20:54 2016 +0200
  12. @@ -54,7 +54,7 @@ The eigenvalues returned by @code{eig} a
  13.  
  14.    if (nargin > 2 || nargin == 0)
  15.      print_usage ();
  16. -
  17. +  
  18.    octave_value_list retval;
  19.  
  20.    octave_value arg_a, arg_b;
  21. @@ -75,26 +75,40 @@ The eigenvalues returned by @code{eig} a
  22.    if (! arg_a.is_double_type () && ! arg_a.is_single_type ())
  23.      err_wrong_type_arg ("eig", arg_a);
  24.  
  25. +  bool balance = true;
  26.    if (nargin == 2)
  27.      {
  28. -      arg_b = args(1);
  29. -      nr_b = arg_b.rows ();
  30. -      nc_b = arg_b.columns ();
  31. +      if(args(1).is_string ())
  32. +        {
  33. +          std::string a1s = args(1).string_value ();
  34. +          if(a1s != "balance" && a1s != "nobalance")
  35. +            error ("eig: unexpected second input: %s", a1s.c_str ());
  36. +          balance = a1s == "balance";
  37. +        }
  38. +      else
  39. +        {
  40. +          arg_b = args(1);
  41. +          nr_b = arg_b.rows ();
  42. +          nc_b = arg_b.columns ();
  43.  
  44. -      arg_is_empty = empty_arg ("eig", nr_b, nc_b);
  45. -      if (arg_is_empty < 0)
  46. -        return retval;
  47. -      else if (arg_is_empty > 0)
  48. -        return ovl (2, Matrix ());
  49. +          arg_is_empty = empty_arg ("eig", nr_b, nc_b);
  50. +          if (arg_is_empty < 0)
  51. +            return retval;
  52. +          else if (arg_is_empty > 0)
  53. +            return ovl (2, Matrix ());
  54.  
  55. -      if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
  56. -        err_wrong_type_arg ("eig", arg_b);
  57. +          if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
  58. +            err_wrong_type_arg ("eig", arg_b);
  59. +        }
  60.      }
  61.  
  62.    if (nr_a != nc_a)
  63.      err_square_matrix_required ("eig", "A");
  64. -
  65. -  if (nargin == 2 && nr_b != nc_b)
  66. +  
  67. +  // determine if it's AEP or GEP
  68. +  bool AEPcase = nargin == 1 || args(1).is_string ();
  69. +  
  70. +  if (!AEPcase && nr_b != nc_b)
  71.      err_square_matrix_required ("eig", "B");
  72.  
  73.    Matrix tmp_a, tmp_b;
  74. @@ -106,22 +120,22 @@ The eigenvalues returned by @code{eig} a
  75.      {
  76.        FloatEIG result;
  77.  
  78. -      if (nargin == 1)
  79. +      if (AEPcase)
  80.          {
  81.            if (arg_a.is_real_type ())
  82.              {
  83.                ftmp_a = arg_a.float_matrix_value ();
  84.  
  85. -              result = FloatEIG (ftmp_a, nargout > 1);
  86. +              result = FloatEIG (ftmp_a, nargout > 1, balance);
  87.              }
  88.            else
  89.              {
  90.                fctmp_a = arg_a.float_complex_matrix_value ();
  91.  
  92. -              result = FloatEIG (fctmp_a, nargout > 1);
  93. +              result = FloatEIG (fctmp_a, nargout > 1, balance);
  94.              }
  95.          }
  96. -      else if (nargin == 2)
  97. +      else
  98.          {
  99.            if (arg_a.is_real_type () && arg_b.is_real_type ())
  100.              {
  101. @@ -155,22 +169,22 @@ The eigenvalues returned by @code{eig} a
  102.      {
  103.        EIG result;
  104.  
  105. -      if (nargin == 1)
  106. +      if (AEPcase)
  107.          {
  108.            if (arg_a.is_real_type ())
  109.              {
  110.                tmp_a = arg_a.matrix_value ();
  111.  
  112. -              result = EIG (tmp_a, nargout > 1);
  113. +              result = EIG (tmp_a, nargout > 1, balance);
  114.              }
  115.            else
  116.              {
  117.                ctmp_a = arg_a.complex_matrix_value ();
  118.  
  119. -              result = EIG (ctmp_a, nargout > 1);
  120. +              result = EIG (ctmp_a, nargout > 1, balance);
  121.              }
  122.          }
  123. -      else if (nargin == 2)
  124. +      else
  125.          {
  126.            if (arg_a.is_real_type () && arg_b.is_real_type ())
  127.              {
  128. diff -r 102b33b53ea4 -r c6383dae0bd5 liboctave/numeric/EIG.cc
  129. --- a/liboctave/numeric/EIG.cc  Mon Mar 28 20:00:37 2016 +1100
  130. +++ b/liboctave/numeric/EIG.cc  Tue Mar 29 19:20:54 2016 +0200
  131. @@ -32,6 +32,24 @@ along with Octave; see the file COPYING.
  132.  extern "C"
  133.  {
  134.    F77_RET_T
  135. +  F77_FUNC (dgeevx, DGEEVX) (F77_CONST_CHAR_ARG_DECL,
  136. +                             F77_CONST_CHAR_ARG_DECL,
  137. +                             F77_CONST_CHAR_ARG_DECL,
  138. +                             F77_CONST_CHAR_ARG_DECL,
  139. +                             const octave_idx_type&, double*,
  140. +                             const octave_idx_type&, double*, double*,
  141. +                             double*, const octave_idx_type&, double*,
  142. +                             const octave_idx_type&, octave_idx_type&,
  143. +                             octave_idx_type&, double*, double&,
  144. +                             double*, double*, double*,
  145. +                             const octave_idx_type&, octave_idx_type*,
  146. +                             octave_idx_type&
  147. +                             F77_CHAR_ARG_LEN_DECL
  148. +                             F77_CHAR_ARG_LEN_DECL
  149. +                             F77_CHAR_ARG_LEN_DECL
  150. +                             F77_CHAR_ARG_LEN_DECL);
  151. +  
  152. +  F77_RET_T
  153.    F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL,
  154.                             F77_CONST_CHAR_ARG_DECL,
  155.                             const octave_idx_type&, double*,
  156. @@ -43,6 +61,24 @@ extern "C"
  157.                             F77_CHAR_ARG_LEN_DECL);
  158.  
  159.    F77_RET_T
  160. +  F77_FUNC (zgeevx, ZGEEVX) (F77_CONST_CHAR_ARG_DECL,
  161. +                             F77_CONST_CHAR_ARG_DECL,
  162. +                             F77_CONST_CHAR_ARG_DECL,
  163. +                             F77_CONST_CHAR_ARG_DECL,
  164. +                             const octave_idx_type&, Complex*,
  165. +                             const octave_idx_type&, Complex*,
  166. +                             Complex*, const octave_idx_type&, Complex*,
  167. +                             const octave_idx_type&, octave_idx_type&,
  168. +                             octave_idx_type&, double*, double&,
  169. +                             double*, double*, Complex*,
  170. +                             const octave_idx_type&, double*,
  171. +                             octave_idx_type&
  172. +                             F77_CHAR_ARG_LEN_DECL
  173. +                             F77_CHAR_ARG_LEN_DECL
  174. +                             F77_CHAR_ARG_LEN_DECL
  175. +                             F77_CHAR_ARG_LEN_DECL);
  176. +
  177. +  F77_RET_T
  178.    F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL,
  179.                             F77_CONST_CHAR_ARG_DECL,
  180.                             const octave_idx_type&, Complex*,
  181. @@ -137,7 +173,7 @@ extern "C"
  182.  }
  183.  
  184.  octave_idx_type
  185. -EIG::init (const Matrix& a, bool calc_ev)
  186. +EIG::init (const Matrix& a, bool calc_ev, bool balance)
  187.  {
  188.    if (a.any_element_is_inf_or_nan ())
  189.      (*current_liboctave_error_handler)
  190. @@ -171,33 +207,61 @@ EIG::init (const Matrix& a, bool calc_ev
  191.  
  192.    double *dummy = 0;
  193.    octave_idx_type idummy = 1;
  194. -
  195. -  F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  196. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  197. -                           n, tmp_data, n, pwr, pwi, dummy,
  198. -                           idummy, pvr, n, &dummy_work, lwork, info
  199. -                           F77_CHAR_ARG_LEN (1)
  200. -                           F77_CHAR_ARG_LEN (1)));
  201. +  
  202. +  octave_idx_type ilo;
  203. +  octave_idx_type ihi;
  204. +  
  205. +  Array<double> scale (dim_vector (n, 1));
  206. +  double *pscale = scale.fortran_vec ();
  207. +  
  208. +  double abnrm;
  209. +  
  210. +  Array<double> rconde (dim_vector (n, 1));
  211. +  double *prconde = rconde.fortran_vec ();
  212. +  
  213. +  Array<double> rcondv (dim_vector (n, 1));
  214. +  double *prcondv = rcondv.fortran_vec ();
  215. +  
  216. +  octave_idx_type dummy_iwork;
  217. +        
  218. +  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  219. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  220. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  221. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  222. +                             n, tmp_data, n, pwr, pwi, dummy,
  223. +                             idummy, pvr, n, ilo, ihi, pscale,
  224. +                             abnrm, prconde, prcondv, &dummy_work,
  225. +                             lwork, &dummy_iwork, info
  226. +                             F77_CHAR_ARG_LEN (1)
  227. +                             F77_CHAR_ARG_LEN (1)
  228. +                             F77_CHAR_ARG_LEN (1)
  229. +                             F77_CHAR_ARG_LEN (1)));
  230.  
  231.    if (info != 0)
  232. -    (*current_liboctave_error_handler) ("dgeev workspace query failed");
  233. +    (*current_liboctave_error_handler) ("dgeevx workspace query failed");
  234.  
  235.    lwork = static_cast<octave_idx_type> (dummy_work);
  236.    Array<double> work (dim_vector (lwork, 1));
  237.    double *pwork = work.fortran_vec ();
  238. -
  239. -  F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  240. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  241. -                           n, tmp_data, n, pwr, pwi, dummy,
  242. -                           idummy, pvr, n, pwork, lwork, info
  243. -                           F77_CHAR_ARG_LEN (1)
  244. -                           F77_CHAR_ARG_LEN (1)));
  245. +  
  246. +  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  247. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  248. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  249. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  250. +                             n, tmp_data, n, pwr, pwi, dummy,
  251. +                             idummy, pvr, n, ilo, ihi, pscale,
  252. +                             abnrm, prconde, prcondv, pwork,
  253. +                             lwork, &dummy_iwork, info
  254. +                             F77_CHAR_ARG_LEN (1)
  255. +                             F77_CHAR_ARG_LEN (1)
  256. +                             F77_CHAR_ARG_LEN (1)
  257. +                             F77_CHAR_ARG_LEN (1)));
  258.  
  259.    if (info < 0)
  260. -    (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
  261. +    (*current_liboctave_error_handler) ("unrecoverable error in dgeevx");
  262.  
  263.    if (info > 0)
  264. -    (*current_liboctave_error_handler) ("dgeev failed to converge");
  265. +    (*current_liboctave_error_handler) ("dgeevx failed to converge");
  266.  
  267.    lambda.resize (n);
  268.    octave_idx_type nvr = calc_ev ? n : 0;
  269. @@ -284,7 +348,7 @@ EIG::symmetric_init (const Matrix& a, bo
  270.  }
  271.  
  272.  octave_idx_type
  273. -EIG::init (const ComplexMatrix& a, bool calc_ev)
  274. +EIG::init (const ComplexMatrix& a, bool calc_ev, bool balance)
  275.  {
  276.    if (a.any_element_is_inf_or_nan ())
  277.      (*current_liboctave_error_handler)
  278. @@ -320,32 +384,59 @@ EIG::init (const ComplexMatrix& a, bool
  279.    Complex *dummy = 0;
  280.    octave_idx_type idummy = 1;
  281.  
  282. -  F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  283. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  284. -                           n, tmp_data, n, pw, dummy, idummy,
  285. -                           pv, n, &dummy_work, lwork, prwork, info
  286. -                           F77_CHAR_ARG_LEN (1)
  287. -                           F77_CHAR_ARG_LEN (1)));
  288. +  octave_idx_type ilo;
  289. +  octave_idx_type ihi;
  290. +  
  291. +  Array<double> scale (dim_vector (n, 1));
  292. +  double *pscale = scale.fortran_vec ();
  293. +  
  294. +  double abnrm;
  295. +  
  296. +  Array<double> rconde (dim_vector (n, 1));
  297. +  double *prconde = rconde.fortran_vec ();
  298. +  
  299. +  Array<double> rcondv (dim_vector (n, 1));
  300. +  double *prcondv = rcondv.fortran_vec ();
  301. +  
  302. +  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  303. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  304. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  305. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  306. +                             n, tmp_data, n, pw, dummy, idummy,
  307. +                             pv, n, ilo, ihi, pscale, abnrm,
  308. +                             prconde, prcondv,
  309. +                             &dummy_work, lwork, prwork, info
  310. +                             F77_CHAR_ARG_LEN (1)
  311. +                             F77_CHAR_ARG_LEN (1)
  312. +                             F77_CHAR_ARG_LEN (1)
  313. +                             F77_CHAR_ARG_LEN (1)));
  314.  
  315.    if (info != 0)
  316. -    (*current_liboctave_error_handler) ("zgeev workspace query failed");
  317. +    (*current_liboctave_error_handler) ("zgeevx workspace query failed");
  318.  
  319.    lwork = static_cast<octave_idx_type> (dummy_work.real ());
  320.    Array<Complex> work (dim_vector (lwork, 1));
  321.    Complex *pwork = work.fortran_vec ();
  322. -
  323. -  F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  324. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  325. -                           n, tmp_data, n, pw, dummy, idummy,
  326. -                           pv, n, pwork, lwork, prwork, info
  327. -                           F77_CHAR_ARG_LEN (1)
  328. -                           F77_CHAR_ARG_LEN (1)));
  329. +  
  330. +  
  331. +  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  332. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  333. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  334. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  335. +                             n, tmp_data, n, pw, dummy, idummy,
  336. +                             pv, n, ilo, ihi, pscale, abnrm,
  337. +                             prconde, prcondv,
  338. +                             pwork, lwork, prwork, info
  339. +                             F77_CHAR_ARG_LEN (1)
  340. +                             F77_CHAR_ARG_LEN (1)
  341. +                             F77_CHAR_ARG_LEN (1)
  342. +                             F77_CHAR_ARG_LEN (1)));
  343.  
  344.    if (info < 0)
  345. -    (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
  346. +    (*current_liboctave_error_handler) ("unrecoverable error in zgeevx");
  347.  
  348.    if (info > 0)
  349. -    (*current_liboctave_error_handler) ("zgeev failed to converge");
  350. +    (*current_liboctave_error_handler) ("zgeevx failed to converge");
  351.  
  352.    lambda = w;
  353.    v = vtmp;
  354. diff -r 102b33b53ea4 -r c6383dae0bd5 liboctave/numeric/EIG.h
  355. --- a/liboctave/numeric/EIG.h   Mon Mar 28 20:00:37 2016 +1100
  356. +++ b/liboctave/numeric/EIG.h   Tue Mar 29 19:20:54 2016 +0200
  357. @@ -42,16 +42,18 @@ public:
  358.  
  359.    EIG (void) : lambda (), v () { }
  360.  
  361. -  EIG (const Matrix& a, bool calc_eigenvectors = true)
  362. +  EIG (const Matrix& a, bool calc_eigenvectors = true,
  363. +       bool balance = true)
  364.      : lambda (), v ()
  365.    {
  366. -    init (a, calc_eigenvectors);
  367. +    init (a, calc_eigenvectors, balance);
  368.    }
  369.  
  370. -  EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
  371. +  EIG (const Matrix& a, octave_idx_type& info,
  372. +       bool calc_eigenvectors = true, bool balance = true)
  373.      : lambda (), v ()
  374.    {
  375. -    info = init (a, calc_eigenvectors);
  376. +    info = init (a, calc_eigenvectors, balance);
  377.    }
  378.  
  379.    EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true)
  380. @@ -67,17 +69,18 @@ public:
  381.      info = init (a, b, calc_eigenvectors);
  382.    }
  383.  
  384. -  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
  385. +  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true,
  386. +       bool balance = true)
  387.      : lambda (), v ()
  388.    {
  389. -    init (a, calc_eigenvectors);
  390. +    init (a, calc_eigenvectors, balance);
  391.    }
  392.  
  393.    EIG (const ComplexMatrix& a, octave_idx_type& info,
  394. -       bool calc_eigenvectors = true)
  395. +       bool calc_eigenvectors = true, bool balance = true)
  396.      : lambda (), v ()
  397.    {
  398. -    info = init (a, calc_eigenvectors);
  399. +    info = init (a, calc_eigenvectors, balance);
  400.    }
  401.  
  402.    EIG (const ComplexMatrix& a, const ComplexMatrix& b,
  403. @@ -120,12 +123,14 @@ private:
  404.    ComplexColumnVector lambda;
  405.    ComplexMatrix v;
  406.  
  407. -  octave_idx_type init (const Matrix& a, bool calc_eigenvectors);
  408. +  octave_idx_type init (const Matrix& a, bool calc_eigenvectors,
  409. +                        bool balance);
  410.  
  411.    octave_idx_type init (const Matrix& a, const Matrix& b,
  412.                          bool calc_eigenvectors);
  413.  
  414. -  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors);
  415. +  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors,
  416. +                        bool balance);
  417.  
  418.    octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b,
  419.                          bool calc_eigenvectors);
  420. diff -r 102b33b53ea4 -r c6383dae0bd5 liboctave/numeric/fEIG.cc
  421. --- a/liboctave/numeric/fEIG.cc Mon Mar 28 20:00:37 2016 +1100
  422. +++ b/liboctave/numeric/fEIG.cc Tue Mar 29 19:20:54 2016 +0200
  423. @@ -32,6 +32,23 @@ along with Octave; see the file COPYING.
  424.  extern "C"
  425.  {
  426.    F77_RET_T
  427. +  F77_FUNC (sgeevx, SGEEVX) (F77_CONST_CHAR_ARG_DECL,
  428. +                             F77_CONST_CHAR_ARG_DECL,
  429. +                             F77_CONST_CHAR_ARG_DECL,
  430. +                             F77_CONST_CHAR_ARG_DECL,
  431. +                             const octave_idx_type&, float*,
  432. +                             const octave_idx_type&, float*, float*, float*,
  433. +                             const octave_idx_type&, float*,
  434. +                             const octave_idx_type&, octave_idx_type&,
  435. +                             octave_idx_type&, float*, float&, float*,
  436. +                             float*, float*, const octave_idx_type&,
  437. +                             octave_idx_type*, octave_idx_type&
  438. +                             F77_CHAR_ARG_LEN_DECL
  439. +                             F77_CHAR_ARG_LEN_DECL
  440. +                             F77_CHAR_ARG_LEN_DECL
  441. +                             F77_CHAR_ARG_LEN_DECL);
  442. +
  443. +  F77_RET_T
  444.    F77_FUNC (sgeev, SGEEV) (F77_CONST_CHAR_ARG_DECL,
  445.                             F77_CONST_CHAR_ARG_DECL,
  446.                             const octave_idx_type&, float*,
  447. @@ -41,6 +58,23 @@ extern "C"
  448.                             const octave_idx_type&, octave_idx_type&
  449.                             F77_CHAR_ARG_LEN_DECL
  450.                             F77_CHAR_ARG_LEN_DECL);
  451. +                          
  452. +  F77_RET_T
  453. +  F77_FUNC (cgeevx, CGEEVX) (F77_CONST_CHAR_ARG_DECL,
  454. +                             F77_CONST_CHAR_ARG_DECL,
  455. +                             F77_CONST_CHAR_ARG_DECL,
  456. +                             F77_CONST_CHAR_ARG_DECL,
  457. +                             const octave_idx_type&, FloatComplex*,
  458. +                             const octave_idx_type&, FloatComplex*, FloatComplex*,
  459. +                             const octave_idx_type&, FloatComplex*,
  460. +                             const octave_idx_type&, octave_idx_type&,
  461. +                             octave_idx_type&, float*, float&, float*,
  462. +                             float*, FloatComplex*, const octave_idx_type&,
  463. +                             float*, octave_idx_type&
  464. +                             F77_CHAR_ARG_LEN_DECL
  465. +                             F77_CHAR_ARG_LEN_DECL
  466. +                             F77_CHAR_ARG_LEN_DECL
  467. +                             F77_CHAR_ARG_LEN_DECL);
  468.  
  469.    F77_RET_T
  470.    F77_FUNC (cgeev, CGEEV) (F77_CONST_CHAR_ARG_DECL,
  471. @@ -134,7 +168,7 @@ extern "C"
  472.  }
  473.  
  474.  octave_idx_type
  475. -FloatEIG::init (const FloatMatrix& a, bool calc_ev)
  476. +FloatEIG::init (const FloatMatrix& a, bool calc_ev, bool balance)
  477.  {
  478.    if (a.any_element_is_inf_or_nan ())
  479.      (*current_liboctave_error_handler)
  480. @@ -169,32 +203,60 @@ FloatEIG::init (const FloatMatrix& a, bo
  481.    float *dummy = 0;
  482.    octave_idx_type idummy = 1;
  483.  
  484. -  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  485. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  486. -                           n, tmp_data, n, pwr, pwi, dummy,
  487. -                           idummy, pvr, n, &dummy_work, lwork, info
  488. -                           F77_CHAR_ARG_LEN (1)
  489. -                           F77_CHAR_ARG_LEN (1)));
  490. +  octave_idx_type ilo;
  491. +  octave_idx_type ihi;
  492. +  
  493. +  Array<float> scale (dim_vector (n, 1));
  494. +  float *pscale = scale.fortran_vec ();
  495. +  
  496. +  float abnrm;
  497. +  
  498. +  Array<float> rconde (dim_vector (n, 1));
  499. +  float *prconde = rconde.fortran_vec ();
  500. +  
  501. +  Array<float> rcondv (dim_vector (n, 1));
  502. +  float *prcondv = rcondv.fortran_vec ();
  503. +  
  504. +  octave_idx_type dummy_iwork;
  505. +  
  506. +  F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  507. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  508. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  509. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  510. +                             n, tmp_data, n, pwr, pwi, dummy,
  511. +                             idummy, pvr, n,
  512. +                             ilo, ihi, pscale, abnrm, prconde, prcondv,
  513. +                             &dummy_work, lwork, &dummy_iwork, info
  514. +                             F77_CHAR_ARG_LEN (1)
  515. +                             F77_CHAR_ARG_LEN (1)
  516. +                             F77_CHAR_ARG_LEN (1)
  517. +                             F77_CHAR_ARG_LEN (1)));                        
  518.  
  519.    if (info != 0)
  520. -    (*current_liboctave_error_handler) ("sgeev workspace query failed");
  521. +    (*current_liboctave_error_handler) ("sgeevx workspace query failed");
  522.  
  523.    lwork = static_cast<octave_idx_type> (dummy_work);
  524.    Array<float> work (dim_vector (lwork, 1));
  525.    float *pwork = work.fortran_vec ();
  526. -
  527. -  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  528. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  529. -                           n, tmp_data, n, pwr, pwi, dummy,
  530. -                           idummy, pvr, n, pwork, lwork, info
  531. -                           F77_CHAR_ARG_LEN (1)
  532. -                           F77_CHAR_ARG_LEN (1)));
  533. +  
  534. +  F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  535. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  536. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  537. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  538. +                             n, tmp_data, n, pwr, pwi, dummy,
  539. +                             idummy, pvr, n,
  540. +                             ilo, ihi, pscale, abnrm, prconde, prcondv,
  541. +                             pwork, lwork, &dummy_iwork, info
  542. +                             F77_CHAR_ARG_LEN (1)
  543. +                             F77_CHAR_ARG_LEN (1)
  544. +                             F77_CHAR_ARG_LEN (1)
  545. +                             F77_CHAR_ARG_LEN (1)));
  546.  
  547.    if (info < 0)
  548. -    (*current_liboctave_error_handler) ("unrecoverable error in sgeev");
  549. +    (*current_liboctave_error_handler) ("unrecoverable error in sgeevx");
  550.  
  551.    if (info > 0)
  552. -    (*current_liboctave_error_handler) ("sgeev failed to converge");
  553. +    (*current_liboctave_error_handler) ("sgeevx failed to converge");
  554.  
  555.    lambda.resize (n);
  556.    v.resize (nvr, nvr);
  557. @@ -280,7 +342,7 @@ FloatEIG::symmetric_init (const FloatMat
  558.  }
  559.  
  560.  octave_idx_type
  561. -FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev)
  562. +FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev, bool balance)
  563.  {
  564.    if (a.any_element_is_inf_or_nan ())
  565.      (*current_liboctave_error_handler)
  566. @@ -316,32 +378,56 @@ FloatEIG::init (const FloatComplexMatrix
  567.    FloatComplex *dummy = 0;
  568.    octave_idx_type idummy = 1;
  569.  
  570. -  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  571. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  572. -                           n, tmp_data, n, pw, dummy, idummy,
  573. -                           pv, n, &dummy_work, lwork, prwork, info
  574. -                           F77_CHAR_ARG_LEN (1)
  575. -                           F77_CHAR_ARG_LEN (1)));
  576. +  octave_idx_type ilo;
  577. +  octave_idx_type ihi;
  578. +  
  579. +  Array<float> scale (dim_vector (n, 1));
  580. +  float *pscale = scale.fortran_vec ();
  581. +  
  582. +  float abnrm;
  583. +  
  584. +  Array<float> rconde (dim_vector (n, 1));
  585. +  float *prconde = rconde.fortran_vec ();
  586. +  
  587. +  Array<float> rcondv (dim_vector (n, 1));
  588. +  float *prcondv = rcondv.fortran_vec ();
  589. +  
  590. +  F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  591. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  592. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  593. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  594. +                             n, tmp_data, n, pw, dummy, idummy,
  595. +                             pv, n, ilo, ihi, pscale, abnrm, prconde, prcondv,
  596. +                             &dummy_work, lwork, prwork, info
  597. +                             F77_CHAR_ARG_LEN (1)
  598. +                             F77_CHAR_ARG_LEN (1)
  599. +                             F77_CHAR_ARG_LEN (1)
  600. +                             F77_CHAR_ARG_LEN (1)));
  601.  
  602.    if (info != 0)
  603. -    (*current_liboctave_error_handler) ("cgeev workspace query failed");
  604. +    (*current_liboctave_error_handler) ("cgeevx workspace query failed");
  605.  
  606.    lwork = static_cast<octave_idx_type> (dummy_work.real ());
  607.    Array<FloatComplex> work (dim_vector (lwork, 1));
  608.    FloatComplex *pwork = work.fortran_vec ();
  609.  
  610. -  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  611. -                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  612. -                           n, tmp_data, n, pw, dummy, idummy,
  613. -                           pv, n, pwork, lwork, prwork, info
  614. -                           F77_CHAR_ARG_LEN (1)
  615. -                           F77_CHAR_ARG_LEN (1)));
  616. +  F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  617. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  618. +                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  619. +                             F77_CONST_CHAR_ARG2 ("N", 1),
  620. +                             n, tmp_data, n, pw, dummy, idummy,
  621. +                             pv, n, ilo, ihi, pscale,abnrm, prconde, prcondv,
  622. +                             pwork, lwork, prwork, info
  623. +                             F77_CHAR_ARG_LEN (1)
  624. +                             F77_CHAR_ARG_LEN (1)
  625. +                             F77_CHAR_ARG_LEN (1)
  626. +                             F77_CHAR_ARG_LEN (1)));  
  627.  
  628.    if (info < 0)
  629. -    (*current_liboctave_error_handler) ("unrecoverable error in cgeev");
  630. +    (*current_liboctave_error_handler) ("unrecoverable error in cgeevx");
  631.  
  632.    if (info > 0)
  633. -    (*current_liboctave_error_handler) ("cgeev failed to converge");
  634. +    (*current_liboctave_error_handler) ("cgeevx failed to converge");
  635.  
  636.    lambda = w;
  637.    v = vtmp;
  638. diff -r 102b33b53ea4 -r c6383dae0bd5 liboctave/numeric/fEIG.h
  639. --- a/liboctave/numeric/fEIG.h  Mon Mar 28 20:00:37 2016 +1100
  640. +++ b/liboctave/numeric/fEIG.h  Tue Mar 29 19:20:54 2016 +0200
  641. @@ -43,17 +43,18 @@ public:
  642.    FloatEIG (void)
  643.      : lambda (), v () { }
  644.  
  645. -  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true)
  646. +  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true,
  647. +            bool balance = true)
  648.      : lambda (), v ()
  649.    {
  650. -    init (a, calc_eigenvectors);
  651. +    init (a, calc_eigenvectors, balance);
  652.    }
  653.  
  654.    FloatEIG (const FloatMatrix& a, octave_idx_type& info,
  655. -            bool calc_eigenvectors = true)
  656. +            bool calc_eigenvectors = true, bool balance = true)
  657.      : lambda (), v ()
  658.    {
  659. -    info = init (a, calc_eigenvectors);
  660. +    info = init (a, calc_eigenvectors, balance);
  661.    }
  662.  
  663.    FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
  664. @@ -63,24 +64,25 @@ public:
  665.      init (a, b, calc_eigenvectors);
  666.    }
  667.  
  668. -  FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info,
  669. -            bool calc_eigenvectors = true)
  670. +  FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
  671. +            octave_idx_type& info, bool calc_eigenvectors = true)
  672.      : lambda (), v ()
  673.    {
  674.      info = init (a, b, calc_eigenvectors);
  675.    }
  676.  
  677. -  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true)
  678. +  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true,
  679. +            bool balance = true)
  680.      : lambda (), v ()
  681.    {
  682. -    init (a, calc_eigenvectors);
  683. +    init (a, calc_eigenvectors, balance);
  684.    }
  685.  
  686.    FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info,
  687. -            bool calc_eigenvectors = true)
  688. +            bool calc_eigenvectors = true, bool balance = true)
  689.      : lambda (), v ()
  690.    {
  691. -    info = init (a, calc_eigenvectors);
  692. +    info = init (a, calc_eigenvectors, balance);
  693.    }
  694.  
  695.    FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
  696. @@ -122,10 +124,12 @@ private:
  697.    FloatComplexColumnVector lambda;
  698.    FloatComplexMatrix v;
  699.  
  700. -  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors);
  701. +  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors,
  702. +                        bool balance);
  703.    octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b,
  704.                          bool calc_eigenvectors);
  705. -  octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors);
  706. +  octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors,
  707. +                        bool balance);
  708.    octave_idx_type init (const FloatComplexMatrix& a,
  709.                          const FloatComplexMatrix& b, bool calc_eigenvectors);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!