Mercurial > octave

comparison liboctave/numeric/gepbalance.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/dbleGEPBAL.cc@f7121e111991
children 40de9f8f23a6
comparison
equal deleted inserted replaced
21267:f5b8c3aca5f8 21268:f08ae27289e4
1 /*
2
3 Copyright (C) 1994-2015 John W. Eaton
4
5 This file is part of Octave.
6
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20
21 */
22
23 #ifdef HAVE_CONFIG_H
24 # include <config.h>
25 #endif
26
27 #include <string>
28 #include <vector>
29
30 #include "Array-util.h"
31 #include "CMatrix.h"
32 #include "dMatrix.h"
33 #include "f77-fcn.h"
34 #include "fCMatrix.h"
35 #include "fMatrix.h"
36 #include "gepbalance.h"
37 #include "oct-locbuf.h"
38
39 extern "C"
40 {
41 F77_RET_T
42 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL,
43 const octave_idx_type& N,
44 double* A, const octave_idx_type& LDA,
45 double* B, const octave_idx_type& LDB,
46 octave_idx_type& ILO, octave_idx_type& IHI,
47 double* LSCALE, double* RSCALE,
48 double* WORK, octave_idx_type& INFO
49 F77_CHAR_ARG_LEN_DECL);
50
51 F77_RET_T
52 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
53 F77_CONST_CHAR_ARG_DECL,
54 const octave_idx_type& N,
55 const octave_idx_type& ILO,
56 const octave_idx_type& IHI,
57 const double* LSCALE, const double* RSCALE,
58 octave_idx_type& M, double* V,
59 const octave_idx_type& LDV, octave_idx_type& INFO
60 F77_CHAR_ARG_LEN_DECL
61 F77_CHAR_ARG_LEN_DECL);
62
63 F77_RET_T
64 F77_FUNC (sggbal, SGGBAL) (F77_CONST_CHAR_ARG_DECL,
65 const octave_idx_type& N, float* A,
66 const octave_idx_type& LDA, float* B,
67 const octave_idx_type& LDB,
68 octave_idx_type& ILO, octave_idx_type& IHI,
69 float* LSCALE, float* RSCALE,
70 float* WORK, octave_idx_type& INFO
71 F77_CHAR_ARG_LEN_DECL);
72
73 F77_RET_T
74 F77_FUNC (sggbak, SGGBAK) (F77_CONST_CHAR_ARG_DECL,
75 F77_CONST_CHAR_ARG_DECL,
76 const octave_idx_type& N,
77 const octave_idx_type& ILO,
78 const octave_idx_type& IHI,
79 const float* LSCALE, const float* RSCALE,
80 octave_idx_type& M, float* V,
81 const octave_idx_type& LDV, octave_idx_type& INFO
82 F77_CHAR_ARG_LEN_DECL
83 F77_CHAR_ARG_LEN_DECL);
84
85 F77_RET_T
86 F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
87 const octave_idx_type& N, Complex* A,
88 const octave_idx_type& LDA, Complex* B,
89 const octave_idx_type& LDB,
90 octave_idx_type& ILO, octave_idx_type& IHI,
91 double* LSCALE, double* RSCALE,
92 double* WORK, octave_idx_type& INFO
93 F77_CHAR_ARG_LEN_DECL);
94
95 F77_RET_T
96 F77_FUNC (cggbal, CGGBAL) (F77_CONST_CHAR_ARG_DECL,
97 const octave_idx_type& N,
98 FloatComplex* A, const octave_idx_type& LDA,
99 FloatComplex* B, const octave_idx_type& LDB,
100 octave_idx_type& ILO, octave_idx_type& IHI,
101 float* LSCALE, float* RSCALE,
102 float* WORK, octave_idx_type& INFO
103 F77_CHAR_ARG_LEN_DECL);
104 }
105
106 template <>
107 octave_idx_type
108 gepbalance<Matrix>::init (const Matrix& a, const Matrix& b,
109 const std::string& balance_job)
110 {
111 octave_idx_type n = a.cols ();
112
113 if (a.rows () != n)
114 (*current_liboctave_error_handler) ("GEPBALANCE requires square matrix");
115
116 if (a.dims () != b.dims ())
117 err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
118
119 octave_idx_type info;
120 octave_idx_type ilo;
121 octave_idx_type ihi;
122
123 OCTAVE_LOCAL_BUFFER (double, plscale, n);
124 OCTAVE_LOCAL_BUFFER (double, prscale, n);
125 OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
126
127 balanced_mat = a;
128 double *p_balanced_mat = balanced_mat.fortran_vec ();
129 balanced_mat2 = b;
130 double *p_balanced_mat2 = balanced_mat2.fortran_vec ();
131
132 char job = balance_job[0];
133
134 F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
135 n, p_balanced_mat, n, p_balanced_mat2,
136 n, ilo, ihi, plscale, prscale, pwork, info
137 F77_CHAR_ARG_LEN (1)));
138
139 balancing_mat = Matrix (n, n, 0.0);
140 balancing_mat2 = Matrix (n, n, 0.0);
141 for (octave_idx_type i = 0; i < n; i++)
142 {
143 octave_quit ();
144 balancing_mat.elem (i ,i) = 1.0;
145 balancing_mat2.elem (i ,i) = 1.0;
146 }
147
148 double *p_balancing_mat = balancing_mat.fortran_vec ();
149 double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
150
151 // first left
152 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
153 F77_CONST_CHAR_ARG2 ("L", 1),
154 n, ilo, ihi, plscale, prscale,
155 n, p_balancing_mat, n, info
156 F77_CHAR_ARG_LEN (1)
157 F77_CHAR_ARG_LEN (1)));
158
159 // then right
160 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
161 F77_CONST_CHAR_ARG2 ("R", 1),
162 n, ilo, ihi, plscale, prscale,
163 n, p_balancing_mat2, n, info
164 F77_CHAR_ARG_LEN (1)
165 F77_CHAR_ARG_LEN (1)));
166
167 return info;
168 }
169
170 template <>
171 octave_idx_type
172 gepbalance<FloatMatrix>::init (const FloatMatrix& a, const FloatMatrix& b,
173 const std::string& balance_job)
174 {
175 octave_idx_type n = a.cols ();
176
177 if (a.rows () != n)
178 (*current_liboctave_error_handler)
179 ("FloatGEPBALANCE requires square matrix");
180
181 if (a.dims () != b.dims ())
182 err_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols());
183
184 octave_idx_type info;
185 octave_idx_type ilo;
186 octave_idx_type ihi;
187
188 OCTAVE_LOCAL_BUFFER (float, plscale, n);
189 OCTAVE_LOCAL_BUFFER (float, prscale, n);
190 OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
191
192 balanced_mat = a;
193 float *p_balanced_mat = balanced_mat.fortran_vec ();
194 balanced_mat2 = b;
195 float *p_balanced_mat2 = balanced_mat2.fortran_vec ();
196
197 char job = balance_job[0];
198
199 F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
200 n, p_balanced_mat, n, p_balanced_mat2,
201 n, ilo, ihi, plscale, prscale, pwork, info
202 F77_CHAR_ARG_LEN (1)));
203
204 balancing_mat = FloatMatrix (n, n, 0.0);
205 balancing_mat2 = FloatMatrix (n, n, 0.0);
206 for (octave_idx_type i = 0; i < n; i++)
207 {
208 octave_quit ();
209 balancing_mat.elem (i ,i) = 1.0;
210 balancing_mat2.elem (i ,i) = 1.0;
211 }
212
213 float *p_balancing_mat = balancing_mat.fortran_vec ();
214 float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
215
216 // first left
217 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
218 F77_CONST_CHAR_ARG2 ("L", 1),
219 n, ilo, ihi, plscale, prscale,
220 n, p_balancing_mat, n, info
221 F77_CHAR_ARG_LEN (1)
222 F77_CHAR_ARG_LEN (1)));
223
224 // then right
225 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
226 F77_CONST_CHAR_ARG2 ("R", 1),
227 n, ilo, ihi, plscale, prscale,
228 n, p_balancing_mat2, n, info
229 F77_CHAR_ARG_LEN (1)
230 F77_CHAR_ARG_LEN (1)));
231
232 return info;
233 }
234
235 template <>
236 octave_idx_type
237 gepbalance<ComplexMatrix>::init (const ComplexMatrix& a,
238 const ComplexMatrix& b,
239 const std::string& balance_job)
240 {
241 octave_idx_type n = a.cols ();
242
243 if (a.rows () != n)
244 (*current_liboctave_error_handler)
245 ("ComplexGEPBALANCE requires square matrix");
246
247 if (a.dims () != b.dims ())
248 err_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols());
249
250 octave_idx_type info;
251 octave_idx_type ilo;
252 octave_idx_type ihi;
253
254 OCTAVE_LOCAL_BUFFER (double, plscale, n);
255 OCTAVE_LOCAL_BUFFER (double, prscale, n);
256 OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
257
258 balanced_mat = a;
259 Complex *p_balanced_mat = balanced_mat.fortran_vec ();
260 balanced_mat2 = b;
261 Complex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
262
263 char job = balance_job[0];
264
265 F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
266 n, p_balanced_mat, n, p_balanced_mat2,
267 n, ilo, ihi, plscale, prscale, pwork, info
268 F77_CHAR_ARG_LEN (1)));
269
270 balancing_mat = Matrix (n, n, 0.0);
271 balancing_mat2 = Matrix (n, n, 0.0);
272 for (octave_idx_type i = 0; i < n; i++)
273 {
274 octave_quit ();
275 balancing_mat.elem (i ,i) = 1.0;
276 balancing_mat2.elem (i ,i) = 1.0;
277 }
278
279 double *p_balancing_mat = balancing_mat.fortran_vec ();
280 double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
281
282 // first left
283 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
284 F77_CONST_CHAR_ARG2 ("L", 1),
285 n, ilo, ihi, plscale, prscale,
286 n, p_balancing_mat, n, info
287 F77_CHAR_ARG_LEN (1)
288 F77_CHAR_ARG_LEN (1)));
289
290 // then right
291 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
292 F77_CONST_CHAR_ARG2 ("R", 1),
293 n, ilo, ihi, plscale, prscale,
294 n, p_balancing_mat2, n, info
295 F77_CHAR_ARG_LEN (1)
296 F77_CHAR_ARG_LEN (1)));
297
298 return info;
299 }
300
301 template <>
302 octave_idx_type
303 gepbalance<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
304 const FloatComplexMatrix& b,
305 const std::string& balance_job)
306 {
307 octave_idx_type n = a.cols ();
308
309 if (a.rows () != n)
310 {
311 (*current_liboctave_error_handler)
312 ("FloatComplexGEPBALANCE requires square matrix");
313 return -1;
314 }
315
316 if (a.dims () != b.dims ())
317 err_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols());
318
319 octave_idx_type info;
320 octave_idx_type ilo;
321 octave_idx_type ihi;
322
323 OCTAVE_LOCAL_BUFFER (float, plscale, n);
324 OCTAVE_LOCAL_BUFFER (float, prscale, n);
325 OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
326
327 balanced_mat = a;
328 FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
329 balanced_mat2 = b;
330 FloatComplex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
331
332 char job = balance_job[0];
333
334 F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
335 n, p_balanced_mat, n, p_balanced_mat2,
336 n, ilo, ihi, plscale,prscale, pwork, info
337 F77_CHAR_ARG_LEN (1)));
338
339 balancing_mat = FloatMatrix (n, n, 0.0);
340 balancing_mat2 = FloatMatrix (n, n, 0.0);
341 for (octave_idx_type i = 0; i < n; i++)
342 {
343 octave_quit ();
344 balancing_mat.elem (i ,i) = 1.0;
345 balancing_mat2.elem (i ,i) = 1.0;
346 }
347
348 float *p_balancing_mat = balancing_mat.fortran_vec ();
349 float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
350
351 // first left
352 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
353 F77_CONST_CHAR_ARG2 ("L", 1),
354 n, ilo, ihi, plscale, prscale,
355 n, p_balancing_mat, n, info
356 F77_CHAR_ARG_LEN (1)
357 F77_CHAR_ARG_LEN (1)));
358
359 // then right
360 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
361 F77_CONST_CHAR_ARG2 ("R", 1),
362 n, ilo, ihi, plscale, prscale,
363 n, p_balancing_mat2, n, info
364 F77_CHAR_ARG_LEN (1)
365 F77_CHAR_ARG_LEN (1)));
366
367 return info;
368 }
369
370 // Instantiations we need.
371
372 template class gepbalance<Matrix>;
373
374 template class gepbalance<FloatMatrix>;
375
376 template class gepbalance<ComplexMatrix>;
377
378 template class gepbalance<FloatComplexMatrix>;