1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1994-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if defined (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include "CMatrix.h" 31 #include "aepbalance.h" 32 #include "dColVector.h" 33 #include "dMatrix.h" 34 #include "fCMatrix.h" 35 #include "fColVector.h" 36 #include "fMatrix.h" 37 #include "lo-error.h" 38 #include "lo-lapack-proto.h" 39 40 namespace octave 41 { 42 static inline char get_job(bool noperm,bool noscal)43 get_job (bool noperm, bool noscal) 44 { 45 return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); 46 } 47 48 namespace math 49 { 50 template <> aepbalance(const Matrix & a,bool noperm,bool noscal)51 aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal) 52 : balanced_mat (a), scale (), ilo (), ihi (), 53 job (get_job (noperm, noscal)) 54 { 55 F77_INT n = to_f77_int (a.cols ()); 56 57 if (a.rows () != n) 58 (*current_liboctave_error_handler) 59 ("aepbalance: requires square matrix"); 60 61 scale = ColumnVector (n); 62 63 F77_INT info, t_ilo, t_ihi; 64 65 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, 66 balanced_mat.fortran_vec (), n, 67 t_ilo, t_ihi, scale.fortran_vec (), info 68 F77_CHAR_ARG_LEN (1))); 69 70 ilo = t_ilo; 71 ihi = t_ihi; 72 } 73 74 template <> 75 Matrix balancing_matrix(void) const76 aepbalance<Matrix>::balancing_matrix (void) const 77 { 78 F77_INT n = to_f77_int (balanced_mat.rows ()); 79 80 Matrix balancing_mat (n, n, 0.0); 81 for (F77_INT i = 0; i < n; i++) 82 balancing_mat.elem (i ,i) = 1.0; 83 84 F77_INT info; 85 F77_INT t_ilo = to_f77_int (ilo); 86 F77_INT t_ihi = to_f77_int (ihi); 87 88 char side = 'R'; 89 90 F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), 91 F77_CONST_CHAR_ARG2 (&side, 1), 92 n, t_ilo, t_ihi, scale.data (), n, 93 balancing_mat.fortran_vec (), n, info 94 F77_CHAR_ARG_LEN (1) 95 F77_CHAR_ARG_LEN (1))); 96 97 return balancing_mat; 98 } 99 100 template <> aepbalance(const FloatMatrix & a,bool noperm,bool noscal)101 aepbalance<FloatMatrix>::aepbalance (const FloatMatrix& a, bool noperm, 102 bool noscal) 103 : balanced_mat (a), scale (), ilo (), ihi (), 104 job (get_job (noperm, noscal)) 105 { 106 F77_INT n = to_f77_int (a.cols ()); 107 108 if (a.rows () != n) 109 (*current_liboctave_error_handler) 110 ("aepbalance: requires square matrix"); 111 112 scale = FloatColumnVector (n); 113 114 F77_INT info, t_ilo, t_ihi; 115 116 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, 117 balanced_mat.fortran_vec (), n, t_ilo, 118 t_ihi, scale.fortran_vec (), info 119 F77_CHAR_ARG_LEN (1))); 120 121 ilo = t_ilo; 122 ihi = t_ihi; 123 } 124 125 template <> 126 FloatMatrix balancing_matrix(void) const127 aepbalance<FloatMatrix>::balancing_matrix (void) const 128 { 129 F77_INT n = to_f77_int (balanced_mat.rows ()); 130 131 FloatMatrix balancing_mat (n, n, 0.0); 132 for (F77_INT i = 0; i < n; i++) 133 balancing_mat.elem (i,i) = 1.0; 134 135 F77_INT info; 136 F77_INT t_ilo = to_f77_int (ilo); 137 F77_INT t_ihi = to_f77_int (ihi); 138 139 char side = 'R'; 140 141 F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), 142 F77_CONST_CHAR_ARG2 (&side, 1), 143 n, t_ilo, t_ihi, scale.data (), n, 144 balancing_mat.fortran_vec (), n, info 145 F77_CHAR_ARG_LEN (1) 146 F77_CHAR_ARG_LEN (1))); 147 148 return balancing_mat; 149 } 150 151 template <> aepbalance(const ComplexMatrix & a,bool noperm,bool noscal)152 aepbalance<ComplexMatrix>::aepbalance (const ComplexMatrix& a, bool noperm, 153 bool noscal) 154 : balanced_mat (a), scale (), ilo (), ihi (), 155 job (get_job (noperm, noscal)) 156 { 157 F77_INT n = to_f77_int (a.cols ()); 158 159 if (a.rows () != n) 160 (*current_liboctave_error_handler) 161 ("aepbalance: requires square matrix"); 162 163 scale = ColumnVector (n); 164 165 F77_INT info, t_ilo, t_ihi; 166 167 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, 168 F77_DBLE_CMPLX_ARG (balanced_mat.fortran_vec ()), 169 n, t_ilo, t_ihi, scale.fortran_vec (), info 170 F77_CHAR_ARG_LEN (1))); 171 172 ilo = t_ilo; 173 ihi = t_ihi; 174 } 175 176 template <> 177 ComplexMatrix balancing_matrix(void) const178 aepbalance<ComplexMatrix>::balancing_matrix (void) const 179 { 180 F77_INT n = to_f77_int (balanced_mat.rows ()); 181 182 ComplexMatrix balancing_mat (n, n, 0.0); 183 for (F77_INT i = 0; i < n; i++) 184 balancing_mat.elem (i, i) = 1.0; 185 186 F77_INT info; 187 F77_INT t_ilo = to_f77_int (ilo); 188 F77_INT t_ihi = to_f77_int (ihi); 189 190 char side = 'R'; 191 192 F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), 193 F77_CONST_CHAR_ARG2 (&side, 1), 194 n, t_ilo, t_ihi, scale.data (), n, 195 F77_DBLE_CMPLX_ARG (balancing_mat.fortran_vec ()), 196 n, info 197 F77_CHAR_ARG_LEN (1) 198 F77_CHAR_ARG_LEN (1))); 199 200 return balancing_mat; 201 } 202 203 template <> aepbalance(const FloatComplexMatrix & a,bool noperm,bool noscal)204 aepbalance<FloatComplexMatrix>::aepbalance (const FloatComplexMatrix& a, 205 bool noperm, bool noscal) 206 : balanced_mat (a), scale (), ilo (), ihi (), 207 job (get_job (noperm, noscal)) 208 { 209 F77_INT n = to_f77_int (a.cols ()); 210 211 if (a.rows () != n) 212 (*current_liboctave_error_handler) 213 ("aepbalance: requires square matrix"); 214 215 scale = FloatColumnVector (n); 216 217 F77_INT info, t_ilo, t_ihi; 218 219 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, 220 F77_CMPLX_ARG (balanced_mat.fortran_vec ()), 221 n, t_ilo, t_ihi, scale.fortran_vec (), info 222 F77_CHAR_ARG_LEN (1))); 223 224 ilo = t_ilo; 225 ihi = t_ihi; 226 } 227 228 template <> 229 FloatComplexMatrix balancing_matrix(void) const230 aepbalance<FloatComplexMatrix>::balancing_matrix (void) const 231 { 232 F77_INT n = to_f77_int (balanced_mat.rows ()); 233 234 FloatComplexMatrix balancing_mat (n, n, 0.0); 235 for (F77_INT i = 0; i < n; i++) 236 balancing_mat.elem (i, i) = 1.0; 237 238 F77_INT info; 239 F77_INT t_ilo = to_f77_int (ilo); 240 F77_INT t_ihi = to_f77_int (ihi); 241 242 char side = 'R'; 243 244 F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), 245 F77_CONST_CHAR_ARG2 (&side, 1), 246 n, t_ilo, t_ihi, scale.data (), n, 247 F77_CMPLX_ARG (balancing_mat.fortran_vec ()), 248 n, info 249 F77_CHAR_ARG_LEN (1) 250 F77_CHAR_ARG_LEN (1))); 251 252 return balancing_mat; 253 } 254 255 // Instantiations we need. 256 257 template class aepbalance<Matrix>; 258 259 template class aepbalance<FloatMatrix>; 260 261 template class aepbalance<ComplexMatrix>; 262 263 template class aepbalance<FloatComplexMatrix>; 264 } 265 } 266