GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qz.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26// Generalized eigenvalue balancing via LAPACK
27
28// Originally written by A. S. Hodel <scotte@eng.auburn.edu>, but is
29// substantially different with the change to use LAPACK.
30
31#undef DEBUG
32#undef DEBUG_SORT
33#undef DEBUG_EIG
34
35#if defined (HAVE_CONFIG_H)
36# include "config.h"
37#endif
38
39#include <cctype>
40#include <cmath>
41
42#if defined (DEBUG_EIG)
43# include <iomanip>
44#endif
45
46#include "f77-fcn.h"
47#include "lo-lapack-proto.h"
48#include "qr.h"
49#include "quit.h"
50
51#include "defun.h"
52#include "error.h"
53#include "errwarn.h"
54#include "ovl.h"
55#if defined (DEBUG) || defined (DEBUG_SORT)
56# include "pager.h"
57# include "pr-output.h"
58#endif
59
60OCTAVE_NAMESPACE_BEGIN
61
62// FIXME: Matlab does not produce lambda as the first output argument.
63// Compatibility problem?
64
65DEFUN (qz, args, nargout,
66 doc: /* -*- texinfo -*-
67@deftypefn {} {@var{lambda} =} qz (@var{A}, @var{B})
68@deftypefnx {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] =} qz (@var{A}, @var{B})
69@deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}] =} qz (@var{A}, @var{B}, @var{opt})
70@deftypefnx {} {[@var{AA}, @var{BB}, @var{Z}, @var{lambda}] =} qz (@var{A}, @var{B}, @var{opt})
71Compute the QZ@tie{}decomposition of a generalized eigenvalue problem.
72
73The generalized eigenvalue problem is defined as
74
75@tex
76$$A x = \lambda B x$$
77@end tex
78@ifnottex
79
80@math{A x = @var{lambda} B x}
81
82@end ifnottex
83
84There are three calling forms of the function:
85
86@enumerate
87@item @code{@var{lambda} = qz (@var{A}, @var{B})}
88
89Compute the generalized eigenvalues
90@tex
91$\lambda.$
92@end tex
93@ifnottex
94@var{lambda}.
95@end ifnottex
96
97@item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] = qz (@var{A}, @var{B})}
98
99Compute QZ@tie{}decomposition, generalized eigenvectors, and generalized
100eigenvalues.
101@tex
102$$ AA = Q^T AZ, BB = Q^T BZ $$
103$$ AV = BV{ \rm diag }(\lambda) $$
104$$ W^T A = { \rm diag }(\lambda)W^T B $$
105@end tex
106@ifnottex
107
108@example
109@group
110
111@var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
112@var{A} * @var{V} = @var{B} * @var{V} * diag (@var{lambda})
113@var{W}' * @var{A} = diag (@var{lambda}) * @var{W}' * @var{B}
114
115@end group
116@end example
117
118@end ifnottex
119with @var{Q} and @var{Z} orthogonal (unitary for complex case).
120
121@item @code{[@var{AA}, @var{BB}, @var{Z} @{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}
122
123As in form 2 above, but allows ordering of generalized eigenpairs for, e.g.,
124solution of discrete time algebraic @nospell{Riccati} equations. Form 3 is not
125available for complex matrices, and does not compute the generalized
126eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.
127
128@table @var
129@item opt
130for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of
131the revised pencil contains all eigenvalues that satisfy:
132
133@table @asis
134@item @qcode{"N"}
135unordered (default)
136
137@item @qcode{"S"}
138small: leading block has all
139@tex
140$|\lambda| < 1$
141@end tex
142@ifnottex
143|@var{lambda}| < 1
144@end ifnottex
145
146@item @qcode{"B"}
147big: leading block has all
148@tex
149$|\lambda| \geq 1$
150@end tex
151@ifnottex
152|@var{lambda}| @geq{} 1
153@end ifnottex
154
155@item @qcode{"-"}
156negative real part: leading block has all eigenvalues in the open left
157half-plane
158
159@item @qcode{"+"}
160non-negative real part: leading block has all eigenvalues in the closed right
161half-plane
162@end table
163@end table
164@end enumerate
165
166Note: @code{qz} performs permutation balancing, but not scaling
167(@pxref{XREFbalance,,@code{balance}}), which may be lead to less accurate
168results than @code{eig}. The order of output arguments was selected for
169compatibility with @sc{matlab}.
170@seealso{eig, gsvd, balance, chol, hess, lu, qr, qzhess, schur}
171@end deftypefn */)
172{
173 int nargin = args.length ();
174
175#if defined (DEBUG)
176 octave_stdout << "qz: nargin = " << nargin
177 << ", nargout = " << nargout << std::endl;
178#endif
179
180 if (nargin < 2 || nargin > 3 || nargout > 7)
181 print_usage ();
182
183 if (nargin == 3 && (nargout < 3 || nargout > 4))
184 error ("qz: invalid number of output arguments for form 3 call");
185
186#if defined (DEBUG)
187 octave_stdout << "qz: determine ordering option" << std::endl;
188#endif
189
190 // Determine ordering option.
191
192 char ord_job = 'N';
193 double safmin = 0.0;
194
195 if (nargin == 3)
196 {
197 std::string opt = args(2).xstring_value ("qz: OPT must be a string");
198
199 if (opt.empty ())
200 error ("qz: OPT must be a non-empty string");
201
202 ord_job = std::toupper (opt[0]);
203
204 std::string valid_opts = "NSB+-";
205
206 if (valid_opts.find_first_of (ord_job) == std::string::npos)
207 error ("qz: invalid order option '%c'", ord_job);
208
209 // overflow constant required by dlag2
210 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
211 safmin
212 F77_CHAR_ARG_LEN (1));
213
214#if defined (DEBUG_EIG)
215 octave_stdout << "qz: initial value of safmin="
216 << setiosflags (std::ios::scientific)
217 << safmin << std::endl;
218#endif
219
220 // Some machines (e.g., DEC alpha) get safmin = 0;
221 // for these, use eps instead to avoid problems in dlag2.
222 if (safmin == 0)
223 {
224#if defined (DEBUG_EIG)
225 octave_stdout << "qz: DANGER WILL ROBINSON: safmin is 0!"
226 << std::endl;
227#endif
228
229 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
230 safmin
231 F77_CHAR_ARG_LEN (1));
232
233#if defined (DEBUG_EIG)
234 octave_stdout << "qz: safmin set to "
235 << setiosflags (std::ios::scientific)
236 << safmin << std::endl;
237#endif
238 }
239 }
240
241#if defined (DEBUG)
242 octave_stdout << "qz: check matrix A" << std::endl;
243#endif
244
245 // Matrix A: check dimensions.
246 F77_INT nn = to_f77_int (args(0).rows ());
247 F77_INT nc = to_f77_int (args(0).columns ());
248
249#if defined (DEBUG)
250 octave_stdout << "Matrix A dimensions: (" << nn << ',' << nc << ')'
251 << std::endl;
252#endif
253
254 if (args(0).isempty ())
255 {
256 warn_empty_arg ("qz: A");
257 return octave_value_list (2, Matrix ());
258 }
259 else if (nc != nn)
260 err_square_matrix_required ("qz", "A");
261
262 // Matrix A: get value.
263 Matrix aa;
264 ComplexMatrix caa;
265
266 if (args(0).iscomplex ())
267 caa = args(0).complex_matrix_value ();
268 else
269 aa = args(0).matrix_value ();
270
271#if defined (DEBUG)
272 octave_stdout << "qz: check matrix B" << std::endl;
273#endif
274
275 // Extract argument 2 (bb, or cbb if complex).
276 F77_INT b_nr = to_f77_int (args(1).rows ());
277 F77_INT b_nc = to_f77_int (args(1).columns ());
278
279 if (nn != b_nc || nn != b_nr)
280 ::err_nonconformant ();
281
282 Matrix bb;
283 ComplexMatrix cbb;
284
285 if (args(1).iscomplex ())
286 cbb = args(1).complex_matrix_value ();
287 else
288 bb = args(1).matrix_value ();
289
290 // Both matrices loaded, now check whether to calculate complex or real val.
291
292 bool complex_case = (args(0).iscomplex () || args(1).iscomplex ());
293
294 if (nargin == 3 && complex_case)
295 error ("qz: cannot re-order complex qz decomposition");
296
297 // First, declare variables used in both the real and complex cases.
298 // FIXME: There are a lot of excess variables declared.
299 // Probably a better way to handle this.
300 Matrix QQ (nn, nn), ZZ (nn, nn), VR (nn, nn), VL (nn, nn);
301 RowVector alphar (nn), alphai (nn), betar (nn);
302 ComplexRowVector xalpha (nn), xbeta (nn);
303 ComplexMatrix CQ (nn, nn), CZ (nn, nn), CVR (nn, nn), CVL (nn, nn);
304 F77_INT ilo, ihi, info;
305 char comp_q = (nargout >= 3 ? 'V' : 'N');
306 char comp_z = ((nargout >= 4 || nargin == 3)? 'V' : 'N');
307
308 // Initialize Q, Z to identity matrix if either is needed
309 if (comp_q == 'V' || comp_z == 'V')
310 {
311 double *QQptr = QQ.fortran_vec ();
312 double *ZZptr = ZZ.fortran_vec ();
313 std::fill_n (QQptr, QQ.numel (), 0.0);
314 std::fill_n (ZZptr, ZZ.numel (), 0.0);
315 for (F77_INT i = 0; i < nn; i++)
316 {
317 QQ(i, i) = 1.0;
318 ZZ(i, i) = 1.0;
319 }
320 }
321
322 // Always perform permutation balancing.
323 const char bal_job = 'P';
324 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
325
326 if (complex_case)
327 {
328#if defined (DEBUG)
329 if (comp_q == 'V')
330 octave_stdout << "qz: performing balancing; CQ =\n" << CQ << std::endl;
331#endif
332 if (args(0).isreal ())
333 caa = ComplexMatrix (aa);
334
335 if (args(1).isreal ())
336 cbb = ComplexMatrix (bb);
337
338 if (comp_q == 'V')
339 CQ = ComplexMatrix (QQ);
340
341 if (comp_z == 'V')
342 CZ = ComplexMatrix (ZZ);
343
344 F77_XFCN (zggbal, ZGGBAL,
345 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
346 nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
347 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
348 nn, ilo, ihi, lscale.fortran_vec (),
349 rscale.fortran_vec (), work.fortran_vec (), info
350 F77_CHAR_ARG_LEN (1)));
351 }
352 else
353 {
354#if defined (DEBUG)
355 if (comp_q == 'V')
356 octave_stdout << "qz: performing balancing; QQ =\n" << QQ << std::endl;
357#endif
358
359 F77_XFCN (dggbal, DGGBAL,
360 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
361 nn, aa.fortran_vec (), nn, bb.fortran_vec (),
362 nn, ilo, ihi, lscale.fortran_vec (),
363 rscale.fortran_vec (), work.fortran_vec (), info
364 F77_CHAR_ARG_LEN (1)));
365 }
366
367 // Only permutation balance above is done. Skip scaling balance.
368
369#if 0
370 // Since we just want the balancing matrices, we can use dggbal
371 // for both the real and complex cases; left first
372
373 if (comp_q == 'V')
374 {
375 F77_XFCN (dggbak, DGGBAK,
376 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
377 F77_CONST_CHAR_ARG2 ("L", 1),
378 nn, ilo, ihi, lscale.data (), rscale.data (),
379 nn, QQ.fortran_vec (), nn, info
380 F77_CHAR_ARG_LEN (1)
381 F77_CHAR_ARG_LEN (1)));
382
383#if defined (DEBUG)
384 if (comp_q == 'V')
385 octave_stdout << "qz: balancing done; QQ =\n" << QQ << std::endl;
386#endif
387 }
388
389 // then right
390 if (comp_z == 'V')
391 {
392 F77_XFCN (dggbak, DGGBAK,
393 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
394 F77_CONST_CHAR_ARG2 ("R", 1),
395 nn, ilo, ihi, lscale.data (), rscale.data (),
396 nn, ZZ.fortran_vec (), nn, info
397 F77_CHAR_ARG_LEN (1)
398 F77_CHAR_ARG_LEN (1)));
399
400#if defined (DEBUG)
401 if (comp_z == 'V')
402 octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl;
403#endif
404 }
405#endif
406
407 char qz_job = (nargout < 2 ? 'E' : 'S');
408
409 if (complex_case)
410 {
411 // Complex case.
412
413 // The QR decomposition of cbb.
414 math::qr<ComplexMatrix> cbqr (cbb);
415 // The R matrix of QR decomposition for cbb.
416 cbb = cbqr.R ();
417 // (Q*)caa for following work.
418 caa = (cbqr.Q ().hermitian ()) * caa;
419 CQ = CQ * cbqr.Q ();
420
421 F77_XFCN (zgghrd, ZGGHRD,
422 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
423 F77_CONST_CHAR_ARG2 (&comp_z, 1),
424 nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()),
425 nn, F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
426 F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
427 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
428 F77_CHAR_ARG_LEN (1)
429 F77_CHAR_ARG_LEN (1)));
430
431 ComplexRowVector cwork (nn);
432
433 F77_XFCN (zhgeqz, ZHGEQZ,
434 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
435 F77_CONST_CHAR_ARG2 (&comp_q, 1),
436 F77_CONST_CHAR_ARG2 (&comp_z, 1),
437 nn, ilo, ihi,
438 F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
439 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
440 F77_DBLE_CMPLX_ARG (xalpha.fortran_vec ()),
441 F77_DBLE_CMPLX_ARG (xbeta.fortran_vec ()),
442 F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn,
443 F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn,
444 F77_DBLE_CMPLX_ARG (cwork.fortran_vec ()), nn,
445 rwork.fortran_vec (), info
446 F77_CHAR_ARG_LEN (1)
447 F77_CHAR_ARG_LEN (1)
448 F77_CHAR_ARG_LEN (1)));
449
450 if (comp_q == 'V')
451 {
452 // Left eigenvector.
453 F77_XFCN (zggbak, ZGGBAK,
454 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
455 F77_CONST_CHAR_ARG2 ("L", 1),
456 nn, ilo, ihi, lscale.data (), rscale.data (),
457 nn, F77_DBLE_CMPLX_ARG (CQ.fortran_vec ()), nn, info
458 F77_CHAR_ARG_LEN (1)
459 F77_CHAR_ARG_LEN (1)));
460 }
461
462 if (comp_z == 'V')
463 {
464 // Right eigenvector.
465 F77_XFCN (zggbak, ZGGBAK,
466 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
467 F77_CONST_CHAR_ARG2 ("R", 1),
468 nn, ilo, ihi, lscale.data (), rscale.data (),
469 nn, F77_DBLE_CMPLX_ARG (CZ.fortran_vec ()), nn, info
470 F77_CHAR_ARG_LEN (1)
471 F77_CHAR_ARG_LEN (1)));
472 }
473
474 }
475 else
476 {
477#if defined (DEBUG)
478 octave_stdout << "qz: performing qr decomposition of bb" << std::endl;
479#endif
480
481 // Compute the QR factorization of bb.
482 math::qr<Matrix> bqr (bb);
483
484#if defined (DEBUG)
485 octave_stdout << "qz: qr (bb) done; now performing qz decomposition"
486 << std::endl;
487#endif
488
489 bb = bqr.R ();
490
491#if defined (DEBUG)
492 octave_stdout << "qz: extracted bb" << std::endl;
493#endif
494
495 aa = (bqr.Q ()).transpose () * aa;
496
497#if defined (DEBUG)
498 octave_stdout << "qz: updated aa " << std::endl;
499 octave_stdout << "bqr.Q () =\n" << bqr.Q () << std::endl;
500
501 if (comp_q == 'V')
502 octave_stdout << "QQ =" << QQ << std::endl;
503#endif
504
505 if (comp_q == 'V')
506 QQ = QQ * bqr.Q ();
507
508#if defined (DEBUG)
509 octave_stdout << "qz: precursors done..." << std::endl;
510#endif
511
512#if defined (DEBUG)
513 octave_stdout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z
514 << std::endl;
515#endif
516
517 // Reduce to generalized Hessenberg form.
518 F77_XFCN (dgghrd, DGGHRD,
519 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
520 F77_CONST_CHAR_ARG2 (&comp_z, 1),
521 nn, ilo, ihi, aa.fortran_vec (),
522 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
523 ZZ.fortran_vec (), nn, info
524 F77_CHAR_ARG_LEN (1)
525 F77_CHAR_ARG_LEN (1)));
526
527 // Check if just computing generalized eigenvalues,
528 // or if we're actually computing the decomposition.
529
530 // Reduce to generalized Schur form.
531 F77_XFCN (dhgeqz, DHGEQZ,
532 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
533 F77_CONST_CHAR_ARG2 (&comp_q, 1),
534 F77_CONST_CHAR_ARG2 (&comp_z, 1),
535 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
536 nn, alphar.fortran_vec (), alphai.fortran_vec (),
537 betar.fortran_vec (), QQ.fortran_vec (), nn,
538 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
539 F77_CHAR_ARG_LEN (1)
540 F77_CHAR_ARG_LEN (1)
541 F77_CHAR_ARG_LEN (1)));
542
543 if (comp_q == 'V')
544 {
545 F77_XFCN (dggbak, DGGBAK,
546 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
547 F77_CONST_CHAR_ARG2 ("L", 1),
548 nn, ilo, ihi, lscale.data (), rscale.data (),
549 nn, QQ.fortran_vec (), nn, info
550 F77_CHAR_ARG_LEN (1)
551 F77_CHAR_ARG_LEN (1)));
552
553#if defined (DEBUG)
554 if (comp_q == 'V')
555 octave_stdout << "qz: balancing done; QQ=\n" << QQ << std::endl;
556#endif
557 }
558
559 // then right
560 if (comp_z == 'V')
561 {
562 F77_XFCN (dggbak, DGGBAK,
563 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
564 F77_CONST_CHAR_ARG2 ("R", 1),
565 nn, ilo, ihi, lscale.data (), rscale.data (),
566 nn, ZZ.fortran_vec (), nn, info
567 F77_CHAR_ARG_LEN (1)
568 F77_CHAR_ARG_LEN (1)));
569
570#if defined (DEBUG)
571 if (comp_z == 'V')
572 octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl;
573#endif
574 }
575
576 }
577
578 // Order the QZ decomposition?
579 if (ord_job != 'N')
580 {
581 if (complex_case)
582 // Probably not needed, but better be safe.
583 error ("qz: cannot re-order complex QZ decomposition");
584
585#if defined (DEBUG_SORT)
586 octave_stdout << "qz: ordering eigenvalues: ord_job = "
587 << ord_job << std::endl;
588#endif
589
590 Array<F77_LOGICAL> select (dim_vector (nn, 1));
591
592 for (int j = 0; j < nn; j++)
593 {
594 switch (ord_job)
595 {
596 case 'S':
597 select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) < betar(j)*betar(j);
598 break;
599
600 case 'B':
601 select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) >= betar(j)*betar(j);
602 break;
603
604 case '+':
605 select(j) = alphar(j) * betar(j) >= 0;
606 break;
607
608 case '-':
609 select(j) = alphar(j) * betar(j) < 0;
610 break;
611
612 default:
613 // Invalid order option
614 // (should never happen since options were checked at the top).
615 panic_impossible ();
616 break;
617 }
618 }
619
620 F77_LOGICAL wantq = 0, wantz = (comp_z == 'V');
621 F77_INT ijob = 0, mm, lrwork3 = 4*nn+16, liwork = nn;
622 F77_DBLE pl, pr;
623 RowVector rwork3(lrwork3);
624 Array<F77_INT> iwork (dim_vector (liwork, 1));
625
626 F77_XFCN (dtgsen, DTGSEN,
627 (ijob, wantq, wantz,
628 select.fortran_vec (), nn,
629 aa.fortran_vec (), nn,
630 bb.fortran_vec (), nn,
631 alphar.fortran_vec (),
632 alphai.fortran_vec (),
633 betar.fortran_vec (),
634 nullptr, nn,
635 ZZ.fortran_vec (), nn,
636 mm,
637 pl, pr,
638 nullptr,
639 rwork3.fortran_vec (), lrwork3,
640 iwork.fortran_vec (), liwork,
641 info));
642
643#if defined (DEBUG_SORT)
644 octave_stdout << "qz: back from dtgsen: aa =\n";
645 octave_print_internal (octave_stdout, aa);
646 octave_stdout << "\nbb =\n";
647 octave_print_internal (octave_stdout, bb);
648 if (comp_z == 'V')
649 {
650 octave_stdout << "\nZZ =\n";
651 octave_print_internal (octave_stdout, ZZ);
652 }
653 octave_stdout << "\nqz: info=" << info;
654 octave_stdout << "\nalphar =\n";
655 octave_print_internal (octave_stdout, Matrix (alphar));
656 octave_stdout << "\nalphai =\n";
657 octave_print_internal (octave_stdout, Matrix (alphai));
658 octave_stdout << "\nbeta =\n";
659 octave_print_internal (octave_stdout, Matrix (betar));
660 octave_stdout << std::endl;
661#endif
662 }
663
664 // Compute the generalized eigenvalues as well?
665 ComplexColumnVector gev;
666
667 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
668 {
669 if (complex_case)
670 {
671 ComplexColumnVector tmp (nn);
672
673 for (F77_INT i = 0; i < nn; i++)
674 tmp(i) = xalpha(i) / xbeta(i);
675
676 gev = tmp;
677 }
678 else
679 {
680#if defined (DEBUG)
681 octave_stdout << "qz: computing generalized eigenvalues" << std::endl;
682#endif
683
684 // Return finite generalized eigenvalues.
685 ComplexColumnVector tmp (nn);
686
687 for (F77_INT i = 0; i < nn; i++)
688 {
689 if (betar(i) != 0)
690 tmp(i) = Complex (alphar(i), alphai(i)) / betar(i);
691 else
692 tmp(i) = numeric_limits<double>::Inf ();
693 }
694
695 gev = tmp;
696 }
697 }
698
699 // Right, left eigenvector matrices.
700 if (nargout >= 5)
701 {
702 // Which side to compute?
703 char side = (nargout == 5 ? 'R' : 'B');
704 // Compute all of them and backtransform
705 char howmany = 'B';
706 // Dummy pointer; select is not used.
707 F77_INT *select = nullptr;
708
709 if (complex_case)
710 {
711 CVL = CQ;
712 CVR = CZ;
713 ComplexRowVector cwork2 (2 * nn);
714 RowVector rwork2 (8 * nn);
715 F77_INT m;
716
717 F77_XFCN (ztgevc, ZTGEVC,
718 (F77_CONST_CHAR_ARG2 (&side, 1),
719 F77_CONST_CHAR_ARG2 (&howmany, 1),
720 select, nn, F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
721 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()),
722 nn, F77_DBLE_CMPLX_ARG (CVL.fortran_vec ()), nn,
723 F77_DBLE_CMPLX_ARG (CVR.fortran_vec ()), nn, nn,
724 m, F77_DBLE_CMPLX_ARG (cwork2.fortran_vec ()),
725 rwork2.fortran_vec (), info
726 F77_CHAR_ARG_LEN (1)
727 F77_CHAR_ARG_LEN (1)));
728 }
729 else
730 {
731#if defined (DEBUG)
732 octave_stdout << "qz: computing generalized eigenvectors" << std::endl;
733#endif
734
735 VL = QQ;
736 VR = ZZ;
737 F77_INT m;
738
739 F77_XFCN (dtgevc, DTGEVC,
740 (F77_CONST_CHAR_ARG2 (&side, 1),
741 F77_CONST_CHAR_ARG2 (&howmany, 1),
742 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
743 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
744 m, work.fortran_vec (), info
745 F77_CHAR_ARG_LEN (1)
746 F77_CHAR_ARG_LEN (1)));
747
748 // Now construct the complex form of VV, WW.
749 F77_INT j = 0;
750
751 while (j < nn)
752 {
753 octave_quit ();
754
755 // See if real or complex eigenvalue.
756
757 // Column increment; assume complex eigenvalue.
758 int cinc = 2;
759
760 if (j == (nn-1))
761 // Single column.
762 cinc = 1;
763 else if (aa(j+1, j) == 0)
764 cinc = 1;
765
766 // Now copy the eigenvector (s) to CVR, CVL.
767 if (cinc == 1)
768 {
769 for (F77_INT i = 0; i < nn; i++)
770 CVR(i, j) = VR(i, j);
771
772 if (side == 'B')
773 for (F77_INT i = 0; i < nn; i++)
774 CVL(i, j) = VL(i, j);
775 }
776 else
777 {
778 // Double column; complex vector.
779
780 for (F77_INT i = 0; i < nn; i++)
781 {
782 CVR(i, j) = Complex (VR(i, j), VR(i, j+1));
783 CVR(i, j+1) = Complex (VR(i, j), -VR(i, j+1));
784 }
785
786 if (side == 'B')
787 for (F77_INT i = 0; i < nn; i++)
788 {
789 CVL(i, j) = Complex (VL(i, j), VL(i, j+1));
790 CVL(i, j+1) = Complex (VL(i, j), -VL(i, j+1));
791 }
792 }
793
794 // Advance to next eigenvectors (if any).
795 j += cinc;
796 }
797 }
798 }
799
800 octave_value_list retval (nargout);
801
802 switch (nargout)
803 {
804 case 7:
805 retval(6) = gev;
806 OCTAVE_FALLTHROUGH;
807
808 case 6:
809 // Return left eigenvectors.
810 retval(5) = CVL;
811 OCTAVE_FALLTHROUGH;
812
813 case 5:
814 // Return right eigenvectors.
815 retval(4) = CVR;
816 OCTAVE_FALLTHROUGH;
817
818 case 4:
819 if (nargin == 3)
820 {
821#if defined (DEBUG)
822 octave_stdout << "qz: sort: retval(3) = gev =\n";
823 octave_print_internal (octave_stdout, ComplexMatrix (gev));
824 octave_stdout << std::endl;
825#endif
826 retval(3) = gev;
827 }
828 else
829 {
830 if (complex_case)
831 retval(3) = CZ;
832 else
833 retval(3) = ZZ;
834 }
835 OCTAVE_FALLTHROUGH;
836
837 case 3:
838 if (nargin == 3)
839 {
840 if (complex_case)
841 retval(2) = CZ;
842 else
843 retval(2) = ZZ;
844 }
845 else
846 {
847 if (complex_case)
848 retval(2) = CQ.hermitian ();
849 else
850 retval(2) = QQ.transpose ();
851 }
852 OCTAVE_FALLTHROUGH;
853
854 case 2:
855 {
856 if (complex_case)
857 {
858#if defined (DEBUG)
859 octave_stdout << "qz: retval(1) = cbb =\n";
860 octave_print_internal (octave_stdout, cbb);
861 octave_stdout << "\nqz: retval(0) = caa =\n";
862 octave_print_internal (octave_stdout, caa);
863 octave_stdout << std::endl;
864#endif
865 retval(1) = cbb;
866 retval(0) = caa;
867 }
868 else
869 {
870#if defined (DEBUG)
871 octave_stdout << "qz: retval(1) = bb =\n";
872 octave_print_internal (octave_stdout, bb);
873 octave_stdout << "\nqz: retval(0) = aa =\n";
874 octave_print_internal (octave_stdout, aa);
875 octave_stdout << std::endl;
876#endif
877 retval(1) = bb;
878 retval(0) = aa;
879 }
880 }
881 break;
882
883 case 1:
884 case 0:
885#if defined (DEBUG)
886 octave_stdout << "qz: retval(0) = gev = " << gev << std::endl;
887#endif
888 retval(0) = gev;
889 break;
890
891 default:
892 error ("qz: too many return arguments");
893 break;
894 }
895
896#if defined (DEBUG)
897 octave_stdout << "qz: exiting (at long last)" << std::endl;
898#endif
899
900 return retval;
901}
902
903/*
904%!test
905%! A = [1 2; 0 3];
906%! B = [1 0; 0 0];
907%! C = [0 1; 0 0];
908%! assert (qz (A,B), [1; Inf]);
909%! assert (qz (A,C), [Inf; Inf]);
910
911## Example 7.7.3 in Golub & Van Loan
912%!test
913%! a = [ 10 1 2;
914%! 1 2 -1;
915%! 1 1 2 ];
916%! b = reshape (1:9, 3,3);
917%! [aa, bb, q, z, v, w, lambda] = qz (a, b);
918%! assert (q * a * z, aa, norm (aa) * 1e-14);
919%! assert (q * b * z, bb, norm (bb) * 1e-14);
920%! is_finite = abs (lambda) < 1 / eps (max (a(:)));
921%! observed = (b * v * diag (lambda))(:,is_finite);
922%! assert (observed, (a*v)(:,is_finite), norm (observed) * 1e-14);
923%! observed = (diag (lambda) * w' * b)(is_finite,:);
924%! assert (observed, (w'*a)(is_finite,:), norm (observed) * 1e-13);
925
926%!test
927%! A = [0, 0, -1, 0; 1, 0, 0, 0; -1, 0, -2, -1; 0, -1, 1, 0];
928%! B = [0, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1];
929%! [AA, BB, Q, Z1] = qz (A, B);
930%! [AA, BB, Z2] = qz (A, B, "-");
931%! assert (Z1, Z2);
932
933%!test
934%! A = [ -1.03428 0.24929 0.43205 -0.12860;
935%! 1.16228 0.27870 2.12954 0.69250;
936%! -0.51524 -0.34939 -0.77820 2.13721;
937%! -1.32941 2.11870 0.72005 1.00835 ];
938%! B = [ 1.407302 -0.632956 -0.360628 0.068534;
939%! 0.149898 0.298248 0.991777 0.023652;
940%! 0.169281 -0.405205 -1.775834 1.511730;
941%! 0.717770 1.291390 -1.766607 -0.531352 ];
942%! [AA, BB, Z, lambda] = qz (A, B, "+");
943%! assert (all (real (lambda(1:3)) >= 0));
944%! assert (real (lambda(4) < 0));
945%! [AA, BB, Z, lambda] = qz (A, B, "-");
946%! assert (real (lambda(1) < 0));
947%! assert (all (real (lambda(2:4)) >= 0));
948%! [AA, BB, Z, lambda] = qz (A, B, "B");
949%! assert (all (abs (lambda(1:3)) >= 1));
950%! assert (abs (lambda(4) < 1));
951%! [AA, BB, Z, lambda] = qz (A, B, "S");
952%! assert (abs (lambda(1) < 1));
953%! assert (all (abs (lambda(2:4)) >= 1));
954*/
955
956OCTAVE_NAMESPACE_END
nn
static F77_INT nn
Definition: DASPK.cc:66
Inf
#define Inf
Definition: Faddeeva.cc:260
Array
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:129
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
Matrix
Definition: dMatrix.h:42
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_value_list
Definition: ovl.h:44
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
error
void error(const char *fmt,...)
Definition: error.cc:980
panic_impossible
#define panic_impossible()
Definition: error.h:411
err_square_matrix_required
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:122
warn_empty_arg
void warn_empty_arg(const char *name)
Definition: errwarn.cc:330
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
F77_DBLE
double F77_DBLE
Definition: f77-fcn.h:302
F77_LOGICAL
octave_f77_int_type F77_LOGICAL
Definition: f77-fcn.h:308
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
VR
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VR
Definition: lo-lapack-proto.h:1794
CZ
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX * CZ
Definition: lo-lapack-proto.h:1401
CQ
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX * CQ
Definition: lo-lapack-proto.h:1400
VL
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VL
Definition: lo-lapack-proto.h:1793
ComplexMatrix
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
octave::err_nonconformant
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-errwarn.cc:71
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
octave_stdout
#define octave_stdout
Definition: pager.h:314
octave_print_internal
void octave_print_internal(std::ostream &os, const float_display_format &fmt, bool d, bool pr_as_read_syntax)
Definition: pr-output.cc:1762
transpose
static void transpose(octave_idx_type N, const octave_idx_type *ridx, const octave_idx_type *cidx, octave_idx_type *ridx2, octave_idx_type *cidx2)
Definition: symrcm.cc:391
xdlamch
subroutine xdlamch(cmach, retval)
Definition: xdlamch.f:2
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg