Mercurial > octave

comparison liboctave/numeric/schur.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
39 39
40 OCTAVE_BEGIN_NAMESPACE(octave) 40 OCTAVE_BEGIN_NAMESPACE(octave)
41 41
42 OCTAVE_BEGIN_NAMESPACE(math) 42 OCTAVE_BEGIN_NAMESPACE(math)
43 43
44 // For real types. 44 // For real types.
45 45
46 static F77_INT 46 static F77_INT
47 select_ana (const double& a, const double&) 47 select_ana (const double& a, const double&)
48 {
49 return (a < 0.0);
50 }
51
52 static F77_INT
53 select_dig (const double& a, const double& b)
54 {
55 return (hypot (a, b) < 1.0);
56 }
57
58 static F77_INT
59 select_ana (const float& a, const float&)
60 {
61 return (a < 0.0);
62 }
63
64 static F77_INT
65 select_dig (const float& a, const float& b)
66 {
67 return (hypot (a, b) < 1.0);
68 }
69
70 // For complex types.
71
72 static F77_INT
73 select_ana (const F77_DBLE_CMPLX& a_arg)
74 {
75 const Complex a = reinterpret_cast<const Complex&> (a_arg);
76 return a.real () < 0.0;
77 }
78
79 static F77_INT
80 select_dig (const F77_DBLE_CMPLX& a_arg)
81 {
82 const Complex& a = reinterpret_cast<const Complex&> (a_arg);
83 return (abs (a) < 1.0);
84 }
85
86 static F77_INT
87 select_ana (const F77_CMPLX& a_arg)
88 {
89 const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg);
90 return a.real () < 0.0;
91 }
92
93 static F77_INT
94 select_dig (const F77_CMPLX& a_arg)
95 {
96 const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg);
97 return (abs (a) < 1.0);
98 }
99
100 template <>
101 OCTAVE_API F77_INT
102 schur<Matrix>::init (const Matrix& a, const std::string& ord,
103 bool calc_unitary)
104 {
105 F77_INT a_nr = to_f77_int (a.rows ());
106 F77_INT a_nc = to_f77_int (a.cols ());
107
108 if (a_nr != a_nc)
109 (*current_liboctave_error_handler) ("schur: requires square matrix");
110
111 if (a_nr == 0)
48 { 112 {
49 return (a < 0.0); 113 m_schur_mat.clear ();
114 m_unitary_schur_mat.clear ();
115 return 0;
50 } 116 }
51 117
52 static F77_INT 118 // Workspace requirements may need to be fixed if any of the
53 select_dig (const double& a, const double& b) 119 // following change.
120
121 char jobvs;
122 char sense = 'N';
123 char sort = 'N';
124
125 if (calc_unitary)
126 jobvs = 'V';
127 else
128 jobvs = 'N';
129
130 char ord_char = (ord.empty () ? 'U' : ord[0]);
131
132 if (ord_char == 'A' || ord_char == 'D'
133 || ord_char == 'a' || ord_char == 'd')
134 sort = 'S';
135
136 volatile double_selector selector = nullptr;
137 if (ord_char == 'A' || ord_char == 'a')
138 selector = select_ana;
139 else if (ord_char == 'D' || ord_char == 'd')
140 selector = select_dig;
141
142 F77_INT n = a_nc;
143 F77_INT lwork = 8 * n;
144 F77_INT liwork = 1;
145 F77_INT info;
146 F77_INT sdim;
147 double rconde;
148 double rcondv;
149
150 m_schur_mat = a;
151
152 if (calc_unitary)
153 m_unitary_schur_mat.clear (n, n);
154
155 double *s = m_schur_mat.fortran_vec ();
156 double *q = m_unitary_schur_mat.fortran_vec ();
157
158 Array<double> wr (dim_vector (n, 1));
159 double *pwr = wr.fortran_vec ();
160
161 Array<double> wi (dim_vector (n, 1));
162 double *pwi = wi.fortran_vec ();
163
164 Array<double> work (dim_vector (lwork, 1));
165 double *pwork = work.fortran_vec ();
166
167 // BWORK is not referenced for the non-ordered Schur routine.
168 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
169 Array<F77_INT> bwork (dim_vector (ntmp, 1));
170 F77_INT *pbwork = bwork.fortran_vec ();
171
172 Array<F77_INT> iwork (dim_vector (liwork, 1));
173 F77_INT *piwork = iwork.fortran_vec ();
174
175 F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
176 F77_CONST_CHAR_ARG2 (&sort, 1),
177 selector,
178 F77_CONST_CHAR_ARG2 (&sense, 1),
179 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
180 pwork, lwork, piwork, liwork, pbwork, info
181 F77_CHAR_ARG_LEN (1)
182 F77_CHAR_ARG_LEN (1)
183 F77_CHAR_ARG_LEN (1)));
184
185 return info;
186 }
187
188 template <>
189 OCTAVE_API F77_INT
190 schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord,
191 bool calc_unitary)
192 {
193 F77_INT a_nr = to_f77_int (a.rows ());
194 F77_INT a_nc = to_f77_int (a.cols ());
195
196 if (a_nr != a_nc)
197 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
198
199 if (a_nr == 0)
54 { 200 {
55 return (hypot (a, b) < 1.0); 201 m_schur_mat.clear ();
202 m_unitary_schur_mat.clear ();
203 return 0;
56 } 204 }
57 205
58 static F77_INT 206 // Workspace requirements may need to be fixed if any of the
59 select_ana (const float& a, const float&) 207 // following change.
208
209 char jobvs;
210 char sense = 'N';
211 char sort = 'N';
212
213 if (calc_unitary)
214 jobvs = 'V';
215 else
216 jobvs = 'N';
217
218 char ord_char = (ord.empty () ? 'U' : ord[0]);
219
220 if (ord_char == 'A' || ord_char == 'D'
221 || ord_char == 'a' || ord_char == 'd')
222 sort = 'S';
223
224 volatile float_selector selector = nullptr;
225 if (ord_char == 'A' || ord_char == 'a')
226 selector = select_ana;
227 else if (ord_char == 'D' || ord_char == 'd')
228 selector = select_dig;
229
230 F77_INT n = a_nc;
231 F77_INT lwork = 8 * n;
232 F77_INT liwork = 1;
233 F77_INT info;
234 F77_INT sdim;
235 float rconde;
236 float rcondv;
237
238 m_schur_mat = a;
239
240 if (calc_unitary)
241 m_unitary_schur_mat.clear (n, n);
242
243 float *s = m_schur_mat.fortran_vec ();
244 float *q = m_unitary_schur_mat.fortran_vec ();
245
246 Array<float> wr (dim_vector (n, 1));
247 float *pwr = wr.fortran_vec ();
248
249 Array<float> wi (dim_vector (n, 1));
250 float *pwi = wi.fortran_vec ();
251
252 Array<float> work (dim_vector (lwork, 1));
253 float *pwork = work.fortran_vec ();
254
255 // BWORK is not referenced for the non-ordered Schur routine.
256 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
257 Array<F77_INT> bwork (dim_vector (ntmp, 1));
258 F77_INT *pbwork = bwork.fortran_vec ();
259
260 Array<F77_INT> iwork (dim_vector (liwork, 1));
261 F77_INT *piwork = iwork.fortran_vec ();
262
263 F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
264 F77_CONST_CHAR_ARG2 (&sort, 1),
265 selector,
266 F77_CONST_CHAR_ARG2 (&sense, 1),
267 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
268 pwork, lwork, piwork, liwork, pbwork, info
269 F77_CHAR_ARG_LEN (1)
270 F77_CHAR_ARG_LEN (1)
271 F77_CHAR_ARG_LEN (1)));
272
273 return info;
274 }
275
276 template <>
277 OCTAVE_API F77_INT
278 schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord,
279 bool calc_unitary)
280 {
281 F77_INT a_nr = to_f77_int (a.rows ());
282 F77_INT a_nc = to_f77_int (a.cols ());
283
284 if (a_nr != a_nc)
285 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
286
287 if (a_nr == 0)
60 { 288 {
61 return (a < 0.0); 289 m_schur_mat.clear ();
290 m_unitary_schur_mat.clear ();
291 return 0;
62 } 292 }
63 293
64 static F77_INT 294 // Workspace requirements may need to be fixed if any of the
65 select_dig (const float& a, const float& b) 295 // following change.
296
297 char jobvs;
298 char sense = 'N';
299 char sort = 'N';
300
301 if (calc_unitary)
302 jobvs = 'V';
303 else
304 jobvs = 'N';
305
306 char ord_char = (ord.empty () ? 'U' : ord[0]);
307
308 if (ord_char == 'A' || ord_char == 'D'
309 || ord_char == 'a' || ord_char == 'd')
310 sort = 'S';
311
312 volatile complex_selector selector = nullptr;
313 if (ord_char == 'A' || ord_char == 'a')
314 selector = select_ana;
315 else if (ord_char == 'D' || ord_char == 'd')
316 selector = select_dig;
317
318 F77_INT n = a_nc;
319 F77_INT lwork = 8 * n;
320 F77_INT info;
321 F77_INT sdim;
322 double rconde;
323 double rcondv;
324
325 m_schur_mat = a;
326 if (calc_unitary)
327 m_unitary_schur_mat.clear (n, n);
328
329 Complex *s = m_schur_mat.fortran_vec ();
330 Complex *q = m_unitary_schur_mat.fortran_vec ();
331
332 Array<double> rwork (dim_vector (n, 1));
333 double *prwork = rwork.fortran_vec ();
334
335 Array<Complex> w (dim_vector (n, 1));
336 Complex *pw = w.fortran_vec ();
337
338 Array<Complex> work (dim_vector (lwork, 1));
339 Complex *pwork = work.fortran_vec ();
340
341 // BWORK is not referenced for non-ordered Schur.
342 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
343 Array<F77_INT> bwork (dim_vector (ntmp, 1));
344 F77_INT *pbwork = bwork.fortran_vec ();
345
346 F77_XFCN (zgeesx, ZGEESX,
347 (F77_CONST_CHAR_ARG2 (&jobvs, 1),
348 F77_CONST_CHAR_ARG2 (&sort, 1),
349 selector,
350 F77_CONST_CHAR_ARG2 (&sense, 1),
351 n, F77_DBLE_CMPLX_ARG (s), n, sdim, F77_DBLE_CMPLX_ARG (pw),
352 F77_DBLE_CMPLX_ARG (q), n, rconde, rcondv,
353 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1)
356 F77_CHAR_ARG_LEN (1)));
357
358 return info;
359 }
360
361 template <>
362 OCTAVE_API schur<ComplexMatrix>
363 rsf2csf<ComplexMatrix, Matrix> (const Matrix& s_arg, const Matrix& u_arg)
364 {
365 ComplexMatrix s (s_arg);
366 ComplexMatrix u (u_arg);
367
368 F77_INT n = to_f77_int (s.rows ());
369
370 if (s.columns () != n || u.rows () != n || u.columns () != n)
371 (*current_liboctave_error_handler)
372 ("rsf2csf: inconsistent matrix dimensions");
373
374 if (n > 0)
66 { 375 {
67 return (hypot (a, b) < 1.0); 376 OCTAVE_LOCAL_BUFFER (double, c, n-1);
377 OCTAVE_LOCAL_BUFFER (double, sx, n-1);
378
379 F77_XFCN (zrsf2csf, ZRSF2CSF,
380 (n, F77_DBLE_CMPLX_ARG (s.fortran_vec ()),
381 F77_DBLE_CMPLX_ARG (u.fortran_vec ()), c, sx));
68 } 382 }
69 383
70 // For complex types. 384 return schur<ComplexMatrix> (s, u);
71 385 }
72 static F77_INT 386
73 select_ana (const F77_DBLE_CMPLX& a_arg) 387 template <>
388 OCTAVE_API F77_INT
389 schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
390 const std::string& ord, bool calc_unitary)
391 {
392 F77_INT a_nr = to_f77_int (a.rows ());
393 F77_INT a_nc = to_f77_int (a.cols ());
394
395 if (a_nr != a_nc)
396 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
397
398 if (a_nr == 0)
74 { 399 {
75 const Complex a = reinterpret_cast<const Complex&> (a_arg); 400 m_schur_mat.clear ();
76 return a.real () < 0.0; 401 m_unitary_schur_mat.clear ();
402 return 0;
77 } 403 }
78 404
79 static F77_INT 405 // Workspace requirements may need to be fixed if any of the
80 select_dig (const F77_DBLE_CMPLX& a_arg) 406 // following change.
407
408 char jobvs;
409 char sense = 'N';
410 char sort = 'N';
411
412 if (calc_unitary)
413 jobvs = 'V';
414 else
415 jobvs = 'N';
416
417 char ord_char = (ord.empty () ? 'U' : ord[0]);
418
419 if (ord_char == 'A' || ord_char == 'D'
420 || ord_char == 'a' || ord_char == 'd')
421 sort = 'S';
422
423 volatile float_complex_selector selector = nullptr;
424 if (ord_char == 'A' || ord_char == 'a')
425 selector = select_ana;
426 else if (ord_char == 'D' || ord_char == 'd')
427 selector = select_dig;
428
429 F77_INT n = a_nc;
430 F77_INT lwork = 8 * n;
431 F77_INT info;
432 F77_INT sdim;
433 float rconde;
434 float rcondv;
435
436 m_schur_mat = a;
437 if (calc_unitary)
438 m_unitary_schur_mat.clear (n, n);
439
440 FloatComplex *s = m_schur_mat.fortran_vec ();
441 FloatComplex *q = m_unitary_schur_mat.fortran_vec ();
442
443 Array<float> rwork (dim_vector (n, 1));
444 float *prwork = rwork.fortran_vec ();
445
446 Array<FloatComplex> w (dim_vector (n, 1));
447 FloatComplex *pw = w.fortran_vec ();
448
449 Array<FloatComplex> work (dim_vector (lwork, 1));
450 FloatComplex *pwork = work.fortran_vec ();
451
452 // BWORK is not referenced for non-ordered Schur.
453 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
454 Array<F77_INT> bwork (dim_vector (ntmp, 1));
455 F77_INT *pbwork = bwork.fortran_vec ();
456
457 F77_XFCN (cgeesx, CGEESX,
458 (F77_CONST_CHAR_ARG2 (&jobvs, 1),
459 F77_CONST_CHAR_ARG2 (&sort, 1),
460 selector,
461 F77_CONST_CHAR_ARG2 (&sense, 1),
462 n, F77_CMPLX_ARG (s), n, sdim, F77_CMPLX_ARG (pw),
463 F77_CMPLX_ARG (q), n,
464 rconde, rcondv,
465 F77_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
466 F77_CHAR_ARG_LEN (1)
467 F77_CHAR_ARG_LEN (1)
468 F77_CHAR_ARG_LEN (1)));
469
470 return info;
471 }
472
473 template <>
474 OCTAVE_API schur<FloatComplexMatrix>
475 rsf2csf<FloatComplexMatrix, FloatMatrix> (const FloatMatrix& s_arg,
476 const FloatMatrix& u_arg)
477 {
478 FloatComplexMatrix s (s_arg);
479 FloatComplexMatrix u (u_arg);
480
481 F77_INT n = to_f77_int (s.rows ());
482
483 if (s.columns () != n || u.rows () != n || u.columns () != n)
484 (*current_liboctave_error_handler)
485 ("rsf2csf: inconsistent matrix dimensions");
486
487 if (n > 0)
81 { 488 {
82 const Complex& a = reinterpret_cast<const Complex&> (a_arg); 489 OCTAVE_LOCAL_BUFFER (float, c, n-1);
83 return (abs (a) < 1.0); 490 OCTAVE_LOCAL_BUFFER (float, sx, n-1);
491
492 F77_XFCN (crsf2csf, CRSF2CSF,
493 (n, F77_CMPLX_ARG (s.fortran_vec ()),
494 F77_CMPLX_ARG (u.fortran_vec ()), c, sx));
84 } 495 }
85 496
86 static F77_INT 497 return schur<FloatComplexMatrix> (s, u);
87 select_ana (const F77_CMPLX& a_arg) 498 }
88 { 499
89 const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg); 500 // Instantiations we need.
90 return a.real () < 0.0; 501
91 } 502 template class schur<ComplexMatrix>;
92 503
93 static F77_INT 504 template class schur<FloatComplexMatrix>;
94 select_dig (const F77_CMPLX& a_arg) 505
95 { 506 template class schur<FloatMatrix>;
96 const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg); 507
97 return (abs (a) < 1.0); 508 template class schur<Matrix>;
98 }
99
100 template <>
101 OCTAVE_API F77_INT
102 schur<Matrix>::init (const Matrix& a, const std::string& ord,
103 bool calc_unitary)
104 {
105 F77_INT a_nr = to_f77_int (a.rows ());
106 F77_INT a_nc = to_f77_int (a.cols ());
107
108 if (a_nr != a_nc)
109 (*current_liboctave_error_handler) ("schur: requires square matrix");
110
111 if (a_nr == 0)
112 {
113 m_schur_mat.clear ();
114 m_unitary_schur_mat.clear ();
115 return 0;
116 }
117
118 // Workspace requirements may need to be fixed if any of the
119 // following change.
120
121 char jobvs;
122 char sense = 'N';
123 char sort = 'N';
124
125 if (calc_unitary)
126 jobvs = 'V';
127 else
128 jobvs = 'N';
129
130 char ord_char = (ord.empty () ? 'U' : ord[0]);
131
132 if (ord_char == 'A' || ord_char == 'D'
133 || ord_char == 'a' || ord_char == 'd')
134 sort = 'S';
135
136 volatile double_selector selector = nullptr;
137 if (ord_char == 'A' || ord_char == 'a')
138 selector = select_ana;
139 else if (ord_char == 'D' || ord_char == 'd')
140 selector = select_dig;
141
142 F77_INT n = a_nc;
143 F77_INT lwork = 8 * n;
144 F77_INT liwork = 1;
145 F77_INT info;
146 F77_INT sdim;
147 double rconde;
148 double rcondv;
149
150 m_schur_mat = a;
151
152 if (calc_unitary)
153 m_unitary_schur_mat.clear (n, n);
154
155 double *s = m_schur_mat.fortran_vec ();
156 double *q = m_unitary_schur_mat.fortran_vec ();
157
158 Array<double> wr (dim_vector (n, 1));
159 double *pwr = wr.fortran_vec ();
160
161 Array<double> wi (dim_vector (n, 1));
162 double *pwi = wi.fortran_vec ();
163
164 Array<double> work (dim_vector (lwork, 1));
165 double *pwork = work.fortran_vec ();
166
167 // BWORK is not referenced for the non-ordered Schur routine.
168 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
169 Array<F77_INT> bwork (dim_vector (ntmp, 1));
170 F77_INT *pbwork = bwork.fortran_vec ();
171
172 Array<F77_INT> iwork (dim_vector (liwork, 1));
173 F77_INT *piwork = iwork.fortran_vec ();
174
175 F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
176 F77_CONST_CHAR_ARG2 (&sort, 1),
177 selector,
178 F77_CONST_CHAR_ARG2 (&sense, 1),
179 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
180 pwork, lwork, piwork, liwork, pbwork, info
181 F77_CHAR_ARG_LEN (1)
182 F77_CHAR_ARG_LEN (1)
183 F77_CHAR_ARG_LEN (1)));
184
185 return info;
186 }
187
188 template <>
189 OCTAVE_API F77_INT
190 schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord,
191 bool calc_unitary)
192 {
193 F77_INT a_nr = to_f77_int (a.rows ());
194 F77_INT a_nc = to_f77_int (a.cols ());
195
196 if (a_nr != a_nc)
197 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
198
199 if (a_nr == 0)
200 {
201 m_schur_mat.clear ();
202 m_unitary_schur_mat.clear ();
203 return 0;
204 }
205
206 // Workspace requirements may need to be fixed if any of the
207 // following change.
208
209 char jobvs;
210 char sense = 'N';
211 char sort = 'N';
212
213 if (calc_unitary)
214 jobvs = 'V';
215 else
216 jobvs = 'N';
217
218 char ord_char = (ord.empty () ? 'U' : ord[0]);
219
220 if (ord_char == 'A' || ord_char == 'D'
221 || ord_char == 'a' || ord_char == 'd')
222 sort = 'S';
223
224 volatile float_selector selector = nullptr;
225 if (ord_char == 'A' || ord_char == 'a')
226 selector = select_ana;
227 else if (ord_char == 'D' || ord_char == 'd')
228 selector = select_dig;
229
230 F77_INT n = a_nc;
231 F77_INT lwork = 8 * n;
232 F77_INT liwork = 1;
233 F77_INT info;
234 F77_INT sdim;
235 float rconde;
236 float rcondv;
237
238 m_schur_mat = a;
239
240 if (calc_unitary)
241 m_unitary_schur_mat.clear (n, n);
242
243 float *s = m_schur_mat.fortran_vec ();
244 float *q = m_unitary_schur_mat.fortran_vec ();
245
246 Array<float> wr (dim_vector (n, 1));
247 float *pwr = wr.fortran_vec ();
248
249 Array<float> wi (dim_vector (n, 1));
250 float *pwi = wi.fortran_vec ();
251
252 Array<float> work (dim_vector (lwork, 1));
253 float *pwork = work.fortran_vec ();
254
255 // BWORK is not referenced for the non-ordered Schur routine.
256 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
257 Array<F77_INT> bwork (dim_vector (ntmp, 1));
258 F77_INT *pbwork = bwork.fortran_vec ();
259
260 Array<F77_INT> iwork (dim_vector (liwork, 1));
261 F77_INT *piwork = iwork.fortran_vec ();
262
263 F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
264 F77_CONST_CHAR_ARG2 (&sort, 1),
265 selector,
266 F77_CONST_CHAR_ARG2 (&sense, 1),
267 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
268 pwork, lwork, piwork, liwork, pbwork, info
269 F77_CHAR_ARG_LEN (1)
270 F77_CHAR_ARG_LEN (1)
271 F77_CHAR_ARG_LEN (1)));
272
273 return info;
274 }
275
276 template <>
277 OCTAVE_API F77_INT
278 schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord,
279 bool calc_unitary)
280 {
281 F77_INT a_nr = to_f77_int (a.rows ());
282 F77_INT a_nc = to_f77_int (a.cols ());
283
284 if (a_nr != a_nc)
285 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
286
287 if (a_nr == 0)
288 {
289 m_schur_mat.clear ();
290 m_unitary_schur_mat.clear ();
291 return 0;
292 }
293
294 // Workspace requirements may need to be fixed if any of the
295 // following change.
296
297 char jobvs;
298 char sense = 'N';
299 char sort = 'N';
300
301 if (calc_unitary)
302 jobvs = 'V';
303 else
304 jobvs = 'N';
305
306 char ord_char = (ord.empty () ? 'U' : ord[0]);
307
308 if (ord_char == 'A' || ord_char == 'D'
309 || ord_char == 'a' || ord_char == 'd')
310 sort = 'S';
311
312 volatile complex_selector selector = nullptr;
313 if (ord_char == 'A' || ord_char == 'a')
314 selector = select_ana;
315 else if (ord_char == 'D' || ord_char == 'd')
316 selector = select_dig;
317
318 F77_INT n = a_nc;
319 F77_INT lwork = 8 * n;
320 F77_INT info;
321 F77_INT sdim;
322 double rconde;
323 double rcondv;
324
325 m_schur_mat = a;
326 if (calc_unitary)
327 m_unitary_schur_mat.clear (n, n);
328
329 Complex *s = m_schur_mat.fortran_vec ();
330 Complex *q = m_unitary_schur_mat.fortran_vec ();
331
332 Array<double> rwork (dim_vector (n, 1));
333 double *prwork = rwork.fortran_vec ();
334
335 Array<Complex> w (dim_vector (n, 1));
336 Complex *pw = w.fortran_vec ();
337
338 Array<Complex> work (dim_vector (lwork, 1));
339 Complex *pwork = work.fortran_vec ();
340
341 // BWORK is not referenced for non-ordered Schur.
342 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
343 Array<F77_INT> bwork (dim_vector (ntmp, 1));
344 F77_INT *pbwork = bwork.fortran_vec ();
345
346 F77_XFCN (zgeesx, ZGEESX,
347 (F77_CONST_CHAR_ARG2 (&jobvs, 1),
348 F77_CONST_CHAR_ARG2 (&sort, 1),
349 selector,
350 F77_CONST_CHAR_ARG2 (&sense, 1),
351 n, F77_DBLE_CMPLX_ARG (s), n, sdim, F77_DBLE_CMPLX_ARG (pw),
352 F77_DBLE_CMPLX_ARG (q), n, rconde, rcondv,
353 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1)
356 F77_CHAR_ARG_LEN (1)));
357
358 return info;
359 }
360
361 template <>
362 OCTAVE_API schur<ComplexMatrix>
363 rsf2csf<ComplexMatrix, Matrix> (const Matrix& s_arg, const Matrix& u_arg)
364 {
365 ComplexMatrix s (s_arg);
366 ComplexMatrix u (u_arg);
367
368 F77_INT n = to_f77_int (s.rows ());
369
370 if (s.columns () != n || u.rows () != n || u.columns () != n)
371 (*current_liboctave_error_handler)
372 ("rsf2csf: inconsistent matrix dimensions");
373
374 if (n > 0)
375 {
376 OCTAVE_LOCAL_BUFFER (double, c, n-1);
377 OCTAVE_LOCAL_BUFFER (double, sx, n-1);
378
379 F77_XFCN (zrsf2csf, ZRSF2CSF,
380 (n, F77_DBLE_CMPLX_ARG (s.fortran_vec ()),
381 F77_DBLE_CMPLX_ARG (u.fortran_vec ()), c, sx));
382 }
383
384 return schur<ComplexMatrix> (s, u);
385 }
386
387 template <>
388 OCTAVE_API F77_INT
389 schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
390 const std::string& ord, bool calc_unitary)
391 {
392 F77_INT a_nr = to_f77_int (a.rows ());
393 F77_INT a_nc = to_f77_int (a.cols ());
394
395 if (a_nr != a_nc)
396 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
397
398 if (a_nr == 0)
399 {
400 m_schur_mat.clear ();
401 m_unitary_schur_mat.clear ();
402 return 0;
403 }
404
405 // Workspace requirements may need to be fixed if any of the
406 // following change.
407
408 char jobvs;
409 char sense = 'N';
410 char sort = 'N';
411
412 if (calc_unitary)
413 jobvs = 'V';
414 else
415 jobvs = 'N';
416
417 char ord_char = (ord.empty () ? 'U' : ord[0]);
418
419 if (ord_char == 'A' || ord_char == 'D'
420 || ord_char == 'a' || ord_char == 'd')
421 sort = 'S';
422
423 volatile float_complex_selector selector = nullptr;
424 if (ord_char == 'A' || ord_char == 'a')
425 selector = select_ana;
426 else if (ord_char == 'D' || ord_char == 'd')
427 selector = select_dig;
428
429 F77_INT n = a_nc;
430 F77_INT lwork = 8 * n;
431 F77_INT info;
432 F77_INT sdim;
433 float rconde;
434 float rcondv;
435
436 m_schur_mat = a;
437 if (calc_unitary)
438 m_unitary_schur_mat.clear (n, n);
439
440 FloatComplex *s = m_schur_mat.fortran_vec ();
441 FloatComplex *q = m_unitary_schur_mat.fortran_vec ();
442
443 Array<float> rwork (dim_vector (n, 1));
444 float *prwork = rwork.fortran_vec ();
445
446 Array<FloatComplex> w (dim_vector (n, 1));
447 FloatComplex *pw = w.fortran_vec ();
448
449 Array<FloatComplex> work (dim_vector (lwork, 1));
450 FloatComplex *pwork = work.fortran_vec ();
451
452 // BWORK is not referenced for non-ordered Schur.
453 F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
454 Array<F77_INT> bwork (dim_vector (ntmp, 1));
455 F77_INT *pbwork = bwork.fortran_vec ();
456
457 F77_XFCN (cgeesx, CGEESX,
458 (F77_CONST_CHAR_ARG2 (&jobvs, 1),
459 F77_CONST_CHAR_ARG2 (&sort, 1),
460 selector,
461 F77_CONST_CHAR_ARG2 (&sense, 1),
462 n, F77_CMPLX_ARG (s), n, sdim, F77_CMPLX_ARG (pw),
463 F77_CMPLX_ARG (q), n,
464 rconde, rcondv,
465 F77_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
466 F77_CHAR_ARG_LEN (1)
467 F77_CHAR_ARG_LEN (1)
468 F77_CHAR_ARG_LEN (1)));
469
470 return info;
471 }
472
473 template <>
474 OCTAVE_API schur<FloatComplexMatrix>
475 rsf2csf<FloatComplexMatrix, FloatMatrix> (const FloatMatrix& s_arg,
476 const FloatMatrix& u_arg)
477 {
478 FloatComplexMatrix s (s_arg);
479 FloatComplexMatrix u (u_arg);
480
481 F77_INT n = to_f77_int (s.rows ());
482
483 if (s.columns () != n || u.rows () != n || u.columns () != n)
484 (*current_liboctave_error_handler)
485 ("rsf2csf: inconsistent matrix dimensions");
486
487 if (n > 0)
488 {
489 OCTAVE_LOCAL_BUFFER (float, c, n-1);
490 OCTAVE_LOCAL_BUFFER (float, sx, n-1);
491
492 F77_XFCN (crsf2csf, CRSF2CSF,
493 (n, F77_CMPLX_ARG (s.fortran_vec ()),
494 F77_CMPLX_ARG (u.fortran_vec ()), c, sx));
495 }
496
497 return schur<FloatComplexMatrix> (s, u);
498 }
499
500 // Instantiations we need.
501
502 template class schur<ComplexMatrix>;
503
504 template class schur<FloatComplexMatrix>;
505
506 template class schur<FloatMatrix>;
507
508 template class schur<Matrix>;
509 509
510 OCTAVE_END_NAMESPACE(math) 510 OCTAVE_END_NAMESPACE(math)
511 OCTAVE_END_NAMESPACE(octave) 511 OCTAVE_END_NAMESPACE(octave)