CmplxSCHUR.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1994-2012 John W. Eaton
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include "CmplxSCHUR.h"
00028 #include "f77-fcn.h"
00029 #include "lo-error.h"
00030 #include "oct-locbuf.h"
00031 
00032 extern "C"
00033 {
00034   F77_RET_T
00035   F77_FUNC (zgeesx, ZGEESX) (F77_CONST_CHAR_ARG_DECL,
00036                              F77_CONST_CHAR_ARG_DECL,
00037                              ComplexSCHUR::select_function,
00038                              F77_CONST_CHAR_ARG_DECL,
00039                              const octave_idx_type&, Complex*,
00040                              const octave_idx_type&, octave_idx_type&,
00041                              Complex*, Complex*, const octave_idx_type&,
00042                              double&, double&, Complex*,
00043                              const octave_idx_type&, double*,
00044                              octave_idx_type*, octave_idx_type&
00045                              F77_CHAR_ARG_LEN_DECL
00046                              F77_CHAR_ARG_LEN_DECL
00047                              F77_CHAR_ARG_LEN_DECL);
00048 
00049   F77_RET_T
00050   F77_FUNC (zrsf2csf, ZRSF2CSF) (const octave_idx_type&, Complex *,
00051                                  Complex *, double *, double *);
00052 }
00053 
00054 static octave_idx_type
00055 select_ana (const Complex& a)
00056 {
00057   return a.real () < 0.0;
00058 }
00059 
00060 static octave_idx_type
00061 select_dig (const Complex& a)
00062 {
00063   return (abs (a) < 1.0);
00064 }
00065 
00066 octave_idx_type
00067 ComplexSCHUR::init (const ComplexMatrix& a, const std::string& ord,
00068                     bool calc_unitary)
00069 {
00070   octave_idx_type a_nr = a.rows ();
00071   octave_idx_type a_nc = a.cols ();
00072 
00073   if (a_nr != a_nc)
00074     {
00075       (*current_liboctave_error_handler)
00076         ("ComplexSCHUR requires square matrix");
00077       return -1;
00078     }
00079   else if (a_nr == 0)
00080     {
00081       schur_mat.clear ();
00082       unitary_mat.clear ();
00083       return 0;
00084     }
00085 
00086   // Workspace requirements may need to be fixed if any of the
00087   // following change.
00088 
00089   char jobvs;
00090   char sense = 'N';
00091   char sort = 'N';
00092 
00093   if (calc_unitary)
00094     jobvs = 'V';
00095   else
00096     jobvs = 'N';
00097 
00098   char ord_char = ord.empty () ? 'U' : ord[0];
00099 
00100   if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
00101     sort = 'S';
00102 
00103   if (ord_char == 'A' || ord_char == 'a')
00104     selector = select_ana;
00105   else if (ord_char == 'D' || ord_char == 'd')
00106     selector = select_dig;
00107   else
00108     selector = 0;
00109 
00110   octave_idx_type n = a_nc;
00111   octave_idx_type lwork = 8 * n;
00112   octave_idx_type info;
00113   octave_idx_type sdim;
00114   double rconde;
00115   double rcondv;
00116 
00117   schur_mat = a;
00118   if (calc_unitary)
00119     unitary_mat.clear (n, n);
00120 
00121   Complex *s = schur_mat.fortran_vec ();
00122   Complex *q = unitary_mat.fortran_vec ();
00123 
00124   Array<double> rwork (dim_vector (n, 1));
00125   double *prwork = rwork.fortran_vec ();
00126 
00127   Array<Complex> w (dim_vector (n, 1));
00128   Complex *pw = w.fortran_vec ();
00129 
00130   Array<Complex> work (dim_vector (lwork, 1));
00131   Complex *pwork = work.fortran_vec ();
00132 
00133   // BWORK is not referenced for non-ordered Schur.
00134   octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
00135   Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
00136   octave_idx_type *pbwork = bwork.fortran_vec ();
00137 
00138   F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
00139                              F77_CONST_CHAR_ARG2 (&sort, 1),
00140                              selector,
00141                              F77_CONST_CHAR_ARG2 (&sense, 1),
00142                              n, s, n, sdim, pw, q, n, rconde, rcondv,
00143                              pwork, lwork, prwork, pbwork, info
00144                              F77_CHAR_ARG_LEN (1)
00145                              F77_CHAR_ARG_LEN (1)
00146                              F77_CHAR_ARG_LEN (1)));
00147 
00148   return info;
00149 }
00150 
00151 ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& s, const ComplexMatrix& u)
00152   : schur_mat (s), unitary_mat (u), selector (0)
00153 {
00154   octave_idx_type n = s.rows ();
00155   if (s.columns () != n || u.rows () != n || u.columns () != n)
00156     (*current_liboctave_error_handler)
00157       ("schur: inconsistent matrix dimensions");
00158 }
00159 
00160 ComplexSCHUR::ComplexSCHUR (const SCHUR& s)
00161   : schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ()),
00162     selector (0)
00163 {
00164   octave_idx_type n = schur_mat.rows ();
00165   if (n > 0)
00166     {
00167       OCTAVE_LOCAL_BUFFER (double, c, n-1);
00168       OCTAVE_LOCAL_BUFFER (double, sx, n-1);
00169 
00170       F77_XFCN (zrsf2csf, ZRSF2CSF, (n, schur_mat.fortran_vec (),
00171                                      unitary_mat.fortran_vec (), c, sx));
00172     }
00173 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines