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