API tools faq
paste
Advertisement
BumiBarbi

EIG.cc

BumiBarbi
Mar 29th, 2016
120
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 29.03 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  
  3. Copyright (C) 1994-2015 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 3 of the License, or (at your
  10. option) any later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, see
  19. <http://www.gnu.org/licenses/>.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #  include "config.h"
  25. #endif
  26.  
  27. #include "EIG.h"
  28. #include "dColVector.h"
  29. #include "f77-fcn.h"
  30. #include "lo-error.h"
  31.  
  32. extern "C"
  33. {
  34.   F77_RET_T
  35.   F77_FUNC (dgeevx, DGEEVX) (F77_CONST_CHAR_ARG_DECL,
  36.                              F77_CONST_CHAR_ARG_DECL,
  37.                              F77_CONST_CHAR_ARG_DECL,
  38.                              F77_CONST_CHAR_ARG_DECL,
  39.                              const octave_idx_type&, double*,
  40.                              const octave_idx_type&, double*, double*,
  41.                              double*, const octave_idx_type&, double*,
  42.                              const octave_idx_type&, octave_idx_type&,
  43.                              octave_idx_type&, double*, double&,
  44.                              double*, double*, double*,
  45.                              const octave_idx_type&, octave_idx_type*,
  46.                              octave_idx_type&
  47.                              F77_CHAR_ARG_LEN_DECL
  48.                              F77_CHAR_ARG_LEN_DECL
  49.                              F77_CHAR_ARG_LEN_DECL
  50.                              F77_CHAR_ARG_LEN_DECL);
  51.  
  52.   F77_RET_T
  53.   F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL,
  54.                            F77_CONST_CHAR_ARG_DECL,
  55.                            const octave_idx_type&, double*,
  56.                            const octave_idx_type&, double*, double*,
  57.                            double*, const octave_idx_type&, double*,
  58.                            const octave_idx_type&, double*,
  59.                            const octave_idx_type&, octave_idx_type&
  60.                            F77_CHAR_ARG_LEN_DECL
  61.                            F77_CHAR_ARG_LEN_DECL);
  62.  
  63.   F77_RET_T
  64.   F77_FUNC (zgeevx, ZGEEVX) (F77_CONST_CHAR_ARG_DECL,
  65.                              F77_CONST_CHAR_ARG_DECL,
  66.                              F77_CONST_CHAR_ARG_DECL,
  67.                              F77_CONST_CHAR_ARG_DECL,
  68.                              const octave_idx_type&, Complex*,
  69.                              const octave_idx_type&, Complex*,
  70.                              Complex*, const octave_idx_type&, Complex*,
  71.                              const octave_idx_type&, octave_idx_type&,
  72.                              octave_idx_type&, double*, double&,
  73.                              double*, double*, Complex*,
  74.                              const octave_idx_type&, double*,
  75.                              octave_idx_type&
  76.                              F77_CHAR_ARG_LEN_DECL
  77.                              F77_CHAR_ARG_LEN_DECL
  78.                              F77_CHAR_ARG_LEN_DECL
  79.                              F77_CHAR_ARG_LEN_DECL);
  80.  
  81.   F77_RET_T
  82.   F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL,
  83.                            F77_CONST_CHAR_ARG_DECL,
  84.                            const octave_idx_type&, Complex*,
  85.                            const octave_idx_type&, Complex*,
  86.                            Complex*, const octave_idx_type&, Complex*,
  87.                            const octave_idx_type&, Complex*,
  88.                            const octave_idx_type&, double*, octave_idx_type&
  89.                            F77_CHAR_ARG_LEN_DECL
  90.                            F77_CHAR_ARG_LEN_DECL);
  91.  
  92.   F77_RET_T
  93.   F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL,
  94.                            F77_CONST_CHAR_ARG_DECL,
  95.                            const octave_idx_type&, double*,
  96.                            const octave_idx_type&, double*, double*,
  97.                            const octave_idx_type&, octave_idx_type&
  98.                            F77_CHAR_ARG_LEN_DECL
  99.                            F77_CHAR_ARG_LEN_DECL);
  100.  
  101.   F77_RET_T
  102.   F77_FUNC (zheev, ZHEEV) (F77_CONST_CHAR_ARG_DECL,
  103.                            F77_CONST_CHAR_ARG_DECL,
  104.                            const octave_idx_type&, Complex*,
  105.                            const octave_idx_type&, double*,
  106.                            Complex*, const octave_idx_type&, double*,
  107.                            octave_idx_type&
  108.                            F77_CHAR_ARG_LEN_DECL
  109.                            F77_CHAR_ARG_LEN_DECL);
  110.  
  111.   F77_RET_T
  112.   F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL,
  113.                              const octave_idx_type&, double*,
  114.                              const octave_idx_type&, octave_idx_type&
  115.                              F77_CHAR_ARG_LEN_DECL
  116.                              F77_CHAR_ARG_LEN_DECL);
  117.  
  118.   F77_RET_T
  119.   F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL,
  120.                              const octave_idx_type&,
  121.                              Complex*, const octave_idx_type&,
  122.                              octave_idx_type&
  123.                              F77_CHAR_ARG_LEN_DECL
  124.                              F77_CHAR_ARG_LEN_DECL);
  125.  
  126.   F77_RET_T
  127.   F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL,
  128.                            F77_CONST_CHAR_ARG_DECL,
  129.                            const octave_idx_type&,
  130.                            double*, const octave_idx_type&,
  131.                            double*, const octave_idx_type&,
  132.                            double*, double*, double *, double*,
  133.                            const octave_idx_type&, double*,
  134.                            const octave_idx_type&, double*,
  135.                            const octave_idx_type&, octave_idx_type&
  136.                            F77_CHAR_ARG_LEN_DECL
  137.                            F77_CHAR_ARG_LEN_DECL);
  138.  
  139.   F77_RET_T
  140.   F77_FUNC (dsygv, DSYGV) (const octave_idx_type&,
  141.                            F77_CONST_CHAR_ARG_DECL,
  142.                            F77_CONST_CHAR_ARG_DECL,
  143.                            const octave_idx_type&, double*,
  144.                            const octave_idx_type&, double*,
  145.                            const octave_idx_type&, double*, double*,
  146.                            const octave_idx_type&, octave_idx_type&
  147.                            F77_CHAR_ARG_LEN_DECL
  148.                            F77_CHAR_ARG_LEN_DECL);
  149.  
  150.   F77_RET_T
  151.   F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL,
  152.                            F77_CONST_CHAR_ARG_DECL,
  153.                            const octave_idx_type&,
  154.                            Complex*, const octave_idx_type&,
  155.                            Complex*, const octave_idx_type&,
  156.                            Complex*, Complex*, Complex*,
  157.                            const octave_idx_type&, Complex*,
  158.                            const octave_idx_type&, Complex*,
  159.                            const octave_idx_type&, double*, octave_idx_type&
  160.                            F77_CHAR_ARG_LEN_DECL
  161.                            F77_CHAR_ARG_LEN_DECL);
  162.  
  163.   F77_RET_T
  164.   F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&,
  165.                            F77_CONST_CHAR_ARG_DECL,
  166.                            F77_CONST_CHAR_ARG_DECL,
  167.                            const octave_idx_type&, Complex*,
  168.                            const octave_idx_type&, Complex*,
  169.                            const octave_idx_type&, double*, Complex*,
  170.                            const octave_idx_type&, double*, octave_idx_type&
  171.                            F77_CHAR_ARG_LEN_DECL
  172.                            F77_CHAR_ARG_LEN_DECL);
  173. }
  174.  
  175. octave_idx_type
  176. EIG::init (const Matrix& a, bool calc_ev, bool balance)
  177. {
  178.   if (a.any_element_is_inf_or_nan ())
  179.     (*current_liboctave_error_handler)
  180.       ("EIG: matrix contains Inf or NaN values");
  181.  
  182.   if (a.is_symmetric ())
  183.     return symmetric_init (a, calc_ev);
  184.  
  185.   octave_idx_type n = a.rows ();
  186.  
  187.   if (n != a.cols ())
  188.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  189.  
  190.   octave_idx_type info = 0;
  191.  
  192.   Matrix atmp = a;
  193.   double *tmp_data = atmp.fortran_vec ();
  194.  
  195.   Array<double> wr (dim_vector (n, 1));
  196.   double *pwr = wr.fortran_vec ();
  197.  
  198.   Array<double> wi (dim_vector (n, 1));
  199.   double *pwi = wi.fortran_vec ();
  200.  
  201.   octave_idx_type tnvr = calc_ev ? n : 0;
  202.   Matrix vr (tnvr, tnvr);
  203.   double *pvr = vr.fortran_vec ();
  204.  
  205.   octave_idx_type lwork = -1;
  206.   double dummy_work;
  207.  
  208.   double *dummy = 0;
  209.   octave_idx_type idummy = 1;
  210.  
  211.   octave_idx_type ilo;
  212.   octave_idx_type ihi;
  213.  
  214.   Array<double> scale (dim_vector (n, 1));
  215.   double *pscale = scale.fortran_vec ();
  216.  
  217.   double abnrm;
  218.  
  219.   Array<double> rconde (dim_vector (n, 1));
  220.   double *prconde = rconde.fortran_vec ();
  221.  
  222.   Array<double> rcondv (dim_vector (n, 1));
  223.   double *prcondv = rcondv.fortran_vec ();
  224.  
  225.   octave_idx_type dummy_iwork;
  226.        
  227.   F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  228.                              F77_CONST_CHAR_ARG2 ("N", 1),
  229.                              F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  230.                              F77_CONST_CHAR_ARG2 ("N", 1),
  231.                              n, tmp_data, n, pwr, pwi, dummy,
  232.                              idummy, pvr, n, ilo, ihi, pscale,
  233.                              abnrm, prconde, prcondv, &dummy_work,
  234.                              lwork, &dummy_iwork, info
  235.                              F77_CHAR_ARG_LEN (1)
  236.                              F77_CHAR_ARG_LEN (1)
  237.                              F77_CHAR_ARG_LEN (1)
  238.                              F77_CHAR_ARG_LEN (1)));
  239.  
  240.   if (info != 0)
  241.     (*current_liboctave_error_handler) ("dgeevx workspace query failed");
  242.  
  243.   lwork = static_cast<octave_idx_type> (dummy_work);
  244.   Array<double> work (dim_vector (lwork, 1));
  245.   double *pwork = work.fortran_vec ();
  246.  
  247.   F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  248.                              F77_CONST_CHAR_ARG2 ("N", 1),
  249.                              F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  250.                              F77_CONST_CHAR_ARG2 ("N", 1),
  251.                              n, tmp_data, n, pwr, pwi, dummy,
  252.                              idummy, pvr, n, ilo, ihi, pscale,
  253.                              abnrm, prconde, prcondv, pwork,
  254.                              lwork, &dummy_iwork, info
  255.                              F77_CHAR_ARG_LEN (1)
  256.                              F77_CHAR_ARG_LEN (1)
  257.                              F77_CHAR_ARG_LEN (1)
  258.                              F77_CHAR_ARG_LEN (1)));
  259.  
  260.   if (info < 0)
  261.     (*current_liboctave_error_handler) ("unrecoverable error in dgeevx");
  262.  
  263.   if (info > 0)
  264.     (*current_liboctave_error_handler) ("dgeevx failed to converge");
  265.  
  266.   lambda.resize (n);
  267.   octave_idx_type nvr = calc_ev ? n : 0;
  268.   v.resize (nvr, nvr);
  269.  
  270.   for (octave_idx_type j = 0; j < n; j++)
  271.     {
  272.       if (wi.elem (j) == 0.0)
  273.         {
  274.           lambda.elem (j) = Complex (wr.elem (j));
  275.           for (octave_idx_type i = 0; i < nvr; i++)
  276.             v.elem (i, j) = vr.elem (i, j);
  277.         }
  278.       else
  279.         {
  280.           if (j+1 >= n)
  281.             (*current_liboctave_error_handler) ("EIG: internal error");
  282.  
  283.           lambda.elem (j) = Complex (wr.elem (j), wi.elem (j));
  284.           lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1));
  285.  
  286.           for (octave_idx_type i = 0; i < nvr; i++)
  287.             {
  288.               double real_part = vr.elem (i, j);
  289.               double imag_part = vr.elem (i, j+1);
  290.               v.elem (i, j) = Complex (real_part, imag_part);
  291.               v.elem (i, j+1) = Complex (real_part, -imag_part);
  292.             }
  293.           j++;
  294.         }
  295.     }
  296.  
  297.   return info;
  298. }
  299.  
  300. octave_idx_type
  301. EIG::symmetric_init (const Matrix& a, bool calc_ev)
  302. {
  303.   octave_idx_type n = a.rows ();
  304.  
  305.   if (n != a.cols ())
  306.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  307.  
  308.   octave_idx_type info = 0;
  309.  
  310.   Matrix atmp = a;
  311.   double *tmp_data = atmp.fortran_vec ();
  312.  
  313.   ColumnVector wr (n);
  314.   double *pwr = wr.fortran_vec ();
  315.  
  316.   octave_idx_type lwork = -1;
  317.   double dummy_work;
  318.  
  319.   F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  320.                            F77_CONST_CHAR_ARG2 ("U", 1),
  321.                            n, tmp_data, n, pwr, &dummy_work, lwork, info
  322.                            F77_CHAR_ARG_LEN (1)
  323.                            F77_CHAR_ARG_LEN (1)));
  324.  
  325.   if (info != 0)
  326.     (*current_liboctave_error_handler) ("dsyev workspace query failed");
  327.  
  328.   lwork = static_cast<octave_idx_type> (dummy_work);
  329.   Array<double> work (dim_vector (lwork, 1));
  330.   double *pwork = work.fortran_vec ();
  331.  
  332.   F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  333.                            F77_CONST_CHAR_ARG2 ("U", 1),
  334.                            n, tmp_data, n, pwr, pwork, lwork, info
  335.                            F77_CHAR_ARG_LEN (1)
  336.                            F77_CHAR_ARG_LEN (1)));
  337.  
  338.   if (info < 0)
  339.     (*current_liboctave_error_handler) ("unrecoverable error in dsyev");
  340.  
  341.   if (info > 0)
  342.     (*current_liboctave_error_handler) ("dsyev failed to converge");
  343.  
  344.   lambda = ComplexColumnVector (wr);
  345.   v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
  346.  
  347.   return info;
  348. }
  349.  
  350. octave_idx_type
  351. EIG::init (const ComplexMatrix& a, bool calc_ev, bool balance)
  352. {
  353.   if (a.any_element_is_inf_or_nan ())
  354.     (*current_liboctave_error_handler)
  355.       ("EIG: matrix contains Inf or NaN values");
  356.  
  357.   if (a.is_hermitian ())
  358.     return hermitian_init (a, calc_ev);
  359.  
  360.   octave_idx_type n = a.rows ();
  361.  
  362.   if (n != a.cols ())
  363.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  364.  
  365.   octave_idx_type info = 0;
  366.  
  367.   ComplexMatrix atmp = a;
  368.   Complex *tmp_data = atmp.fortran_vec ();
  369.  
  370.   ComplexColumnVector w (n);
  371.   Complex *pw = w.fortran_vec ();
  372.  
  373.   octave_idx_type nvr = calc_ev ? n : 0;
  374.   ComplexMatrix vtmp (nvr, nvr);
  375.   Complex *pv = vtmp.fortran_vec ();
  376.  
  377.   octave_idx_type lwork = -1;
  378.   Complex dummy_work;
  379.  
  380.   octave_idx_type lrwork = 2*n;
  381.   Array<double> rwork (dim_vector (lrwork, 1));
  382.   double *prwork = rwork.fortran_vec ();
  383.  
  384.   Complex *dummy = 0;
  385.   octave_idx_type idummy = 1;
  386.  
  387.   octave_idx_type ilo;
  388.   octave_idx_type ihi;
  389.  
  390.   Array<double> scale (dim_vector (n, 1));
  391.   double *pscale = scale.fortran_vec ();
  392.  
  393.   double abnrm;
  394.  
  395.   Array<double> rconde (dim_vector (n, 1));
  396.   double *prconde = rconde.fortran_vec ();
  397.  
  398.   Array<double> rcondv (dim_vector (n, 1));
  399.   double *prcondv = rcondv.fortran_vec ();
  400.    
  401.   F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  402.                              F77_CONST_CHAR_ARG2 ("N", 1),
  403.                              F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  404.                              F77_CONST_CHAR_ARG2 ("N", 1),
  405.                              n, tmp_data, n, pw, dummy, idummy,
  406.                              pv, n, ilo, ihi, pscale, abnrm,
  407.                              prconde, prcondv,
  408.                              &dummy_work, lwork, prwork, info
  409.                              F77_CHAR_ARG_LEN (1)
  410.                              F77_CHAR_ARG_LEN (1)
  411.                              F77_CHAR_ARG_LEN (1)
  412.                              F77_CHAR_ARG_LEN (1)));
  413.  
  414.   if (info != 0)
  415.     (*current_liboctave_error_handler) ("zgeevx workspace query failed");
  416.  
  417.   lwork = static_cast<octave_idx_type> (dummy_work.real ());
  418.   Array<Complex> work (dim_vector (lwork, 1));
  419.   Complex *pwork = work.fortran_vec ();
  420.  
  421.   F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
  422.                              F77_CONST_CHAR_ARG2 ("N", 1),
  423.                              F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  424.                              F77_CONST_CHAR_ARG2 ("N", 1),
  425.                              n, tmp_data, n, pw, dummy, idummy,
  426.                              pv, n, ilo, ihi, pscale, abnrm,
  427.                              prconde, prcondv,
  428.                              pwork, lwork, prwork, info
  429.                              F77_CHAR_ARG_LEN (1)
  430.                              F77_CHAR_ARG_LEN (1)
  431.                              F77_CHAR_ARG_LEN (1)
  432.                              F77_CHAR_ARG_LEN (1)));
  433.  
  434.   if (info < 0)
  435.     (*current_liboctave_error_handler) ("unrecoverable error in zgeevx");
  436.  
  437.   if (info > 0)
  438.     (*current_liboctave_error_handler) ("zgeevx failed to converge");
  439.  
  440.   lambda = w;
  441.   v = vtmp;
  442.  
  443.   return info;
  444. }
  445.  
  446. octave_idx_type
  447. EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev)
  448. {
  449.   octave_idx_type n = a.rows ();
  450.  
  451.   if (n != a.cols ())
  452.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  453.  
  454.   octave_idx_type info = 0;
  455.  
  456.   ComplexMatrix atmp = a;
  457.   Complex *tmp_data = atmp.fortran_vec ();
  458.  
  459.   ColumnVector wr (n);
  460.   double *pwr = wr.fortran_vec ();
  461.  
  462.   octave_idx_type lwork = -1;
  463.   Complex dummy_work;
  464.  
  465.   octave_idx_type lrwork = 3*n;
  466.   Array<double> rwork (dim_vector (lrwork, 1));
  467.   double *prwork = rwork.fortran_vec ();
  468.  
  469.   F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  470.                            F77_CONST_CHAR_ARG2 ("U", 1),
  471.                            n, tmp_data, n, pwr, &dummy_work, lwork,
  472.                            prwork, info
  473.                            F77_CHAR_ARG_LEN (1)
  474.                            F77_CHAR_ARG_LEN (1)));
  475.  
  476.   if (info != 0)
  477.     (*current_liboctave_error_handler) ("zheev workspace query failed");
  478.  
  479.   lwork = static_cast<octave_idx_type> (dummy_work.real ());
  480.   Array<Complex> work (dim_vector (lwork, 1));
  481.   Complex *pwork = work.fortran_vec ();
  482.  
  483.   F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  484.                            F77_CONST_CHAR_ARG2 ("U", 1),
  485.                            n, tmp_data, n, pwr, pwork, lwork, prwork, info
  486.                            F77_CHAR_ARG_LEN (1)
  487.                            F77_CHAR_ARG_LEN (1)));
  488.  
  489.   if (info < 0)
  490.     (*current_liboctave_error_handler) ("unrecoverable error in zheev");
  491.  
  492.   if (info > 0)
  493.     (*current_liboctave_error_handler) ("zheev failed to converge");
  494.  
  495.   lambda = ComplexColumnVector (wr);
  496.   v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
  497.  
  498.   return info;
  499. }
  500.  
  501. octave_idx_type
  502. EIG::init (const Matrix& a, const Matrix& b, bool calc_ev)
  503. {
  504.   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
  505.     (*current_liboctave_error_handler)
  506.       ("EIG: matrix contains Inf or NaN values");
  507.  
  508.   octave_idx_type n = a.rows ();
  509.   octave_idx_type nb = b.rows ();
  510.  
  511.   if (n != a.cols () || nb != b.cols ())
  512.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  513.  
  514.   if (n != nb)
  515.     (*current_liboctave_error_handler) ("EIG requires same size matrices");
  516.  
  517.   octave_idx_type info = 0;
  518.  
  519.   Matrix tmp = b;
  520.   double *tmp_data = tmp.fortran_vec ();
  521.  
  522.   F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
  523.                              n, tmp_data, n,
  524.                              info
  525.                              F77_CHAR_ARG_LEN (1)
  526.                              F77_CHAR_ARG_LEN (1)));
  527.  
  528.   if (a.is_symmetric () && b.is_symmetric () && info == 0)
  529.     return symmetric_init (a, b, calc_ev);
  530.  
  531.   Matrix atmp = a;
  532.   double *atmp_data = atmp.fortran_vec ();
  533.  
  534.   Matrix btmp = b;
  535.   double *btmp_data = btmp.fortran_vec ();
  536.  
  537.   Array<double> ar (dim_vector (n, 1));
  538.   double *par = ar.fortran_vec ();
  539.  
  540.   Array<double> ai (dim_vector (n, 1));
  541.   double *pai = ai.fortran_vec ();
  542.  
  543.   Array<double> beta (dim_vector (n, 1));
  544.   double *pbeta = beta.fortran_vec ();
  545.  
  546.   octave_idx_type tnvr = calc_ev ? n : 0;
  547.   Matrix vr (tnvr, tnvr);
  548.   double *pvr = vr.fortran_vec ();
  549.  
  550.   octave_idx_type lwork = -1;
  551.   double dummy_work;
  552.  
  553.   double *dummy = 0;
  554.   octave_idx_type idummy = 1;
  555.  
  556.   F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  557.                            F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  558.                            n, atmp_data, n, btmp_data, n,
  559.                            par, pai, pbeta,
  560.                            dummy, idummy, pvr, n,
  561.                            &dummy_work, lwork, info
  562.                            F77_CHAR_ARG_LEN (1)
  563.                            F77_CHAR_ARG_LEN (1)));
  564.  
  565.   if (info != 0)
  566.     (*current_liboctave_error_handler) ("dggev workspace query failed");
  567.  
  568.   lwork = static_cast<octave_idx_type> (dummy_work);
  569.   Array<double> work (dim_vector (lwork, 1));
  570.   double *pwork = work.fortran_vec ();
  571.  
  572.   F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  573.                            F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  574.                            n, atmp_data, n, btmp_data, n,
  575.                            par, pai, pbeta,
  576.                            dummy, idummy, pvr, n,
  577.                            pwork, lwork, info
  578.                            F77_CHAR_ARG_LEN (1)
  579.                            F77_CHAR_ARG_LEN (1)));
  580.  
  581.   if (info < 0)
  582.     (*current_liboctave_error_handler) ("unrecoverable error in dggev");
  583.  
  584.   if (info > 0)
  585.     (*current_liboctave_error_handler) ("dggev failed to converge");
  586.  
  587.   lambda.resize (n);
  588.   octave_idx_type nvr = calc_ev ? n : 0;
  589.   v.resize (nvr, nvr);
  590.  
  591.   for (octave_idx_type j = 0; j < n; j++)
  592.     {
  593.       if (ai.elem (j) == 0.0)
  594.         {
  595.           lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j));
  596.           for (octave_idx_type i = 0; i < nvr; i++)
  597.             v.elem (i, j) = vr.elem (i, j);
  598.         }
  599.       else
  600.         {
  601.           if (j+1 >= n)
  602.             (*current_liboctave_error_handler) ("EIG: internal error");
  603.  
  604.           lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j),
  605.                                      ai.elem (j) / beta.elem (j));
  606.           lambda.elem (j+1) = Complex (ar.elem (j+1) / beta.elem (j+1),
  607.                                        ai.elem (j+1) / beta.elem (j+1));
  608.  
  609.           for (octave_idx_type i = 0; i < nvr; i++)
  610.             {
  611.               double real_part = vr.elem (i, j);
  612.               double imag_part = vr.elem (i, j+1);
  613.               v.elem (i, j) = Complex (real_part, imag_part);
  614.               v.elem (i, j+1) = Complex (real_part, -imag_part);
  615.             }
  616.           j++;
  617.         }
  618.     }
  619.  
  620.   return info;
  621. }
  622.  
  623. octave_idx_type
  624. EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev)
  625. {
  626.   octave_idx_type n = a.rows ();
  627.   octave_idx_type nb = b.rows ();
  628.  
  629.   if (n != a.cols () || nb != b.cols ())
  630.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  631.  
  632.   if (n != nb)
  633.     (*current_liboctave_error_handler) ("EIG requires same size matrices");
  634.  
  635.   octave_idx_type info = 0;
  636.  
  637.   Matrix atmp = a;
  638.   double *atmp_data = atmp.fortran_vec ();
  639.  
  640.   Matrix btmp = b;
  641.   double *btmp_data = btmp.fortran_vec ();
  642.  
  643.   ColumnVector wr (n);
  644.   double *pwr = wr.fortran_vec ();
  645.  
  646.   octave_idx_type lwork = -1;
  647.   double dummy_work;
  648.  
  649.   F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  650.                            F77_CONST_CHAR_ARG2 ("U", 1),
  651.                            n, atmp_data, n,
  652.                            btmp_data, n,
  653.                            pwr, &dummy_work, lwork, info
  654.                            F77_CHAR_ARG_LEN (1)
  655.                            F77_CHAR_ARG_LEN (1)));
  656.  
  657.   if (info != 0)
  658.     (*current_liboctave_error_handler) ("dsygv workspace query failed");
  659.  
  660.   lwork = static_cast<octave_idx_type> (dummy_work);
  661.   Array<double> work (dim_vector (lwork, 1));
  662.   double *pwork = work.fortran_vec ();
  663.  
  664.   F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  665.                            F77_CONST_CHAR_ARG2 ("U", 1),
  666.                            n, atmp_data, n,
  667.                            btmp_data, n,
  668.                            pwr, pwork, lwork, info
  669.                            F77_CHAR_ARG_LEN (1)
  670.                            F77_CHAR_ARG_LEN (1)));
  671.  
  672.   if (info < 0)
  673.     (*current_liboctave_error_handler) ("unrecoverable error in dsygv");
  674.  
  675.   if (info > 0)
  676.     (*current_liboctave_error_handler) ("dsygv failed to converge");
  677.  
  678.   lambda = ComplexColumnVector (wr);
  679.   v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
  680.  
  681.   return info;
  682. }
  683.  
  684. octave_idx_type
  685. EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev)
  686. {
  687.   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
  688.     (*current_liboctave_error_handler)
  689.       ("EIG: matrix contains Inf or NaN values");
  690.  
  691.   octave_idx_type n = a.rows ();
  692.   octave_idx_type nb = b.rows ();
  693.  
  694.   if (n != a.cols () || nb != b.cols ())
  695.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  696.  
  697.   if (n != nb)
  698.     (*current_liboctave_error_handler) ("EIG requires same size matrices");
  699.  
  700.   octave_idx_type info = 0;
  701.  
  702.   ComplexMatrix tmp = b;
  703.   Complex*tmp_data = tmp.fortran_vec ();
  704.  
  705.   F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
  706.                              n, tmp_data, n,
  707.                              info
  708.                              F77_CHAR_ARG_LEN (1)
  709.                              F77_CHAR_ARG_LEN (1)));
  710.  
  711.   if (a.is_hermitian () && b.is_hermitian () && info == 0)
  712.     return hermitian_init (a, b, calc_ev);
  713.  
  714.   ComplexMatrix atmp = a;
  715.   Complex *atmp_data = atmp.fortran_vec ();
  716.  
  717.   ComplexMatrix btmp = b;
  718.   Complex *btmp_data = btmp.fortran_vec ();
  719.  
  720.   ComplexColumnVector alpha (n);
  721.   Complex *palpha = alpha.fortran_vec ();
  722.  
  723.   ComplexColumnVector beta (n);
  724.   Complex *pbeta = beta.fortran_vec ();
  725.  
  726.   octave_idx_type nvr = calc_ev ? n : 0;
  727.   ComplexMatrix vtmp (nvr, nvr);
  728.   Complex *pv = vtmp.fortran_vec ();
  729.  
  730.   octave_idx_type lwork = -1;
  731.   Complex dummy_work;
  732.  
  733.   octave_idx_type lrwork = 8*n;
  734.   Array<double> rwork (dim_vector (lrwork, 1));
  735.   double *prwork = rwork.fortran_vec ();
  736.  
  737.   Complex *dummy = 0;
  738.   octave_idx_type idummy = 1;
  739.  
  740.   F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  741.                            F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  742.                            n, atmp_data, n, btmp_data, n,
  743.                            palpha, pbeta, dummy, idummy,
  744.                            pv, n, &dummy_work, lwork, prwork, info
  745.                            F77_CHAR_ARG_LEN (1)
  746.                            F77_CHAR_ARG_LEN (1)));
  747.  
  748.   if (info != 0)
  749.     (*current_liboctave_error_handler) ("zggev workspace query failed");
  750.  
  751.   lwork = static_cast<octave_idx_type> (dummy_work.real ());
  752.   Array<Complex> work (dim_vector (lwork, 1));
  753.   Complex *pwork = work.fortran_vec ();
  754.  
  755.   F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
  756.                            F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  757.                            n, atmp_data, n, btmp_data, n,
  758.                            palpha, pbeta, dummy, idummy,
  759.                            pv, n, pwork, lwork, prwork, info
  760.                            F77_CHAR_ARG_LEN (1)
  761.                            F77_CHAR_ARG_LEN (1)));
  762.  
  763.   if (info < 0)
  764.     (*current_liboctave_error_handler) ("unrecoverable error in zggev");
  765.  
  766.   if (info > 0)
  767.     (*current_liboctave_error_handler) ("zggev failed to converge");
  768.  
  769.   lambda.resize (n);
  770.  
  771.   for (octave_idx_type j = 0; j < n; j++)
  772.     lambda.elem (j) = alpha.elem (j) / beta.elem (j);
  773.  
  774.   v = vtmp;
  775.  
  776.   return info;
  777. }
  778.  
  779. octave_idx_type
  780. EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b,
  781.                      bool calc_ev)
  782. {
  783.   octave_idx_type n = a.rows ();
  784.   octave_idx_type nb = b.rows ();
  785.  
  786.   if (n != a.cols () || nb != b.cols ())
  787.     (*current_liboctave_error_handler) ("EIG requires square matrix");
  788.  
  789.   if (n != nb)
  790.     (*current_liboctave_error_handler) ("EIG requires same size matrices");
  791.  
  792.   octave_idx_type info = 0;
  793.  
  794.   ComplexMatrix atmp = a;
  795.   Complex *atmp_data = atmp.fortran_vec ();
  796.  
  797.   ComplexMatrix btmp = b;
  798.   Complex *btmp_data = btmp.fortran_vec ();
  799.  
  800.   ColumnVector wr (n);
  801.   double *pwr = wr.fortran_vec ();
  802.  
  803.   octave_idx_type lwork = -1;
  804.   Complex dummy_work;
  805.  
  806.   octave_idx_type lrwork = 3*n;
  807.   Array<double> rwork (dim_vector (lrwork, 1));
  808.   double *prwork = rwork.fortran_vec ();
  809.  
  810.   F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  811.                            F77_CONST_CHAR_ARG2 ("U", 1),
  812.                            n, atmp_data, n,
  813.                            btmp_data, n,
  814.                            pwr, &dummy_work, lwork,
  815.                            prwork, info
  816.                            F77_CHAR_ARG_LEN (1)
  817.                            F77_CHAR_ARG_LEN (1)));
  818.  
  819.   if (info != 0)
  820.     (*current_liboctave_error_handler) ("zhegv workspace query failed");
  821.  
  822.   lwork = static_cast<octave_idx_type> (dummy_work.real ());
  823.   Array<Complex> work (dim_vector (lwork, 1));
  824.   Complex *pwork = work.fortran_vec ();
  825.  
  826.   F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
  827.                            F77_CONST_CHAR_ARG2 ("U", 1),
  828.                            n, atmp_data, n,
  829.                            btmp_data, n,
  830.                            pwr, pwork, lwork, prwork, info
  831.                            F77_CHAR_ARG_LEN (1)
  832.                            F77_CHAR_ARG_LEN (1)));
  833.  
  834.   if (info < 0)
  835.     (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
  836.  
  837.   if (info > 0)
  838.     (*current_liboctave_error_handler) ("zhegv failed to converge");
  839.  
  840.   lambda = ComplexColumnVector (wr);
  841.   v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
  842.  
  843.   return info;
  844. }
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!