Mercurial > octave
annotate libinterp/corefcn/qz.cc @ 21200:fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
* builtins.h, Cell.cc, __contourc__.cc, __dispatch__.cc, __dsearchn__.cc,
__ichol__.cc, __ilu__.cc, __lin_interpn__.cc, __pchip_deriv__.cc, __qp__.cc,
balance.cc, besselj.cc, betainc.cc, bitfcns.cc, bsxfun.cc,
c-file-ptr-stream.cc, c-file-ptr-stream.h, cellfun.cc, colloc.cc,
comment-list.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc, debug.cc,
defaults.cc, defaults.in.h, defun-dld.h, defun.cc, defun.h, det.cc, dirfns.cc,
display.cc, dlmread.cc, dot.cc, dynamic-ld.cc, eig.cc, ellipj.cc, error.cc,
errwarn.cc, event-queue.cc, fft.cc, fft2.cc, fftn.cc, file-io.cc, filter.cc,
find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc,
givens.cc, gl-render.cc, gl2ps-print.cc, graphics.cc, graphics.in.h, gripes.cc,
hash.cc, help.cc, hess.cc, hex2num.cc, input.cc, inv.cc, jit-ir.cc,
jit-typeinfo.cc, jit-util.cc, jit-util.h, kron.cc, load-path.cc, load-save.cc,
lookup.cc, ls-ascii-helper.cc, ls-hdf5.cc, ls-mat-ascii.cc, ls-mat4.cc,
ls-mat5.cc, ls-oct-binary.cc, ls-oct-text.cc, ls-oct-text.h, ls-utils.cc,
ls-utils.h, lsode.cc, lu.cc, luinc.cc, mappers.cc, matrix_type.cc, max.cc,
mex.h, mexproto.h, mgorth.cc, nproc.cc, oct-errno.in.cc, oct-fstrm.cc,
oct-hdf5-types.cc, oct-hdf5.h, oct-hist.cc, oct-iostrm.cc, oct-lvalue.cc,
oct-map.cc, oct-prcstrm.cc, oct-procbuf.cc, oct-stream.cc, oct-strstrm.cc,
octave-link.cc, ordschur.cc, pager.cc, pinv.cc, pr-output.cc, procstream.cc,
profiler.cc, psi.cc, pt-jit.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc,
regexp.cc, schur.cc, sighandlers.cc, sparse-xdiv.cc, sparse-xpow.cc, sparse.cc,
spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc,
sylvester.cc, symtab.cc, syscalls.cc, sysdep.cc, sysdep.h, time.cc, toplev.cc,
tril.cc, tsearch.cc, txt-eng-ft.cc, txt-eng.cc, typecast.cc, urlwrite.cc,
utils.cc, variables.cc, xdiv.cc, xnorm.cc, xpow.cc, zfstream.cc,
__delaunayn__.cc, __eigs__.cc, __fltk_uigetfile__.cc, __glpk__.cc,
__init_fltk__.cc, __init_gnuplot__.cc, __magick_read__.cc, __osmesa_print__.cc,
__voronoi__.cc, amd.cc, audiodevinfo.cc, audioread.cc, ccolamd.cc, chol.cc,
colamd.cc, convhulln.cc, dmperm.cc, fftw.cc, oct-qhull.h, qr.cc, symbfact.cc,
symrcm.cc, oct-conf.in.cc, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc,
ov-base-scalar.cc, ov-base-sparse.cc, ov-base.cc, ov-bool-mat.cc,
ov-bool-sparse.cc, ov-bool.cc, ov-builtin.cc, ov-cell.cc, ov-ch-mat.cc,
ov-class.cc, ov-classdef.cc, ov-colon.cc, ov-complex.cc, ov-cs-list.cc,
ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-fcn-handle.cc,
ov-fcn-inline.cc, ov-fcn.cc, 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-int16.cc,
ov-int32.cc, ov-int64.cc, ov-int8.cc, ov-java.cc, ov-lazy-idx.cc,
ov-mex-fcn.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-perm.cc, ov-range.cc,
ov-re-diag.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc,
ov-struct.cc, ov-typeinfo.cc, ov-uint16.cc, ov-uint32.cc, ov-uint64.cc,
ov-uint8.cc, ov-usr-fcn.cc, ov.cc, ovl.cc, octave.cc, op-b-b.cc, op-b-bm.cc,
op-b-sbm.cc, op-bm-b.cc, op-bm-bm.cc, op-bm-sbm.cc, op-cdm-cdm.cc, op-cell.cc,
op-chm.cc, op-class.cc, op-cm-cm.cc, op-cm-cs.cc, op-cm-m.cc, op-cm-s.cc,
op-cm-scm.cc, op-cm-sm.cc, op-cs-cm.cc, op-cs-cs.cc, op-cs-m.cc, op-cs-s.cc,
op-cs-scm.cc, op-cs-sm.cc, op-dm-dm.cc, op-dm-scm.cc, op-dm-sm.cc,
op-dm-template.cc, op-dms-template.cc, op-double-conv.cc, op-fcdm-fcdm.cc,
op-fcdm-fdm.cc, op-fcm-fcm.cc, op-fcm-fcs.cc, op-fcm-fm.cc, op-fcm-fs.cc,
op-fcn.cc, op-fcs-fcm.cc, op-fcs-fcs.cc, op-fcs-fm.cc, op-fcs-fs.cc,
op-fdm-fdm.cc, op-float-conv.cc, op-fm-fcm.cc, op-fm-fcs.cc, op-fm-fm.cc,
op-fm-fs.cc, op-fs-fcm.cc, op-fs-fcs.cc, op-fs-fm.cc, op-fs-fs.cc,
op-i16-i16.cc, op-i32-i32.cc, op-i64-i64.cc, op-i8-i8.cc, op-int-concat.cc,
op-int-conv.cc, op-m-cm.cc, op-m-cs.cc, op-m-m.cc, op-m-s.cc, op-m-scm.cc,
op-m-sm.cc, op-pm-pm.cc, op-pm-scm.cc, op-pm-sm.cc, op-pm-template.cc,
op-range.cc, op-s-cm.cc, op-s-cs.cc, op-s-m.cc, op-s-s.cc, op-s-scm.cc,
op-s-sm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-sbm-sbm.cc, op-scm-cm.cc,
op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc, op-scm-sm.cc,
op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc, op-sm-sm.cc,
op-str-m.cc, op-str-s.cc, op-str-str.cc, op-struct.cc, op-ui16-ui16.cc,
op-ui32-ui32.cc, op-ui64-ui64.cc, op-ui8-ui8.cc, pt-arg-list.cc,
pt-array-list.cc, pt-assign.cc, pt-binop.cc, pt-bp.cc, pt-cbinop.cc,
pt-cell.cc, pt-check.cc, pt-classdef.cc, pt-cmd.cc, pt-colon.cc, pt-colon.h,
pt-const.cc, pt-decl.cc, pt-eval.cc, pt-except.cc, pt-exp.cc, pt-fcn-handle.cc,
pt-funcall.cc, pt-id.cc, pt-idx.cc, pt-jump.cc, pt-loop.cc, pt-mat.cc,
pt-misc.cc, pt-pr-code.cc, pt-select.cc, pt-stmt.cc, pt-unop.cc, pt.cc,
token.cc, Array-jit.cc, Array-os.cc, Array-sym.cc, Array-tc.cc, version.cc:
Indent #ifdef blocks in libinterp.
| author | Rik <rik@octave.org> |
|---|---|
| date | Fri, 05 Feb 2016 16:29:08 -0800 |
| parents | 3d0d84305600 |
| children | eb1524b07fe3 |
| rev | line source |
|---|---|
| 3183 | 1 /* |
| 2 | |
|
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19403
diff
changeset
|
3 Copyright (C) 1998-2015 A. S. Hodel |
| 3183 | 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 | |
| 7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
| 10 option) any later version. | |
| 3183 | 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 | |
| 7016 | 18 along with Octave; see the file COPYING. If not, see |
| 19 <http://www.gnu.org/licenses/>. | |
| 3183 | 20 |
| 21 */ | |
| 22 | |
| 23 // Generalized eigenvalue balancing via LAPACK | |
| 3911 | 24 |
| 25 // Author: A. S. Hodel <scotte@eng.auburn.edu> | |
| 3183 | 26 |
| 27 #undef DEBUG | |
| 28 #undef DEBUG_SORT | |
| 29 #undef DEBUG_EIG | |
| 30 | |
|
9786
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
31 #ifdef HAVE_CONFIG_H |
|
21200
fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
Rik <rik@octave.org>
parents:
21110
diff
changeset
|
32 # include <config.h> |
|
9786
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
33 #endif |
| 3183 | 34 |
| 35 #include <cfloat> | |
| 4051 | 36 |
| 3523 | 37 #include <iostream> |
| 4051 | 38 #include <iomanip> |
| 3183 | 39 |
| 40 #include "CmplxQRP.h" | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
41 #include "CmplxQR.h" |
| 3183 | 42 #include "dbleQR.h" |
| 4153 | 43 #include "f77-fcn.h" |
| 7231 | 44 #include "lo-math.h" |
| 4153 | 45 #include "quit.h" |
| 46 | |
|
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
47 #include "defun.h" |
| 3183 | 48 #include "error.h" |
|
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
49 #include "errwarn.h" |
|
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
20892
diff
changeset
|
50 #include "ovl.h" |
| 3183 | 51 #include "oct-map.h" |
| 52 #include "ov.h" | |
| 53 #include "pager.h" | |
| 3185 | 54 #if defined (DEBUG) || defined (DEBUG_SORT) |
|
21200
fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
Rik <rik@octave.org>
parents:
21110
diff
changeset
|
55 # include "pr-output.h" |
| 3183 | 56 #endif |
| 57 #include "symtab.h" | |
| 58 #include "utils.h" | |
| 59 #include "variables.h" | |
| 60 | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
61 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE, |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
62 const double& ALPHA, |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
63 const double& BETA, const double& S, |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
64 const double& P); |
| 3183 | 65 |
| 66 extern "C" | |
| 67 { | |
| 4552 | 68 F77_RET_T |
| 69 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, | |
| 10311 | 70 const octave_idx_type& N, double* A, |
| 71 const octave_idx_type& LDA, double* B, | |
| 72 const octave_idx_type& LDB, octave_idx_type& ILO, | |
| 73 octave_idx_type& IHI, double* LSCALE, | |
| 74 double* RSCALE, double* WORK, | |
| 75 octave_idx_type& INFO | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
76 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 77 |
|
19861
19755f4fc851
maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
78 F77_RET_T |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
79 F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL, |
| 10311 | 80 const octave_idx_type& N, Complex* A, |
| 81 const octave_idx_type& LDA, Complex* B, | |
| 82 const octave_idx_type& LDB, octave_idx_type& ILO, | |
| 83 octave_idx_type& IHI, double* LSCALE, | |
| 84 double* RSCALE, double* WORK, | |
| 85 octave_idx_type& INFO | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
86 F77_CHAR_ARG_LEN_DECL); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
87 |
| 4552 | 88 F77_RET_T |
| 89 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
90 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 91 const octave_idx_type& N, |
| 92 const octave_idx_type& ILO, | |
| 93 const octave_idx_type& IHI, | |
| 94 const double* LSCALE, const double* RSCALE, | |
| 95 octave_idx_type& M, double* V, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
96 const octave_idx_type& LDV, octave_idx_type& INFO |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
97 F77_CHAR_ARG_LEN_DECL |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
98 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 99 |
|
19861
19755f4fc851
maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
100 F77_RET_T |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
101 F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
102 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 103 const octave_idx_type& N, |
| 104 const octave_idx_type& ILO, | |
| 105 const octave_idx_type& IHI, | |
| 106 const double* LSCALE, const double* RSCALE, | |
| 107 octave_idx_type& M, Complex* V, | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
108 const octave_idx_type& LDV, octave_idx_type& INFO |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
109 F77_CHAR_ARG_LEN_DECL |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
110 F77_CHAR_ARG_LEN_DECL); |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
111 |
| 4552 | 112 F77_RET_T |
| 113 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
114 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 115 const octave_idx_type& N, |
| 116 const octave_idx_type& ILO, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
117 const octave_idx_type& IHI, double* A, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
118 const octave_idx_type& LDA, double* B, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
119 const octave_idx_type& LDB, double* Q, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
120 const octave_idx_type& LDQ, double* Z, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
121 const octave_idx_type& LDZ, octave_idx_type& INFO |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
122 F77_CHAR_ARG_LEN_DECL |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
123 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 124 |
|
19861
19755f4fc851
maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
125 F77_RET_T |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
126 F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
127 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 128 const octave_idx_type& N, |
| 129 const octave_idx_type& ILO, | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
130 const octave_idx_type& IHI, Complex* A, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
131 const octave_idx_type& LDA, Complex* B, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
132 const octave_idx_type& LDB, Complex* Q, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
133 const octave_idx_type& LDQ, Complex* Z, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
134 const octave_idx_type& LDZ, octave_idx_type& INFO |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
135 F77_CHAR_ARG_LEN_DECL |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
136 F77_CHAR_ARG_LEN_DECL); |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
137 |
| 4552 | 138 F77_RET_T |
| 139 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
140 F77_CONST_CHAR_ARG_DECL, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
141 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 142 const octave_idx_type& N, |
| 143 const octave_idx_type& ILO, | |
| 144 const octave_idx_type& IHI, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
145 double* A, const octave_idx_type& LDA, double* B, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
146 const octave_idx_type& LDB, double* ALPHAR, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
147 double* ALPHAI, double* BETA, double* Q, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
148 const octave_idx_type& LDQ, double* Z, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
149 const octave_idx_type& LDZ, double* WORK, |
| 10311 | 150 const octave_idx_type& LWORK, |
| 151 octave_idx_type& INFO | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
152 F77_CHAR_ARG_LEN_DECL |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
153 F77_CHAR_ARG_LEN_DECL |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
154 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 155 |
|
19861
19755f4fc851
maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
156 F77_RET_T |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
157 F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
158 F77_CONST_CHAR_ARG_DECL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
159 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 160 const octave_idx_type& N, |
| 161 const octave_idx_type& ILO, | |
| 162 const octave_idx_type& IHI, | |
| 163 Complex* A, const octave_idx_type& LDA, | |
| 164 Complex* B, const octave_idx_type& LDB, | |
| 165 Complex* ALPHA, Complex* BETA, Complex* CQ, | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
166 const octave_idx_type& LDQ, |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
167 Complex* CZ, const octave_idx_type& LDZ, |
| 10311 | 168 Complex* WORK, const octave_idx_type& LWORK, |
| 169 double* RWORK, octave_idx_type& INFO | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
170 F77_CHAR_ARG_LEN_DECL |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
171 F77_CHAR_ARG_LEN_DECL |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
172 F77_CHAR_ARG_LEN_DECL); |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
173 |
| 4552 | 174 F77_RET_T |
| 10311 | 175 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, |
| 176 const double* B, const octave_idx_type& LDB, | |
| 177 const double& SAFMIN, double& SCALE1, | |
| 178 double& SCALE2, double& WR1, double& WR2, | |
| 179 double& WI); | |
| 3183 | 180 |
| 181 // Van Dooren's code (netlib.org: toms/590) for reordering | |
| 182 // GEP. Only processes Z, not Q. | |
| 4552 | 183 F77_RET_T |
| 10311 | 184 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX, |
| 185 const octave_idx_type& N, double* A, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
186 double* B, double* Z, sort_function, |
| 10311 | 187 const double& EPS, octave_idx_type& NDIM, |
| 188 octave_idx_type& FAIL, octave_idx_type* IND); | |
| 3183 | 189 |
| 10311 | 190 // Documentation for DTGEVC incorrectly states that VR, VL are |
| 3183 | 191 // complex*16; they are declared in DTGEVC as double precision |
| 10311 | 192 // (probably a cut and paste problem fro ZTGEVC). |
| 4552 | 193 F77_RET_T |
| 194 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
195 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 196 octave_idx_type* SELECT, |
| 197 const octave_idx_type& N, double* A, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
198 const octave_idx_type& LDA, double* B, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
199 const octave_idx_type& LDB, double* VL, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
200 const octave_idx_type& LDVL, double* VR, |
| 10311 | 201 const octave_idx_type& LDVR, |
| 202 const octave_idx_type& MM, octave_idx_type& M, | |
| 203 double* WORK, octave_idx_type& INFO | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
204 F77_CHAR_ARG_LEN_DECL |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
205 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 206 |
|
19861
19755f4fc851
maint: Cleanup C++ code to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
207 F77_RET_T |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
208 F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
209 F77_CONST_CHAR_ARG_DECL, |
| 10311 | 210 octave_idx_type* SELECT, |
| 211 const octave_idx_type& N, const Complex* A, | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
212 const octave_idx_type& LDA,const Complex* B, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
213 const octave_idx_type& LDB, Complex* xVL, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
214 const octave_idx_type& LDVL, Complex* xVR, |
| 10311 | 215 const octave_idx_type& LDVR, |
| 216 const octave_idx_type& MM, octave_idx_type& M, | |
| 217 Complex* CWORK, double* RWORK, | |
| 218 octave_idx_type& INFO | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
219 F77_CHAR_ARG_LEN_DECL |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
220 F77_CHAR_ARG_LEN_DECL); |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
221 |
| 4552 | 222 F77_RET_T |
| 223 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
224 double& retval |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
225 F77_CHAR_ARG_LEN_DECL); |
| 3185 | 226 |
| 4552 | 227 F77_RET_T |
| 228 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL, | |
| 10311 | 229 const octave_idx_type&, |
| 230 const octave_idx_type&, const double*, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
231 const octave_idx_type&, double*, double& |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
232 F77_CHAR_ARG_LEN_DECL); |
| 3183 | 233 } |
| 234 | |
| 235 // fcrhp, fin, fout, folhp: | |
| 236 // routines for ordering of generalized eigenvalues | |
| 237 // return 1 if test is passed, 0 otherwise | |
| 238 // fin: |lambda| < 1 | |
| 239 // fout: |lambda| >= 1 | |
| 240 // fcrhp: real(lambda) >= 0 | |
| 241 // folhp: real(lambda) < 0 | |
| 242 | |
| 5275 | 243 static octave_idx_type |
| 244 fcrhp (const octave_idx_type& lsize, const double& alpha, | |
| 3185 | 245 const double& beta, const double& s, const double&) |
| 3183 | 246 { |
| 3185 | 247 if (lsize == 1) |
| 10311 | 248 return (alpha * beta >= 0 ? 1 : -1); |
| 3185 | 249 else |
| 3183 | 250 return (s >= 0 ? 1 : -1); |
| 251 } | |
| 3185 | 252 |
| 5275 | 253 static octave_idx_type |
| 254 fin (const octave_idx_type& lsize, const double& alpha, | |
| 3185 | 255 const double& beta, const double&, const double& p) |
| 3183 | 256 { |
| 5275 | 257 octave_idx_type retval; |
| 3183 | 258 |
| 3185 | 259 if (lsize == 1) |
| 260 retval = (fabs (alpha) < fabs (beta) ? 1 : -1); | |
| 261 else | |
| 262 retval = (fabs (p) < 1 ? 1 : -1); | |
| 263 | |
| 264 #ifdef DEBUG | |
| 3538 | 265 std::cout << "qz: fin: retval=" << retval << std::endl; |
| 3185 | 266 #endif |
| 267 | |
| 3183 | 268 return retval; |
| 269 } | |
| 3185 | 270 |
| 5275 | 271 static octave_idx_type |
| 272 folhp (const octave_idx_type& lsize, const double& alpha, | |
| 3185 | 273 const double& beta, const double& s, const double&) |
| 3183 | 274 { |
| 3185 | 275 if (lsize == 1) |
| 10311 | 276 return (alpha * beta < 0 ? 1 : -1); |
| 3185 | 277 else |
| 3183 | 278 return (s < 0 ? 1 : -1); |
| 279 } | |
| 3185 | 280 |
| 5275 | 281 static octave_idx_type |
| 282 fout (const octave_idx_type& lsize, const double& alpha, | |
| 3185 | 283 const double& beta, const double&, const double& p) |
| 3183 | 284 { |
| 3185 | 285 if (lsize == 1) |
| 286 return (fabs (alpha) >= fabs (beta) ? 1 : -1); | |
| 287 else | |
| 288 return (fabs (p) >= 1 ? 1 : -1); | |
| 3183 | 289 } |
| 290 | |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
291 |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
292 //FIXME: Matlab does not produce lambda as the first output argument. |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
293 // Compatibility problem? |
|
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14876
diff
changeset
|
294 DEFUN (qz, args, nargout, |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
295 "-*- texinfo -*-\n\ |
|
20853
1142cf6abc0d
2015 Code Sprint: remove class of function from docstring for all C++ files.
Rik <rik@octave.org>
parents:
20831
diff
changeset
|
296 @deftypefn {} {@var{lambda} =} qz (@var{A}, @var{B})\n\ |
|
1142cf6abc0d
2015 Code Sprint: remove class of function from docstring for all C++ files.
Rik <rik@octave.org>
parents:
20831
diff
changeset
|
297 @deftypefnx {} {@var{lambda} =} qz (@var{A}, @var{B}, @var{opt})\n\ |
|
11593
1577c6f80926
Use non-breaking spaces between certain adjectives and their nouns in docstrings.
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
298 QZ@tie{}decomposition of the generalized eigenvalue problem\n\ |
|
20172
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
299 (@math{A x = s B x}).\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
300 \n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
301 There are three ways to call this function:\n\ |
| 3372 | 302 @enumerate\n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
303 @item @code{@var{lambda} = qz (@var{A}, @var{B})}\n\ |
| 3185 | 304 \n\ |
| 5016 | 305 Computes the generalized eigenvalues\n\ |
| 306 @tex\n\ | |
| 307 $\\lambda$\n\ | |
| 308 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
309 @ifnottex\n\ |
| 5016 | 310 @var{lambda}\n\ |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
311 @end ifnottex\n\ |
| 5016 | 312 of @math{(A - s B)}.\n\ |
| 10840 | 313 \n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
314 @item @code{[AA, BB, Q, Z, V, W, @var{lambda}] = qz (@var{A}, @var{B})}\n\ |
| 3185 | 315 \n\ |
|
20172
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
316 Computes QZ@tie{}decomposition, generalized eigenvectors, and generalized\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
317 eigenvalues of @math{(A - s B)}\n\ |
| 5016 | 318 @tex\n\ |
| 319 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\ | |
| 320 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\ | |
| 321 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\ | |
| 322 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
323 @ifnottex\n\ |
| 10840 | 324 \n\ |
| 3372 | 325 @example\n\ |
| 326 @group\n\ | |
| 5481 | 327 \n\ |
|
14360
97883071e8e4
doc: Correct off-by-1 spacings in all .cc docstrings
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
328 A * V = B * V * diag (@var{lambda})\n\ |
|
97883071e8e4
doc: Correct off-by-1 spacings in all .cc docstrings
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
329 W' * A = diag (@var{lambda}) * W' * B\n\ |
|
97883071e8e4
doc: Correct off-by-1 spacings in all .cc docstrings
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
330 AA = Q * A * Z, BB = Q * B * Z\n\ |
| 5481 | 331 \n\ |
| 3372 | 332 @end group\n\ |
| 333 @end example\n\ | |
| 10840 | 334 \n\ |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8286
diff
changeset
|
335 @end ifnottex\n\ |
| 5016 | 336 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\ |
| 3185 | 337 \n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
338 @item @code{[AA,BB,Z@{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}\n\ |
| 3185 | 339 \n\ |
|
20172
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
340 As in form [2], but allows ordering of generalized eigenpairs for, e.g.,\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
341 solution of discrete time algebraic Riccati equations. Form 3 is not\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
342 available for complex matrices, and does not compute the generalized\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
343 eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.\n\ |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10350
diff
changeset
|
344 \n\ |
| 3372 | 345 @table @var\n\ |
| 346 @item opt\n\ | |
|
20172
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
347 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of\n\ |
|
4f45eaf83908
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19861
diff
changeset
|
348 the revised pencil contains all eigenvalues that satisfy:\n\ |
|
14360
97883071e8e4
doc: Correct off-by-1 spacings in all .cc docstrings
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
349 \n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
350 @table @asis\n\ |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
17097
diff
changeset
|
351 @item @qcode{\"N\"}\n\ |
|
12642
f96b9b9f141b
doc: Periodic grammarcheck and spellcheck of documentation.
Rik <octave@nomad.inbox5.com>
parents:
11593
diff
changeset
|
352 = unordered (default)\n\ |
| 3183 | 353 \n\ |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
17097
diff
changeset
|
354 @item @qcode{\"S\"}\n\ |
|
12642
f96b9b9f141b
doc: Periodic grammarcheck and spellcheck of documentation.
Rik <octave@nomad.inbox5.com>
parents:
11593
diff
changeset
|
355 = small: leading block has all |lambda| @leq{} 1\n\ |
| 3185 | 356 \n\ |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
17097
diff
changeset
|
357 @item @qcode{\"B\"}\n\ |
|
12642
f96b9b9f141b
doc: Periodic grammarcheck and spellcheck of documentation.
Rik <octave@nomad.inbox5.com>
parents:
11593
diff
changeset
|
358 = big: leading block has all |lambda| @geq{} 1\n\ |
| 3372 | 359 \n\ |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
17097
diff
changeset
|
360 @item @qcode{\"-\"}\n\ |
| 5481 | 361 = negative real part: leading block has all eigenvalues\n\ |
| 362 in the open left half-plane\n\ | |
| 3372 | 363 \n\ |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
17097
diff
changeset
|
364 @item @qcode{\"+\"}\n\ |
|
9209
923c7cb7f13f
Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction.
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
365 = non-negative real part: leading block has all eigenvalues\n\ |
| 5481 | 366 in the closed right half-plane\n\ |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8927
diff
changeset
|
367 @end table\n\ |
| 3372 | 368 @end table\n\ |
| 369 @end enumerate\n\ | |
| 3183 | 370 \n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
371 Note: @code{qz} performs permutation balancing, but not scaling\n\ |
|
17097
e7a059a9a644
doc: Use XREF as anchor prefix in documentation for clearer results in Info viewer.
Rik <rik@octave.org>
parents:
16920
diff
changeset
|
372 (@pxref{XREFbalance}). The order of output arguments was selected for\n\ |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
373 compatibility with @sc{matlab}.\n\ |
|
16920
53eaa83e4181
doc: Add seealso links between various factorization forms.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
374 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd}\n\ |
| 3372 | 375 @end deftypefn") |
| 3183 | 376 { |
| 377 int nargin = args.length (); | |
| 378 | |
| 3185 | 379 #ifdef DEBUG |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
380 std::cout << "qz: nargin = " << nargin |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
381 << ", nargout = " << nargout << std::endl; |
| 3185 | 382 #endif |
| 3183 | 383 |
| 3185 | 384 if (nargin < 2 || nargin > 3 || nargout > 7) |
|
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
385 print_usage (); |
|
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
386 |
|
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
387 if (nargin == 3 && (nargout < 3 || nargout > 4)) |
|
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20700
diff
changeset
|
388 error ("qz: invalid number of output arguments for form [3] call"); |
| 3183 | 389 |
| 3185 | 390 #ifdef DEBUG |
| 3538 | 391 std::cout << "qz: determine ordering option" << std::endl; |
| 3185 | 392 #endif |
| 3183 | 393 |
| 10311 | 394 // Determine ordering option. |
|
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
395 volatile char ord_job = 0; |
| 3183 | 396 static double safmin; |
| 3185 | 397 |
| 398 if (nargin == 2) | |
|
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
399 ord_job = 'N'; |
| 3185 | 400 else |
| 401 { | |
|
20700
68e3a747ca02
rename octave_value value extractors that accept error message args
John W. Eaton <jwe@octave.org>
parents:
20581
diff
changeset
|
402 std::string tmp = args(2).xstring_value ("qz: OPT must be a string"); |
|
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
403 |
|
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
404 if (! tmp.empty ()) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
405 ord_job = tmp[0]; |
| 3183 | 406 |
|
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
407 if (! (ord_job == 'N' || ord_job == 'n' |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
408 || ord_job == 'S' || ord_job == 's' |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
409 || ord_job == 'B' || ord_job == 'b' |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
410 || ord_job == '+' || ord_job == '-')) |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
411 error ("qz: invalid order option"); |
| 3185 | 412 |
| 413 // overflow constant required by dlag2 | |
| 4552 | 414 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1), |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
415 safmin |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
416 F77_CHAR_ARG_LEN (1)); |
| 3183 | 417 |
| 3185 | 418 #ifdef DEBUG_EIG |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
419 std::cout << "qz: initial value of safmin=" |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
420 << setiosflags (std::ios::scientific) |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
421 << safmin << std::endl; |
| 3185 | 422 #endif |
| 3183 | 423 |
| 10311 | 424 // Some machines (e.g., DEC alpha) get safmin = 0; |
| 425 // for these, use eps instead to avoid problems in dlag2. | |
| 3185 | 426 if (safmin == 0) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
427 { |
| 3185 | 428 #ifdef DEBUG_EIG |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
429 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; |
| 3185 | 430 #endif |
| 3183 | 431 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
432 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
433 safmin |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
434 F77_CHAR_ARG_LEN (1)); |
| 3185 | 435 |
| 436 #ifdef DEBUG_EIG | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
437 std::cout << "qz: safmin set to " |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
438 << setiosflags (std::ios::scientific) |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
439 << safmin << std::endl; |
| 3185 | 440 #endif |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
441 } |
| 3183 | 442 } |
| 443 | |
| 3185 | 444 #ifdef DEBUG |
| 3538 | 445 std::cout << "qz: check argument 1" << std::endl; |
| 3185 | 446 #endif |
| 3183 | 447 |
| 10311 | 448 // Argument 1: check if it's o.k. dimensioned. |
| 5275 | 449 octave_idx_type nn = args(0).rows (); |
| 3185 | 450 |
| 451 #ifdef DEBUG | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
452 std::cout << "argument 1 dimensions: (" |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
453 << nn << "," << args(0).columns () << ")" |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
454 << std::endl; |
| 3185 | 455 #endif |
| 456 | |
| 20892 | 457 octave_value_list retval; |
| 458 | |
| 3185 | 459 int arg_is_empty = empty_arg ("qz", nn, args(0).columns ()); |
| 460 | |
| 3183 | 461 if (arg_is_empty < 0) |
| 3185 | 462 { |
|
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
463 warn_empty_arg ("qz: parameter 1"); |
| 3185 | 464 return retval; |
| 465 } | |
| 3183 | 466 else if (arg_is_empty > 0) |
| 3185 | 467 { |
|
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
468 warn_empty_arg ("qz: parameter 1; continuing"); |
| 3185 | 469 return octave_value_list (2, Matrix ()); |
| 470 } | |
| 471 else if (args(0).columns () != nn) | |
|
21110
3d0d84305600
Use err_square_matrix_required more widely.
Rik <rik@octave.org>
parents:
21100
diff
changeset
|
472 err_square_matrix_required ("qz", "A"); |
| 3183 | 473 |
| 10311 | 474 // Argument 1: dimensions look good; get the value. |
| 3183 | 475 Matrix aa; |
| 476 ComplexMatrix caa; | |
| 3185 | 477 |
| 478 if (args(0).is_complex_type ()) | |
| 3183 | 479 caa = args(0).complex_matrix_value (); |
| 3185 | 480 else |
| 3183 | 481 aa = args(0).matrix_value (); |
| 3185 | 482 |
| 483 #ifdef DEBUG | |
| 3538 | 484 std::cout << "qz: check argument 2" << std::endl; |
| 3185 | 485 #endif |
| 3183 | 486 |
| 10311 | 487 // Extract argument 2 (bb, or cbb if complex). |
| 3185 | 488 if ((nn != args(1).columns ()) || (nn != args(1).rows ())) |
|
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
489 err_nonconformant (); |
| 3185 | 490 |
| 3183 | 491 Matrix bb; |
| 492 ComplexMatrix cbb; | |
| 3185 | 493 |
| 494 if (args(1).is_complex_type ()) | |
| 3183 | 495 cbb = args(1).complex_matrix_value (); |
| 496 else | |
| 497 bb = args(1).matrix_value (); | |
| 3185 | 498 |
| 3183 | 499 // Both matrices loaded, now let's check what kind of arithmetic: |
| 10311 | 500 // declared volatile to avoid compiler warnings about long jumps, |
| 501 // vforks. | |
| 3185 | 502 |
| 10311 | 503 volatile int complex_case |
| 3185 | 504 = (args(0).is_complex_type () || args(1).is_complex_type ()); |
| 10311 | 505 |
| 3185 | 506 if (nargin == 3 && complex_case) |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
507 error ("qz: cannot re-order complex qz decomposition"); |
| 3183 | 508 |
| 10311 | 509 // First, declare variables used in both the real and complex case. |
| 3183 | 510 Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); |
| 511 RowVector alphar(nn), alphai(nn), betar(nn); | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
512 ComplexRowVector xalpha(nn), xbeta(nn); |
| 3185 | 513 ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn); |
| 5275 | 514 octave_idx_type ilo, ihi, info; |
| 3185 | 515 char compq = (nargout >= 3 ? 'V' : 'N'); |
|
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
516 char compz = ((nargout >= 4 || nargin == 3)? 'V' : 'N'); |
| 3183 | 517 |
| 10311 | 518 // Initialize Q, Z to identity if we need either of them. |
| 3185 | 519 if (compq == 'V' || compz == 'V') |
| 5275 | 520 for (octave_idx_type ii = 0; ii < nn; ii++) |
| 521 for (octave_idx_type jj = 0; jj < nn; jj++) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
522 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
523 OCTAVE_QUIT; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
524 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
525 } |
| 3183 | 526 |
| 10311 | 527 // Always perform permutation balancing. |
| 4552 | 528 const char bal_job = 'P'; |
| 10311 | 529 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn); |
| 3183 | 530 |
| 3185 | 531 if (complex_case) |
| 532 { | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
533 #ifdef DEBUG |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
534 if (compq == 'V') |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
535 std::cout << "qz: performing balancing; CQ=" << std::endl |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
536 << CQ << std::endl; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
537 #endif |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
538 if (args(0).is_real_type ()) |
| 10311 | 539 caa = ComplexMatrix (aa); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
540 |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
541 if (args(1).is_real_type ()) |
| 10311 | 542 cbb = ComplexMatrix (bb); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
543 |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
544 if (compq == 'V') |
| 10311 | 545 CQ = ComplexMatrix (QQ); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
546 |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
547 if (compz == 'V') |
| 10311 | 548 CZ = ComplexMatrix (ZZ); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
549 |
| 10311 | 550 F77_XFCN (zggbal, ZGGBAL, |
| 551 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
| 552 nn, caa.fortran_vec (), nn, cbb.fortran_vec (), | |
| 553 nn, ilo, ihi, lscale.fortran_vec (), | |
| 554 rscale.fortran_vec (), work.fortran_vec (), info | |
| 555 F77_CHAR_ARG_LEN (1))); | |
| 3185 | 556 } |
| 3183 | 557 else |
| 3185 | 558 { |
| 559 #ifdef DEBUG | |
| 560 if (compq == 'V') | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
561 std::cout << "qz: performing balancing; QQ=" << std::endl |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
562 << QQ << std::endl; |
| 3185 | 563 #endif |
| 3183 | 564 |
| 3185 | 565 F77_XFCN (dggbal, DGGBAL, |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
566 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
567 nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
568 nn, ilo, ihi, lscale.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
569 rscale.fortran_vec (), work.fortran_vec (), info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
570 F77_CHAR_ARG_LEN (1))); |
| 3185 | 571 } |
| 3183 | 572 |
| 573 // Since we just want the balancing matrices, we can use dggbal | |
| 10311 | 574 // for both the real and complex cases; left first |
| 3185 | 575 |
| 10311 | 576 #if 0 |
| 577 if (compq == 'V') | |
| 3185 | 578 { |
| 579 F77_XFCN (dggbak, DGGBAK, | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
580 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
581 F77_CONST_CHAR_ARG2 ("L", 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
582 nn, ilo, ihi, lscale.data (), rscale.data (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
583 nn, QQ.fortran_vec (), nn, info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
584 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
585 F77_CHAR_ARG_LEN (1))); |
| 3183 | 586 |
| 3185 | 587 #ifdef DEBUG |
| 588 if (compq == 'V') | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
589 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; |
| 3185 | 590 #endif |
| 3183 | 591 } |
| 592 | |
| 593 // then right | |
| 3185 | 594 if (compz == 'V') |
| 595 { | |
| 4552 | 596 F77_XFCN (dggbak, DGGBAK, |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
597 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
598 F77_CONST_CHAR_ARG2 ("R", 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
599 nn, ilo, ihi, lscale.data (), rscale.data (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
600 nn, ZZ.fortran_vec (), nn, info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
601 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
602 F77_CHAR_ARG_LEN (1))); |
| 3183 | 603 |
| 3185 | 604 #ifdef DEBUG |
| 605 if (compz == 'V') | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
606 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; |
| 3185 | 607 #endif |
| 10311 | 608 } |
| 609 #endif | |
| 3183 | 610 |
| 611 static char qz_job; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
612 qz_job = (nargout < 2 ? 'E' : 'S'); |
| 3185 | 613 |
| 3183 | 614 if (complex_case) |
| 3185 | 615 { |
| 10311 | 616 // Complex case. |
| 617 | |
| 618 // The QR decomposition of cbb. | |
| 619 ComplexQR cbqr (cbb); | |
| 620 // The R matrix of QR decomposition for cbb. | |
| 621 cbb = cbqr.R (); | |
| 622 // (Q*)caa for following work. | |
| 623 caa = (cbqr.Q ().hermitian ()) * caa; | |
| 624 CQ = CQ * cbqr.Q (); | |
| 625 | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
626 F77_XFCN (zgghrd, ZGGHRD, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
627 (F77_CONST_CHAR_ARG2 (&compq, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
628 F77_CONST_CHAR_ARG2 (&compz, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
629 nn, ilo, ihi, caa.fortran_vec (), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
630 nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
631 CZ.fortran_vec (), nn, info |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
632 F77_CHAR_ARG_LEN (1) |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
633 F77_CHAR_ARG_LEN (1))); |
| 10311 | 634 |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
635 ComplexRowVector cwork (1 * nn); |
| 10311 | 636 |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
637 F77_XFCN (zhgeqz, ZHGEQZ, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
638 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
639 F77_CONST_CHAR_ARG2 (&compq, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
640 F77_CONST_CHAR_ARG2 (&compz, 1), |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
641 nn, ilo, ihi, |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
642 caa.fortran_vec (), nn, |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
643 cbb.fortran_vec (),nn, |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
644 xalpha.fortran_vec (), xbeta.fortran_vec (), |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
645 CQ.fortran_vec (), nn, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
646 CZ.fortran_vec (), nn, |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
647 cwork.fortran_vec (), nn, rwork.fortran_vec (), info |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
648 F77_CHAR_ARG_LEN (1) |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
649 F77_CHAR_ARG_LEN (1) |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
650 F77_CHAR_ARG_LEN (1))); |
| 3185 | 651 |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
652 if (compq == 'V') |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
653 { |
| 10311 | 654 // Left eigenvector. |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
655 F77_XFCN (zggbak, ZGGBAK, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
656 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
657 F77_CONST_CHAR_ARG2 ("L", 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
658 nn, ilo, ihi, lscale.data (), rscale.data (), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
659 nn, CQ.fortran_vec (), nn, info |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
660 F77_CHAR_ARG_LEN (1) |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
661 F77_CHAR_ARG_LEN (1))); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
662 } |
| 3185 | 663 |
| 10311 | 664 // Right eigenvector. |
| 3185 | 665 if (compz == 'V') |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
666 { |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
667 F77_XFCN (zggbak, ZGGBAK, |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
668 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
669 F77_CONST_CHAR_ARG2 ("R", 1), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
670 nn, ilo, ihi, lscale.data (), rscale.data (), |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
671 nn, CZ.fortran_vec (), nn, info |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
672 F77_CHAR_ARG_LEN (1) |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
673 F77_CHAR_ARG_LEN (1))); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
674 } |
| 3185 | 675 |
| 676 } | |
| 10311 | 677 else |
| 3185 | 678 { |
| 679 #ifdef DEBUG | |
| 3538 | 680 std::cout << "qz: peforming qr decomposition of bb" << std::endl; |
| 3185 | 681 #endif |
| 3183 | 682 |
| 10311 | 683 // Compute the QR factorization of bb. |
| 3185 | 684 QR bqr (bb); |
| 685 | |
| 686 #ifdef DEBUG | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
687 std::cout << "qz: qr (bb) done; now peforming qz decomposition" |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
688 << std::endl; |
| 3185 | 689 #endif |
| 3183 | 690 |
| 3185 | 691 bb = bqr.R (); |
| 692 | |
| 693 #ifdef DEBUG | |
| 3538 | 694 std::cout << "qz: extracted bb" << std::endl; |
| 3185 | 695 #endif |
| 3183 | 696 |
| 10311 | 697 aa = (bqr.Q ()).transpose () * aa; |
| 3185 | 698 |
| 699 #ifdef DEBUG | |
| 3538 | 700 std::cout << "qz: updated aa " << std::endl; |
| 701 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; | |
| 3183 | 702 |
| 3185 | 703 if (compq == 'V') |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
704 std::cout << "QQ =" << QQ << std::endl; |
| 3185 | 705 #endif |
| 3183 | 706 |
| 3185 | 707 if (compq == 'V') |
| 10311 | 708 QQ = QQ * bqr.Q (); |
| 3183 | 709 |
| 3185 | 710 #ifdef DEBUG |
| 3538 | 711 std::cout << "qz: precursors done..." << std::endl; |
| 3185 | 712 #endif |
| 3183 | 713 |
| 3185 | 714 #ifdef DEBUG |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
715 std::cout << "qz: compq = " << compq << ", compz = " << compz |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
716 << std::endl; |
| 3185 | 717 #endif |
| 3183 | 718 |
| 10311 | 719 // Reduce to generalized hessenberg form. |
| 3185 | 720 F77_XFCN (dgghrd, DGGHRD, |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
721 (F77_CONST_CHAR_ARG2 (&compq, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
722 F77_CONST_CHAR_ARG2 (&compz, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
723 nn, ilo, ihi, aa.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
724 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
725 ZZ.fortran_vec (), nn, info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
726 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
727 F77_CHAR_ARG_LEN (1))); |
| 3183 | 728 |
| 10311 | 729 // Check if just computing generalized eigenvalues or if we're |
| 730 // actually computing the decomposition. | |
| 3183 | 731 |
| 10311 | 732 // Reduce to generalized Schur form. |
| 3185 | 733 F77_XFCN (dhgeqz, DHGEQZ, |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
734 (F77_CONST_CHAR_ARG2 (&qz_job, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
735 F77_CONST_CHAR_ARG2 (&compq, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
736 F77_CONST_CHAR_ARG2 (&compz, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
737 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
738 nn, alphar.fortran_vec (), alphai.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
739 betar.fortran_vec (), QQ.fortran_vec (), nn, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
740 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
741 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
742 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
743 F77_CHAR_ARG_LEN (1))); |
| 10311 | 744 |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
745 if (compq == 'V') |
| 10311 | 746 { |
| 747 F77_XFCN (dggbak, DGGBAK, | |
| 748 (F77_CONST_CHAR_ARG2 (&bal_job, 1), | |
| 749 F77_CONST_CHAR_ARG2 ("L", 1), | |
| 750 nn, ilo, ihi, lscale.data (), rscale.data (), | |
| 751 nn, QQ.fortran_vec (), nn, info | |
| 752 F77_CHAR_ARG_LEN (1) | |
| 753 F77_CHAR_ARG_LEN (1))); | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
754 |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
755 #ifdef DEBUG |
| 10311 | 756 if (compq == 'V') |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
757 std::cout << "qz: balancing done; QQ=" << std::endl |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
758 << QQ << std::endl; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
759 #endif |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
760 } |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
761 |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
762 // then right |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
763 if (compz == 'V') |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
764 { |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
765 F77_XFCN (dggbak, DGGBAK, |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
766 (F77_CONST_CHAR_ARG2 (&bal_job, 1), |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
767 F77_CONST_CHAR_ARG2 ("R", 1), |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
768 nn, ilo, ihi, lscale.data (), rscale.data (), |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
769 nn, ZZ.fortran_vec (), nn, info |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
770 F77_CHAR_ARG_LEN (1) |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
771 F77_CHAR_ARG_LEN (1))); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
772 |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
773 #ifdef DEBUG |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
774 if (compz == 'V') |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
775 std::cout << "qz: balancing done; ZZ=" << std::endl |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
776 << ZZ << std::endl; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
777 #endif |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
778 } |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
779 |
| 3185 | 780 } |
| 3183 | 781 |
| 10311 | 782 // Order the QZ decomposition? |
|
8927
d75f4ee0538d
qz.cc (Fqz): avoid maybe clobbred by vfork warning from GCC
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
783 if (! (ord_job == 'N' || ord_job == 'n')) |
| 3183 | 784 { |
| 3185 | 785 if (complex_case) |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
786 // Probably not needed, but better be safe. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
787 error ("qz: cannot re-order complex qz decomposition"); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
788 |
| 3185 | 789 #ifdef DEBUG_SORT |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
790 std::cout << "qz: ordering eigenvalues: ord_job = " |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
791 << ord_job << std::endl; |
| 3185 | 792 #endif |
| 3183 | 793 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
794 // Declared static to avoid vfork/long jump compiler complaints. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
795 static sort_function sort_test; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
796 sort_test = 0; |
| 3183 | 797 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
798 switch (ord_job) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
799 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
800 case 'S': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
801 case 's': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
802 sort_test = &fin; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
803 break; |
| 3183 | 804 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
805 case 'B': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
806 case 'b': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
807 sort_test = &fout; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
808 break; |
| 3183 | 809 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
810 case '+': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
811 sort_test = &fcrhp; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
812 break; |
| 3183 | 813 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
814 case '-': |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
815 sort_test = &folhp; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
816 break; |
| 3185 | 817 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
818 default: |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
819 // Invalid order option (should never happen, since we |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
820 // checked the options at the top). |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
821 panic_impossible (); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
822 break; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
823 } |
| 3183 | 824 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
825 octave_idx_type ndim, fail; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
826 double inf_norm; |
| 3185 | 827 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
828 F77_XFCN (xdlange, XDLANGE, |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
829 (F77_CONST_CHAR_ARG2 ("I", 1), |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
830 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
831 F77_CHAR_ARG_LEN (1))); |
| 3185 | 832 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
833 double eps = std::numeric_limits<double>::epsilon () * inf_norm * nn; |
| 3185 | 834 |
| 835 #ifdef DEBUG_SORT | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
836 std::cout << "qz: calling dsubsp: aa=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
837 octave_print_internal (std::cout, aa, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
838 std::cout << std::endl << "bb=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
839 octave_print_internal (std::cout, bb, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
840 if (compz == 'V') |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
841 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
842 std::cout << std::endl << "ZZ=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
843 octave_print_internal (std::cout, ZZ, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
844 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
845 std::cout << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
846 std::cout << "alphar = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
847 octave_print_internal (std::cout, (Matrix) alphar, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
848 std::cout << std::endl << "alphai = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
849 octave_print_internal (std::cout, (Matrix) alphai, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
850 std::cout << std::endl << "beta = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
851 octave_print_internal (std::cout, (Matrix) betar, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
852 std::cout << std::endl; |
| 3185 | 853 #endif |
| 854 | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
855 Array<octave_idx_type> ind (dim_vector (nn, 1)); |
| 3550 | 856 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
857 F77_XFCN (dsubsp, DSUBSP, |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
858 (nn, nn, aa.fortran_vec (), bb.fortran_vec (), |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
859 ZZ.fortran_vec (), sort_test, eps, ndim, fail, |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
860 ind.fortran_vec ())); |
| 3185 | 861 |
| 862 #ifdef DEBUG | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
863 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
864 octave_print_internal (std::cout, aa, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
865 std::cout << std::endl << "bb=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
866 octave_print_internal (std::cout, bb, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
867 if (compz == 'V') |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
868 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
869 std::cout << std::endl << "ZZ=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
870 octave_print_internal (std::cout, ZZ, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
871 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
872 std::cout << std::endl; |
| 3185 | 873 #endif |
| 874 | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
875 // Manually update alphar, alphai, betar. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
876 static int jj; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
877 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
878 jj = 0; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
879 while (jj < nn) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
880 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
881 #ifdef DEBUG_EIG |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
882 std::cout << "computing gen eig #" << jj << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
883 #endif |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
884 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
885 // Number of zeros in this block. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
886 static int zcnt; |
| 3185 | 887 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
888 if (jj == (nn-1)) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
889 zcnt = 1; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
890 else if (aa(jj+1,jj) == 0) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
891 zcnt = 1; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
892 else zcnt = 2; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
893 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
894 if (zcnt == 1) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
895 { |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
896 // Real zero. |
| 3185 | 897 #ifdef DEBUG_EIG |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
898 std::cout << " single gen eig:" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
899 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
900 << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
901 std::cout << " betar(" << jj << ") = " << bb(jj,jj) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
902 << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
903 std::cout << " alphai(" << jj << ") = 0" << std::endl; |
| 3185 | 904 #endif |
| 905 | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
906 alphar(jj) = aa(jj,jj); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
907 alphai(jj) = 0; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
908 betar(jj) = bb(jj,jj); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
909 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
910 else |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
911 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
912 // Complex conjugate pair. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
913 #ifdef DEBUG_EIG |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
914 std::cout << "qz: calling dlag2:" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
915 std::cout << "safmin=" |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
916 << setiosflags (std::ios::scientific) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
917 << safmin << std::endl; |
| 3185 | 918 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
919 for (int idr = jj; idr <= jj+1; idr++) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
920 { |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
921 for (int idc = jj; idc <= jj+1; idc++) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
922 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
923 std::cout << "aa(" << idr << "," << idc << ")=" |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
924 << aa(idr,idc) << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
925 std::cout << "bb(" << idr << "," << idc << ")=" |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
926 << bb(idr,idc) << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
927 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
928 } |
| 3185 | 929 #endif |
| 930 | |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
931 // FIXME: probably should be using |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
932 // fortran_vec instead of &aa(jj,jj) here. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
933 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
934 double scale1, scale2, wr1, wr2, wi; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
935 const double *aa_ptr = aa.data () + jj * nn + jj; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
936 const double *bb_ptr = bb.data () + jj * nn + jj; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
937 F77_XFCN (dlag2, DLAG2, |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
938 (aa_ptr, nn, bb_ptr, nn, safmin, |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
939 scale1, scale2, wr1, wr2, wi)); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
940 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
941 #ifdef DEBUG_EIG |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
942 std::cout << "dlag2 returns: scale1=" << scale1 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
943 << "\tscale2=" << scale2 << std::endl |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
944 << "\twr1=" << wr1 << "\twr2=" << wr2 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
945 << "\twi=" << wi << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
946 #endif |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
947 |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
948 // Just to be safe, check if it's a real pair. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
949 if (wi == 0) |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
950 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
951 alphar(jj) = wr1; |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
952 alphai(jj) = 0; |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
953 betar(jj) = scale1; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
954 alphar(jj+1) = wr2; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
955 alphai(jj+1) = 0; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
956 betar(jj+1) = scale2; |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
957 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
958 else |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
959 { |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
960 alphar(jj) = alphar(jj+1) = wr1; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
961 alphai(jj) = -(alphai(jj+1) = wi); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
962 betar(jj) = betar(jj+1) = scale1; |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
963 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
964 } |
| 3185 | 965 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
966 // Advance past this block. |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
967 jj += zcnt; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
968 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
969 |
| 3185 | 970 #ifdef DEBUG_SORT |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
971 std::cout << "qz: back from dsubsp: aa=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
972 octave_print_internal (std::cout, aa, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
973 std::cout << std::endl << "bb=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
974 octave_print_internal (std::cout, bb, 0); |
| 3185 | 975 |
|
20831
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
976 if (compz == 'V') |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
977 { |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
978 std::cout << std::endl << "ZZ=" << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
979 octave_print_internal (std::cout, ZZ, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
980 } |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
981 std::cout << std::endl << "qz: ndim=" << ndim << std::endl |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
982 << "fail=" << fail << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
983 std::cout << "alphar = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
984 octave_print_internal (std::cout, (Matrix) alphar, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
985 std::cout << std::endl << "alphai = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
986 octave_print_internal (std::cout, (Matrix) alphai, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
987 std::cout << std::endl << "beta = " << std::endl; |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
988 octave_print_internal (std::cout, (Matrix) betar, 0); |
|
35241c4b696c
eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents:
20802
diff
changeset
|
989 std::cout << std::endl; |
| 3185 | 990 #endif |
| 3183 | 991 } |
| 3185 | 992 |
| 10311 | 993 // Compute generalized eigenvalues? |
| 3183 | 994 ComplexColumnVector gev; |
| 3185 | 995 |
| 996 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) | |
| 3183 | 997 { |
| 3185 | 998 if (complex_case) |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
999 { |
| 10311 | 1000 int cnt = 0; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1001 |
| 10311 | 1002 for (int ii = 0; ii < nn; ii++) |
| 1003 cnt++; | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1004 |
| 10311 | 1005 ComplexColumnVector tmp (cnt); |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1006 |
| 10311 | 1007 cnt = 0; |
| 1008 for (int ii = 0; ii < nn; ii++) | |
| 1009 tmp(cnt++) = xalpha(ii) / xbeta(ii); | |
| 1010 | |
| 1011 gev = tmp; | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1012 } |
| 3185 | 1013 else |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1014 { |
| 3185 | 1015 #ifdef DEBUG |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1016 std::cout << "qz: computing generalized eigenvalues" << std::endl; |
| 3185 | 1017 #endif |
| 3183 | 1018 |
| 10311 | 1019 // Return finite generalized eigenvalues. |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1020 int cnt = 0; |
| 3185 | 1021 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1022 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1023 if (betar(ii) != 0) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1024 cnt++; |
| 3185 | 1025 |
| 10311 | 1026 ComplexColumnVector tmp (cnt); |
| 3185 | 1027 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1028 cnt = 0; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1029 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1030 if (betar(ii) != 0) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1031 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii); |
| 10311 | 1032 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1033 gev = tmp; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1034 } |
| 3183 | 1035 } |
| 1036 | |
| 10311 | 1037 // Right, left eigenvector matrices. |
| 3185 | 1038 if (nargout >= 5) |
| 3183 | 1039 { |
| 10311 | 1040 // Which side to compute? |
| 1041 char side = (nargout == 5 ? 'R' : 'B'); | |
| 1042 // Compute all of them and backtransform | |
| 1043 char howmny = 'B'; | |
| 1044 // Dummy pointer; select is not used. | |
| 1045 octave_idx_type *select = 0; | |
| 3185 | 1046 |
| 1047 if (complex_case) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1048 { |
| 10311 | 1049 CVL = CQ; |
| 1050 CVR = CZ; | |
| 1051 ComplexRowVector cwork2 (2 * nn); | |
| 1052 RowVector rwork2 (8 * nn); | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1053 octave_idx_type m; |
|
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1054 |
| 10311 | 1055 F77_XFCN (ztgevc, ZTGEVC, |
| 1056 (F77_CONST_CHAR_ARG2 (&side, 1), | |
| 1057 F77_CONST_CHAR_ARG2 (&howmny, 1), | |
| 1058 select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (), | |
| 1059 nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn, | |
| 1060 m, cwork2.fortran_vec (), rwork2.fortran_vec (), info | |
| 1061 F77_CHAR_ARG_LEN (1) | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
1062 F77_CHAR_ARG_LEN (1))); |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1063 } |
| 3185 | 1064 else |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1065 { |
| 3185 | 1066 #ifdef DEBUG |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1067 std::cout << "qz: computing generalized eigenvectors" << std::endl; |
| 3185 | 1068 #endif |
| 1069 | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1070 VL = QQ; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1071 VR = ZZ; |
| 10311 | 1072 octave_idx_type m; |
| 1073 | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1074 F77_XFCN (dtgevc, DTGEVC, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1075 (F77_CONST_CHAR_ARG2 (&side, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1076 F77_CONST_CHAR_ARG2 (&howmny, 1), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1077 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1078 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1079 m, work.fortran_vec (), info |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1080 F77_CHAR_ARG_LEN (1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1081 F77_CHAR_ARG_LEN (1))); |
| 3185 | 1082 |
| 10311 | 1083 // Now construct the complex form of VV, WW. |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1084 int jj = 0; |
| 3185 | 1085 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1086 while (jj < nn) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1087 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1088 OCTAVE_QUIT; |
| 4153 | 1089 |
| 10311 | 1090 // See if real or complex eigenvalue. |
| 1091 | |
| 1092 // Column increment; assume complex eigenvalue. | |
| 1093 int cinc = 2; | |
| 3185 | 1094 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1095 if (jj == (nn-1)) |
| 10311 | 1096 // Single column. |
| 1097 cinc = 1; | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1098 else if (aa(jj+1,jj) == 0) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1099 cinc = 1; |
| 3185 | 1100 |
| 10311 | 1101 // Now copy the eigenvector (s) to CVR, CVL. |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1102 if (cinc == 1) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1103 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1104 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1105 CVR(ii,jj) = VR(ii,jj); |
| 3185 | 1106 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1107 if (side == 'B') |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1108 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1109 CVL(ii,jj) = VL(ii,jj); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1110 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1111 else |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1112 { |
| 10311 | 1113 // Double column; complex vector. |
| 3185 | 1114 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1115 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1116 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1117 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1118 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1119 } |
| 3183 | 1120 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1121 if (side == 'B') |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1122 for (int ii = 0; ii < nn; ii++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1123 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1124 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1125 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1126 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1127 } |
| 3185 | 1128 |
| 10311 | 1129 // Advance to next eigenvectors (if any). |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1130 jj += cinc; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1131 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1132 } |
| 3183 | 1133 } |
| 3185 | 1134 |
| 1135 switch (nargout) | |
| 1136 { | |
| 1137 case 7: | |
| 1138 retval(6) = gev; | |
| 1139 | |
| 10311 | 1140 case 6: |
| 1141 // Return eigenvectors. | |
| 3185 | 1142 retval(5) = CVL; |
| 1143 | |
| 10311 | 1144 case 5: |
| 1145 // Return eigenvectors. | |
| 3185 | 1146 retval(4) = CVR; |
| 1147 | |
| 1148 case 4: | |
| 1149 if (nargin == 3) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1150 { |
| 3185 | 1151 #ifdef DEBUG |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1152 std::cout << "qz: sort: retval(3) = gev = " << std::endl; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1153 octave_print_internal (std::cout, gev); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1154 std::cout << std::endl; |
| 3185 | 1155 #endif |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1156 retval(3) = gev; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
1157 } |
| 3185 | 1158 else |
| 10311 | 1159 { |
| 1160 if (complex_case) | |
| 1161 retval(3) = CZ; | |
| 1162 else | |
| 1163 retval(3) = ZZ; | |
| 1164 } | |
| 3185 | 1165 |
| 1166 case 3: | |
| 1167 if (nargin == 3) | |
|
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1168 { |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1169 if (complex_case) |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1170 retval(2) = CZ; |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1171 else |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1172 retval(2) = ZZ; |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1173 } |
| 3185 | 1174 else |
| 10311 | 1175 { |
| 1176 if (complex_case) | |
| 1177 retval(2) = CQ.hermitian (); | |
| 1178 else | |
| 1179 retval(2) = QQ.transpose (); | |
| 1180 } | |
| 1181 | |
| 3185 | 1182 case 2: |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1183 { |
| 10311 | 1184 if (complex_case) |
| 1185 { | |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1186 #ifdef DEBUG |
|
14844
5bc9b9cb4362
maint: Use Octave coding conventions for cuddled parenthesis in retval assignments.
Rik <octave@nomad.inbox5.com>
parents:
14501
diff
changeset
|
1187 std::cout << "qz: retval(1) = cbb = " << std::endl; |
| 10311 | 1188 octave_print_internal (std::cout, cbb, 0); |
| 1189 std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl; | |
| 1190 octave_print_internal (std::cout, caa, 0); | |
| 1191 std::cout << std::endl; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
1192 #endif |
| 10311 | 1193 retval(1) = cbb; |
| 1194 retval(0) = caa; | |
| 1195 } | |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1196 else |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1197 { |
| 3185 | 1198 #ifdef DEBUG |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1199 std::cout << "qz: retval(1) = bb = " << std::endl; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1200 octave_print_internal (std::cout, bb, 0); |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1201 std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1202 octave_print_internal (std::cout, aa, 0); |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1203 std::cout << std::endl; |
| 3185 | 1204 #endif |
|
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1205 retval(1) = bb; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1206 retval(0) = aa; |
|
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
1207 } |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
1208 } |
| 10311 | 1209 break; |
|
10307
4e4270ab70d6
support complex case in qz
Jyh-miin Lin <jyhmiinlin@gmail.com>
parents:
10155
diff
changeset
|
1210 |
| 3185 | 1211 case 1: |
| 1212 case 0: | |
| 1213 #ifdef DEBUG | |
| 3538 | 1214 std::cout << "qz: retval(0) = gev = " << gev << std::endl; |
| 3185 | 1215 #endif |
| 1216 retval(0) = gev; | |
| 1217 break; | |
| 1218 | |
| 1219 default: | |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
1220 error ("qz: too many return arguments"); |
| 3185 | 1221 break; |
| 3183 | 1222 } |
| 1223 | |
| 3185 | 1224 #ifdef DEBUG |
| 3538 | 1225 std::cout << "qz: exiting (at long last)" << std::endl; |
| 3185 | 1226 #endif |
| 3183 | 1227 |
| 1228 return retval; | |
| 1229 } | |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1230 |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1231 /* |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1232 %!shared a, b, c |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1233 %! a = [1 2; 0 3]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1234 %! b = [1 0; 0 0]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1235 %! c = [0 1; 0 0]; |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1236 %!assert (qz (a,b), 1) |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1237 %!assert (isempty (qz (a,c))) |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1238 |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1239 ## Exaple 7.7.3 in Golub & Van Loan |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1240 %!test |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1241 %! a = [ 10 1 2; |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1242 %! 1 2 -1; |
|
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1243 %! 1 1 2]; |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1244 %! b = reshape (1:9,3,3); |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1245 %! [aa, bb, q, z, v, w, lambda] = qz (a, b); |
|
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1246 %! sz = length (lambda); |
|
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14360
diff
changeset
|
1247 %! observed = (b * v * diag ([lambda;0])) (:, 1:sz); |
|
18678
6113e0c6920b
maint: Clean up extra spaces before/after parentheses.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
1248 %! assert ((a*v)(:, 1:sz), observed, norm (observed) * 1e-14); |
|
13088
4ffea87ad71b
codesprint: Fix tolerance for qz.cc tests.
Ben Abbott <bpabbott@mac.com>
parents:
13074
diff
changeset
|
1249 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :); |
|
18678
6113e0c6920b
maint: Clean up extra spaces before/after parentheses.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
1250 %! assert ((w'*a)(1:sz, :) , observed, norm (observed) * 1e-13); |
|
13088
4ffea87ad71b
codesprint: Fix tolerance for qz.cc tests.
Ben Abbott <bpabbott@mac.com>
parents:
13074
diff
changeset
|
1251 %! assert (q * a * z, aa, norm (aa) * 1e-14); |
|
4ffea87ad71b
codesprint: Fix tolerance for qz.cc tests.
Ben Abbott <bpabbott@mac.com>
parents:
13074
diff
changeset
|
1252 %! assert (q * b * z, bb, norm (bb) * 1e-14); |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1253 |
|
14876
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1254 %!test |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1255 %! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0]; |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1256 %! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1]; |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1257 %! [AA, BB, Q, Z1] = qz (A, B); |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1258 %! [AA, BB, Z2] = qz (A, B, '-'); |
|
54a386f2ac4e
correcly compute Z for 3-output call to qz (bug #36728)
John W. Eaton <jwe@octave.org>
parents:
14844
diff
changeset
|
1259 %! assert (Z1, Z2); |
|
13074
0b789a03bde1
codesprint: Add 3 tests for qz.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
12642
diff
changeset
|
1260 */ |