Mercurial > octave

comparison liboctave/numeric/aepbalance.cc @ 21268:f08ae27289e4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
better use of templates for balance classes * aepbalance.h, aepbalance.cc: New files generated from base-aepbal.h, CmplxAEPBAL.cc, CmplxAEPBAL.h, dbleAEPBAL.cc, dbleAEPBAL.h, fCmplxAEPBAL.cc, fCmplxAEPBAL.h, floatAEPBAL.cc, and floatAEPBAL.h and making them templates. * gepbalance.h, gepbalance.cc: New files generate from CmplxGEPBAL.cc, CmplxGEPBAL.h, dbleGEPBAL.cc, dbleGEPBAL.h, fCmplxGEPBAL.cc, fCmplxGEPBAL.h, floatGEPBAL.cc, and floatGEPBAL.h and making them templates. * liboctave/numeric/module.mk: Update. * balance.cc, CMatrix.h, dMatrix.h, fCMatrix.h, fMatrix.h, mx-defs.h, mx-ext.h: Use new classes.
author John W. Eaton <jwe@octave.org>
date Tue, 16 Feb 2016 00:32:29 -0500
parents liboctave/numeric/dbleAEPBAL.cc@f7121e111991
children 40de9f8f23a6
comparison
equal deleted inserted replaced
21267:f5b8c3aca5f8 21268:f08ae27289e4
1 /*
2
3 Copyright (C) 1994-2015 John W. Eaton
4 Copyright (C) 2008 Jaroslav Hajek
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 # include <config.h>
26 #endif
27
28 #include <string>
29
30 #include "CColVector.h"
31 #include "CMatrix.h"
32 #include "aepbalance.h"
33 #include "dColVector.h"
34 #include "dMatrix.h"
35 #include "f77-fcn.h"
36 #include "fCColVector.h"
37 #include "fCMatrix.h"
38 #include "fColVector.h"
39 #include "fMatrix.h"
40
41 extern "C"
42 {
43 F77_RET_T
44 F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
45 const octave_idx_type&, double*,
46 const octave_idx_type&, octave_idx_type&,
47 octave_idx_type&, double*, octave_idx_type&
48 F77_CHAR_ARG_LEN_DECL);
49
50 F77_RET_T
51 F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
52 F77_CONST_CHAR_ARG_DECL,
53 const octave_idx_type&, const octave_idx_type&,
54 const octave_idx_type&, const double*,
55 const octave_idx_type&, double*,
56 const octave_idx_type&, octave_idx_type&
57 F77_CHAR_ARG_LEN_DECL
58 F77_CHAR_ARG_LEN_DECL);
59
60 F77_RET_T
61 F77_FUNC (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL,
62 const octave_idx_type&, float*,
63 const octave_idx_type&, octave_idx_type&,
64 octave_idx_type&, float*, octave_idx_type&
65 F77_CHAR_ARG_LEN_DECL);
66
67 F77_RET_T
68 F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL,
69 F77_CONST_CHAR_ARG_DECL,
70 const octave_idx_type&, const octave_idx_type&,
71 const octave_idx_type&, const float*,
72 const octave_idx_type&, float*,
73 const octave_idx_type&, octave_idx_type&
74 F77_CHAR_ARG_LEN_DECL
75 F77_CHAR_ARG_LEN_DECL);
76
77 F77_RET_T
78 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
79 const octave_idx_type&, Complex*,
80 const octave_idx_type&, octave_idx_type&,
81 octave_idx_type&, double*, octave_idx_type&
82 F77_CHAR_ARG_LEN_DECL);
83
84 F77_RET_T
85 F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL,
86 F77_CONST_CHAR_ARG_DECL,
87 const octave_idx_type&, const octave_idx_type&,
88 const octave_idx_type&, const double*,
89 const octave_idx_type&, Complex*,
90 const octave_idx_type&, octave_idx_type&
91 F77_CHAR_ARG_LEN_DECL
92 F77_CHAR_ARG_LEN_DECL);
93
94 F77_RET_T
95 F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL,
96 const octave_idx_type&, FloatComplex*,
97 const octave_idx_type&, octave_idx_type&,
98 octave_idx_type&, float*, octave_idx_type&
99 F77_CHAR_ARG_LEN_DECL);
100
101 F77_RET_T
102 F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL,
103 F77_CONST_CHAR_ARG_DECL,
104 const octave_idx_type&, const octave_idx_type&,
105 const octave_idx_type&, const float*,
106 const octave_idx_type&, FloatComplex*,
107 const octave_idx_type&, octave_idx_type&
108 F77_CHAR_ARG_LEN_DECL
109 F77_CHAR_ARG_LEN_DECL);
110 }
111
112 static inline char
113 get_job (bool noperm, bool noscal)
114 {
115 return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
116 }
117
118 template <>
119 aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal)
120 : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
121 {
122 octave_idx_type n = a.cols ();
123
124 if (a.rows () != n)
125 (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
126
127 scale = ColumnVector (n);
128
129 octave_idx_type info;
130
131 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
132 balanced_mat.fortran_vec (), n, ilo, ihi,
133 scale.fortran_vec (), info
134 F77_CHAR_ARG_LEN (1)));
135 }
136
137 template <>
138 Matrix
139 aepbalance<Matrix>::balancing_matrix (void) const
140 {
141 octave_idx_type n = balanced_mat.rows ();
142 Matrix balancing_mat (n, n, 0.0);
143 for (octave_idx_type i = 0; i < n; i++)
144 balancing_mat.elem (i ,i) = 1.0;
145
146 octave_idx_type info;
147
148 char side = 'R';
149
150 F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
151 F77_CONST_CHAR_ARG2 (&side, 1),
152 n, ilo, ihi, scale.data (), n,
153 balancing_mat.fortran_vec (), n, info
154 F77_CHAR_ARG_LEN (1)
155 F77_CHAR_ARG_LEN (1)));
156
157 return balancing_mat;
158 }
159
160 template <>
161 aepbalance<FloatMatrix>::aepbalance (const FloatMatrix& a, bool noperm,
162 bool noscal)
163 : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
164 {
165 octave_idx_type n = a.cols ();
166
167 if (a.rows () != n)
168 (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
169
170 scale = FloatColumnVector (n);
171
172 octave_idx_type info;
173
174 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
175 balanced_mat.fortran_vec (), n, ilo, ihi,
176 scale.fortran_vec (), info
177 F77_CHAR_ARG_LEN (1)));
178 }
179
180 template <>
181 FloatMatrix
182 aepbalance<FloatMatrix>::balancing_matrix (void) const
183 {
184 octave_idx_type n = balanced_mat.rows ();
185 FloatMatrix balancing_mat (n, n, 0.0);
186 for (octave_idx_type i = 0; i < n; i++)
187 balancing_mat.elem (i ,i) = 1.0;
188
189 octave_idx_type info;
190
191 char side = 'R';
192
193 F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
194 F77_CONST_CHAR_ARG2 (&side, 1),
195 n, ilo, ihi, scale.data (), n,
196 balancing_mat.fortran_vec (), n, info
197 F77_CHAR_ARG_LEN (1)
198 F77_CHAR_ARG_LEN (1)));
199
200 return balancing_mat;
201 }
202
203 template <>
204 aepbalance<ComplexMatrix>::aepbalance (const ComplexMatrix& a, bool noperm,
205 bool noscal)
206 : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
207 {
208 octave_idx_type n = a.cols ();
209
210 if (a.rows () != n)
211 (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
212
213 scale = ColumnVector (n);
214
215 octave_idx_type info;
216
217 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
218 balanced_mat.fortran_vec (), n, ilo, ihi,
219 scale.fortran_vec (), info
220 F77_CHAR_ARG_LEN (1)));
221 }
222
223 template <>
224 ComplexMatrix
225 aepbalance<ComplexMatrix>::balancing_matrix (void) const
226 {
227 octave_idx_type n = balanced_mat.rows ();
228 ComplexMatrix balancing_mat (n, n, 0.0);
229 for (octave_idx_type i = 0; i < n; i++)
230 balancing_mat.elem (i, i) = 1.0;
231
232 octave_idx_type info;
233
234 char side = 'R';
235
236 F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
237 F77_CONST_CHAR_ARG2 (&side, 1),
238 n, ilo, ihi, scale.data (), n,
239 balancing_mat.fortran_vec (), n, info
240 F77_CHAR_ARG_LEN (1)
241 F77_CHAR_ARG_LEN (1)));
242
243 return balancing_mat;
244 }
245
246 template <>
247 aepbalance<FloatComplexMatrix>::aepbalance (const FloatComplexMatrix& a,
248 bool noperm, bool noscal)
249 : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
250 {
251 octave_idx_type n = a.cols ();
252
253 if (a.rows () != n)
254 (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
255
256 scale = FloatColumnVector (n);
257
258 octave_idx_type info;
259
260 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
261 balanced_mat.fortran_vec (), n, ilo, ihi,
262 scale.fortran_vec (), info
263 F77_CHAR_ARG_LEN (1)));
264 }
265
266 template <>
267 FloatComplexMatrix
268 aepbalance<FloatComplexMatrix>::balancing_matrix (void) const
269 {
270 octave_idx_type n = balanced_mat.rows ();
271 FloatComplexMatrix balancing_mat (n, n, 0.0);
272 for (octave_idx_type i = 0; i < n; i++)
273 balancing_mat.elem (i, i) = 1.0;
274
275 octave_idx_type info;
276
277 char side = 'R';
278
279 F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
280 F77_CONST_CHAR_ARG2 (&side, 1),
281 n, ilo, ihi, scale.data (), n,
282 balancing_mat.fortran_vec (), n, info
283 F77_CHAR_ARG_LEN (1)
284 F77_CHAR_ARG_LEN (1)));
285
286 return balancing_mat;
287 }
288
289 // Instantiations we need.
290
291 template class aepbalance<Matrix>;
292
293 template class aepbalance<FloatMatrix>;
294
295 template class aepbalance<ComplexMatrix>;
296
297 template class aepbalance<FloatComplexMatrix>;