136 int nargin = args.length ();
140 <<
", nargout = " << nargout << std::endl;
143 if (nargin < 2 || nargin > 3 || nargout > 6)
146 bool complex_case =
true;
150 std::string opt = args(2).xstring_value (
"qz: OPT must be a string");
154 if (args(0).iscomplex () || args(1).iscomplex ())
156 warning (
"qz: ignoring 'real' option with complex matrices");
160 complex_case =
false;
162 else if (opt ==
"complex")
165 error (
"qz: OPT must be 'real' or 'complex'");
173 F77_INT nn = to_f77_int (args(0).rows ());
174 F77_INT nc = to_f77_int (args(0).columns ());
177 octave_stdout <<
"Matrix A dimensions: (" << nn <<
',' << nc <<
')'
188 if (args(0).iscomplex ())
189 caa = args(0).complex_matrix_value ();
191 aa = args(0).matrix_value ();
198 F77_INT b_nr = to_f77_int (args(1).rows ());
199 F77_INT b_nc = to_f77_int (args(1).columns ());
201 if (nn != b_nc || nn != b_nr)
207 if (args(1).iscomplex ())
208 cbb = args(1).complex_matrix_value ();
210 bb = args(1).matrix_value ();
215 Matrix QQ (nn, nn), ZZ (nn, nn),
VR (nn, nn),
VL (nn, nn);
216 RowVector alphar (nn), alphai (nn), betar (nn);
220 char comp_q = (nargout >= 3 ?
'V' :
'N');
221 char comp_z = (nargout >= 4 ?
'V' :
'N');
224 if (comp_q ==
'V' || comp_z ==
'V')
226 double *QQptr = QQ.rwdata ();
227 double *ZZptr = ZZ.rwdata ();
228 std::fill_n (QQptr, QQ.numel (), 0.0);
229 std::fill_n (ZZptr, ZZ.numel (), 0.0);
230 for (
F77_INT i = 0; i < nn; i++)
238 const char bal_job =
'P';
239 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
245 octave_stdout <<
"qz: performing balancing; CQ =\n" <<
CQ << std::endl;
247 if (args(0).isreal ())
250 if (args(1).isreal ())
260 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
263 nn, ilo, ihi, lscale.rwdata (),
264 rscale.rwdata (), work.rwdata (), info
265 F77_CHAR_ARG_LEN (1)));
271 octave_stdout <<
"qz: performing balancing; QQ =\n" << QQ << std::endl;
275 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
277 nn, ilo, ihi, lscale.rwdata (),
278 rscale.rwdata (), work.rwdata (), info
279 F77_CHAR_ARG_LEN (1)));
284 char qz_job = (nargout < 2 ?
'E' :
'S');
291 math::qr<ComplexMatrix> cbqr (cbb);
295 caa = (cbqr.Q ().hermitian ()) * caa;
299 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
300 F77_CONST_CHAR_ARG2 (&comp_z, 1),
306 F77_CHAR_ARG_LEN (1)));
311 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
312 F77_CONST_CHAR_ARG2 (&comp_q, 1),
313 F77_CONST_CHAR_ARG2 (&comp_z, 1),
325 F77_CHAR_ARG_LEN (1)));
331 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
332 F77_CONST_CHAR_ARG2 (
"L", 1),
333 nn, ilo, ihi, lscale.data (), rscale.data (),
336 F77_CHAR_ARG_LEN (1)));
343 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
344 F77_CONST_CHAR_ARG2 (
"R", 1),
345 nn, ilo, ihi, lscale.data (), rscale.data (),
348 F77_CHAR_ARG_LEN (1)));
355 octave_stdout <<
"qz: performing qr decomposition of bb" << std::endl;
359 math::qr<Matrix> bqr (bb);
362 octave_stdout <<
"qz: qr (bb) done; now performing qz decomposition"
372 aa = (bqr.Q ()).transpose () * aa;
390 octave_stdout <<
"qz: comp_q = " << comp_q <<
", comp_z = " << comp_z
396 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
397 F77_CONST_CHAR_ARG2 (&comp_z, 1),
398 nn, ilo, ihi, aa.
rwdata (),
399 nn, bb.
rwdata (), nn, QQ.rwdata (), nn,
400 ZZ.rwdata (), nn, info
402 F77_CHAR_ARG_LEN (1)));
409 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
410 F77_CONST_CHAR_ARG2 (&comp_q, 1),
411 F77_CONST_CHAR_ARG2 (&comp_z, 1),
413 nn, alphar.rwdata (), alphai.rwdata (),
414 betar.
rwdata (), QQ.rwdata (), nn,
415 ZZ.rwdata (), nn, work.rwdata (), nn, info
418 F77_CHAR_ARG_LEN (1)));
423 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
424 F77_CONST_CHAR_ARG2 (
"L", 1),
425 nn, ilo, ihi, lscale.data (), rscale.data (),
426 nn, QQ.rwdata (), nn, info
428 F77_CHAR_ARG_LEN (1)));
432 octave_stdout <<
"qz: balancing done; QQ=\n" << QQ << std::endl;
440 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
441 F77_CONST_CHAR_ARG2 (
"R", 1),
442 nn, ilo, ihi, lscale.data (), rscale.data (),
443 nn, ZZ.rwdata (), nn, info
445 F77_CHAR_ARG_LEN (1)));
449 octave_stdout <<
"qz: balancing done; ZZ=\n" << ZZ << std::endl;
459 char side = (nargout == 5 ?
'R' :
'B');
474 (F77_CONST_CHAR_ARG2 (&side, 1),
475 F77_CONST_CHAR_ARG2 (&howmany, 1),
483 F77_CHAR_ARG_LEN (1)));
488 octave_stdout <<
"qz: computing generalized eigenvectors" << std::endl;
496 (F77_CONST_CHAR_ARG2 (&side, 1),
497 F77_CONST_CHAR_ARG2 (&howmany, 1),
499 nn,
VL.rwdata (), nn,
VR.rwdata (), nn, nn,
500 m, work.rwdata (), info
502 F77_CHAR_ARG_LEN (1)));
535 retval(2) =
CQ.hermitian ();
537 retval(2) = QQ.transpose ();
578 "qz: requesting a single output argument no longer returns eigenvalues since version 9");
582 if (nargin == 2 && args(0).isreal () && args(1).isreal ()
583 && retval(0).iscomplex ())
586 "qz: returns the complex QZ by default on real matrices since version 9");