CmplxGEPBAL.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 <string>
00028 #include <vector>
00029 
00030 #include "CmplxGEPBAL.h"
00031 #include "Array-util.h"
00032 #include "f77-fcn.h"
00033 #include "oct-locbuf.h"
00034 
00035 extern "C"
00036 {
00037   F77_RET_T
00038   F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
00039                              const octave_idx_type& N, Complex* A,
00040                              const octave_idx_type& LDA, Complex* B,
00041                              const octave_idx_type& LDB,
00042                              octave_idx_type& ILO, octave_idx_type& IHI,
00043                              double* LSCALE, double* RSCALE,
00044                              double* WORK, octave_idx_type& INFO
00045                              F77_CHAR_ARG_LEN_DECL);
00046 
00047   F77_RET_T
00048   F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
00049                              F77_CONST_CHAR_ARG_DECL,
00050                              const octave_idx_type& N,
00051                              const octave_idx_type& ILO,
00052                              const octave_idx_type& IHI,
00053                              const double* LSCALE, const double* RSCALE,
00054                              octave_idx_type& M, double* V,
00055                              const octave_idx_type& LDV, octave_idx_type& INFO
00056                              F77_CHAR_ARG_LEN_DECL
00057                              F77_CHAR_ARG_LEN_DECL);
00058 
00059 }
00060 
00061 octave_idx_type
00062 ComplexGEPBALANCE::init (const ComplexMatrix& a, const ComplexMatrix& b,
00063                   const std::string& balance_job)
00064 {
00065   octave_idx_type n = a.cols ();
00066 
00067   if (a.rows () != n)
00068     {
00069       (*current_liboctave_error_handler) ("ComplexGEPBALANCE requires square matrix");
00070       return -1;
00071     }
00072 
00073   if (a.dims() != b.dims ())
00074     {
00075       gripe_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols());
00076       return -1;
00077     }
00078 
00079   octave_idx_type info;
00080   octave_idx_type ilo;
00081   octave_idx_type ihi;
00082 
00083   OCTAVE_LOCAL_BUFFER (double, plscale, n);
00084   OCTAVE_LOCAL_BUFFER (double, prscale,  n);
00085   OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
00086 
00087   balanced_mat = a;
00088   Complex *p_balanced_mat = balanced_mat.fortran_vec ();
00089   balanced_mat2 = b;
00090   Complex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
00091 
00092   char job = balance_job[0];
00093 
00094   F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
00095                              n, p_balanced_mat, n, p_balanced_mat2,
00096                              n, ilo, ihi, plscale, prscale, pwork, info
00097                              F77_CHAR_ARG_LEN (1)));
00098 
00099   balancing_mat = Matrix (n, n, 0.0);
00100   balancing_mat2 = Matrix (n, n, 0.0);
00101   for (octave_idx_type i = 0; i < n; i++)
00102     {
00103       octave_quit ();
00104       balancing_mat.elem (i ,i) = 1.0;
00105       balancing_mat2.elem (i ,i) = 1.0;
00106     }
00107 
00108   double *p_balancing_mat = balancing_mat.fortran_vec ();
00109   double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
00110 
00111   // first left
00112   F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
00113                              F77_CONST_CHAR_ARG2 ("L", 1),
00114                              n, ilo, ihi, plscale, prscale,
00115                              n, p_balancing_mat, n, info
00116                              F77_CHAR_ARG_LEN (1)
00117                              F77_CHAR_ARG_LEN (1)));
00118 
00119   // then right
00120   F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
00121                              F77_CONST_CHAR_ARG2 ("R", 1),
00122                              n, ilo, ihi, plscale, prscale,
00123                              n, p_balancing_mat2, n, info
00124                              F77_CHAR_ARG_LEN (1)
00125                              F77_CHAR_ARG_LEN (1)));
00126 
00127   return info;
00128 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines