Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/qz.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 2c279308f6ab
children d0ce5e973937
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
57 #include "symtab.h" 57 #include "symtab.h"
58 #include "utils.h" 58 #include "utils.h"
59 #include "variables.h" 59 #include "variables.h"
60 60
61 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE, const double& ALPHA, 61 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE, const double& ALPHA,
62 const double& BETA, const double& S, 62 const double& BETA, const double& S,
63 const double& P); 63 const double& P);
64 64
65 extern "C" 65 extern "C"
66 { 66 {
67 F77_RET_T 67 F77_RET_T
68 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, 68 F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL,
69 const octave_idx_type& N, double* A, const octave_idx_type& LDA, 69 const octave_idx_type& N, double* A, const octave_idx_type& LDA,
70 double* B, const octave_idx_type& LDB, octave_idx_type& ILO, 70 double* B, const octave_idx_type& LDB, octave_idx_type& ILO,
71 octave_idx_type& IHI, double* LSCALE, double* RSCALE, 71 octave_idx_type& IHI, double* LSCALE, double* RSCALE,
72 double* WORK, octave_idx_type& INFO 72 double* WORK, octave_idx_type& INFO
73 F77_CHAR_ARG_LEN_DECL); 73 F77_CHAR_ARG_LEN_DECL);
74 74
75 F77_RET_T 75 F77_RET_T
76 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, 76 F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
77 F77_CONST_CHAR_ARG_DECL, 77 F77_CONST_CHAR_ARG_DECL,
78 const octave_idx_type& N, const octave_idx_type& ILO, 78 const octave_idx_type& N, const octave_idx_type& ILO,
79 const octave_idx_type& IHI, const double* LSCALE, 79 const octave_idx_type& IHI, const double* LSCALE,
80 const double* RSCALE, octave_idx_type& M, double* V, 80 const double* RSCALE, octave_idx_type& M, double* V,
81 const octave_idx_type& LDV, octave_idx_type& INFO 81 const octave_idx_type& LDV, octave_idx_type& INFO
82 F77_CHAR_ARG_LEN_DECL 82 F77_CHAR_ARG_LEN_DECL
83 F77_CHAR_ARG_LEN_DECL); 83 F77_CHAR_ARG_LEN_DECL);
84 84
85 F77_RET_T 85 F77_RET_T
86 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL, 86 F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL,
87 F77_CONST_CHAR_ARG_DECL, 87 F77_CONST_CHAR_ARG_DECL,
88 const octave_idx_type& N, const octave_idx_type& ILO, 88 const octave_idx_type& N, const octave_idx_type& ILO,
89 const octave_idx_type& IHI, double* A, 89 const octave_idx_type& IHI, double* A,
90 const octave_idx_type& LDA, double* B, 90 const octave_idx_type& LDA, double* B,
91 const octave_idx_type& LDB, double* Q, 91 const octave_idx_type& LDB, double* Q,
92 const octave_idx_type& LDQ, double* Z, 92 const octave_idx_type& LDQ, double* Z,
93 const octave_idx_type& LDZ, octave_idx_type& INFO 93 const octave_idx_type& LDZ, octave_idx_type& INFO
94 F77_CHAR_ARG_LEN_DECL 94 F77_CHAR_ARG_LEN_DECL
95 F77_CHAR_ARG_LEN_DECL); 95 F77_CHAR_ARG_LEN_DECL);
96 96
97 F77_RET_T 97 F77_RET_T
98 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL, 98 F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL,
99 F77_CONST_CHAR_ARG_DECL, 99 F77_CONST_CHAR_ARG_DECL,
100 F77_CONST_CHAR_ARG_DECL, 100 F77_CONST_CHAR_ARG_DECL,
101 const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI, 101 const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI,
102 double* A, const octave_idx_type& LDA, double* B, 102 double* A, const octave_idx_type& LDA, double* B,
103 const octave_idx_type& LDB, double* ALPHAR, 103 const octave_idx_type& LDB, double* ALPHAR,
104 double* ALPHAI, double* BETA, double* Q, 104 double* ALPHAI, double* BETA, double* Q,
105 const octave_idx_type& LDQ, double* Z, 105 const octave_idx_type& LDQ, double* Z,
106 const octave_idx_type& LDZ, double* WORK, 106 const octave_idx_type& LDZ, double* WORK,
107 const octave_idx_type& LWORK, octave_idx_type& INFO 107 const octave_idx_type& LWORK, octave_idx_type& INFO
108 F77_CHAR_ARG_LEN_DECL 108 F77_CHAR_ARG_LEN_DECL
109 F77_CHAR_ARG_LEN_DECL 109 F77_CHAR_ARG_LEN_DECL
110 F77_CHAR_ARG_LEN_DECL); 110 F77_CHAR_ARG_LEN_DECL);
111 111
112 F77_RET_T 112 F77_RET_T
113 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, const double* B, 113 F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, const double* B,
114 const octave_idx_type& LDB, const double& SAFMIN, 114 const octave_idx_type& LDB, const double& SAFMIN,
115 double& SCALE1, double& SCALE2, 115 double& SCALE1, double& SCALE2,
116 double& WR1, double& WR2, double& WI); 116 double& WR1, double& WR2, double& WI);
117 117
118 // Van Dooren's code (netlib.org: toms/590) for reordering 118 // Van Dooren's code (netlib.org: toms/590) for reordering
119 // GEP. Only processes Z, not Q. 119 // GEP. Only processes Z, not Q.
120 F77_RET_T 120 F77_RET_T
121 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX, const octave_idx_type& N, double* A, 121 F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX, const octave_idx_type& N, double* A,
122 double* B, double* Z, sort_function, 122 double* B, double* Z, sort_function,
123 const double& EPS, octave_idx_type& NDIM, octave_idx_type& FAIL, 123 const double& EPS, octave_idx_type& NDIM, octave_idx_type& FAIL,
124 octave_idx_type* IND); 124 octave_idx_type* IND);
125 125
126 // documentation for DTGEVC incorrectly states that VR, VL are 126 // documentation for DTGEVC incorrectly states that VR, VL are
127 // complex*16; they are declared in DTGEVC as double precision 127 // complex*16; they are declared in DTGEVC as double precision
128 // (probably a cut and paste problem fro ZTGEVC) 128 // (probably a cut and paste problem fro ZTGEVC)
129 F77_RET_T 129 F77_RET_T
130 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL, 130 F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL,
131 F77_CONST_CHAR_ARG_DECL, 131 F77_CONST_CHAR_ARG_DECL,
132 octave_idx_type* SELECT, const octave_idx_type& N, double* A, 132 octave_idx_type* SELECT, const octave_idx_type& N, double* A,
133 const octave_idx_type& LDA, double* B, 133 const octave_idx_type& LDA, double* B,
134 const octave_idx_type& LDB, double* VL, 134 const octave_idx_type& LDB, double* VL,
135 const octave_idx_type& LDVL, double* VR, 135 const octave_idx_type& LDVL, double* VR,
136 const octave_idx_type& LDVR, const octave_idx_type& MM, 136 const octave_idx_type& LDVR, const octave_idx_type& MM,
137 octave_idx_type& M, double* WORK, octave_idx_type& INFO 137 octave_idx_type& M, double* WORK, octave_idx_type& INFO
138 F77_CHAR_ARG_LEN_DECL 138 F77_CHAR_ARG_LEN_DECL
139 F77_CHAR_ARG_LEN_DECL); 139 F77_CHAR_ARG_LEN_DECL);
140 140
141 F77_RET_T 141 F77_RET_T
142 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL, 142 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL,
143 double& retval 143 double& retval
144 F77_CHAR_ARG_LEN_DECL); 144 F77_CHAR_ARG_LEN_DECL);
145 145
146 F77_RET_T 146 F77_RET_T
147 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL, 147 F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL,
148 const octave_idx_type&, const octave_idx_type&, const double*, 148 const octave_idx_type&, const octave_idx_type&, const double*,
149 const octave_idx_type&, double*, double& 149 const octave_idx_type&, double*, double&
150 F77_CHAR_ARG_LEN_DECL); 150 F77_CHAR_ARG_LEN_DECL);
151 } 151 }
152 152
153 // fcrhp, fin, fout, folhp: 153 // fcrhp, fin, fout, folhp:
154 // routines for ordering of generalized eigenvalues 154 // routines for ordering of generalized eigenvalues
155 // return 1 if test is passed, 0 otherwise 155 // return 1 if test is passed, 0 otherwise
317 else 317 else
318 { 318 {
319 std::string tmp = args(2).string_value (); 319 std::string tmp = args(2).string_value ();
320 320
321 if (! tmp.empty ()) 321 if (! tmp.empty ())
322 ord_job = tmp[0]; 322 ord_job = tmp[0];
323 323
324 if (! (ord_job == 'N' || ord_job == 'n' 324 if (! (ord_job == 'N' || ord_job == 'n'
325 || ord_job == 'S' || ord_job == 's' 325 || ord_job == 'S' || ord_job == 's'
326 || ord_job == 'B' || ord_job == 'b' 326 || ord_job == 'B' || ord_job == 'b'
327 || ord_job == '+' || ord_job == '-')) 327 || ord_job == '+' || ord_job == '-'))
328 { 328 {
329 error ("qz: invalid order option"); 329 error ("qz: invalid order option");
330 return retval; 330 return retval;
331 } 331 }
332 332
333 // overflow constant required by dlag2 333 // overflow constant required by dlag2
334 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1), 334 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
335 safmin 335 safmin
336 F77_CHAR_ARG_LEN (1)); 336 F77_CHAR_ARG_LEN (1));
337 337
338 #ifdef DEBUG_EIG 338 #ifdef DEBUG_EIG
339 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific) 339 std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific)
340 << safmin << std::endl; 340 << safmin << std::endl;
341 #endif 341 #endif
342 342
343 // some machines (e.g., DEC alpha) get safmin = 0; 343 // some machines (e.g., DEC alpha) get safmin = 0;
344 // for these, use eps instead to avoid problems in dlag2 344 // for these, use eps instead to avoid problems in dlag2
345 if (safmin == 0) 345 if (safmin == 0)
346 { 346 {
347 #ifdef DEBUG_EIG 347 #ifdef DEBUG_EIG
348 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl; 348 std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
349 #endif 349 #endif
350 350
351 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1), 351 F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
352 safmin 352 safmin
353 F77_CHAR_ARG_LEN (1)); 353 F77_CHAR_ARG_LEN (1));
354 354
355 #ifdef DEBUG_EIG 355 #ifdef DEBUG_EIG
356 std::cout << "qz: safmin set to " << setiosflags (std::ios::scientific) 356 std::cout << "qz: safmin set to " << setiosflags (std::ios::scientific)
357 << safmin << std::endl; 357 << safmin << std::endl;
358 #endif 358 #endif
359 } 359 }
360 } 360 }
361 361
362 #ifdef DEBUG 362 #ifdef DEBUG
363 std::cout << "qz: check argument 1" << std::endl; 363 std::cout << "qz: check argument 1" << std::endl;
364 #endif 364 #endif
446 446
447 // initialize Q, Z to identity if we need either of them 447 // initialize Q, Z to identity if we need either of them
448 if (compq == 'V' || compz == 'V') 448 if (compq == 'V' || compz == 'V')
449 for (octave_idx_type ii = 0; ii < nn; ii++) 449 for (octave_idx_type ii = 0; ii < nn; ii++)
450 for (octave_idx_type jj = 0; jj < nn; jj++) 450 for (octave_idx_type jj = 0; jj < nn; jj++)
451 { 451 {
452 OCTAVE_QUIT; 452 OCTAVE_QUIT;
453 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); 453 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
454 } 454 }
455 455
456 // always perform permutation balancing 456 // always perform permutation balancing
457 const char bal_job = 'P'; 457 const char bal_job = 'P';
458 RowVector lscale(nn), rscale(nn), work(6*nn); 458 RowVector lscale(nn), rscale(nn), work(6*nn);
459 459
464 } 464 }
465 else 465 else
466 { 466 {
467 #ifdef DEBUG 467 #ifdef DEBUG
468 if (compq == 'V') 468 if (compq == 'V')
469 std::cout << "qz: performing balancing; QQ=" << std::endl << QQ << std::endl; 469 std::cout << "qz: performing balancing; QQ=" << std::endl << QQ << std::endl;
470 #endif 470 #endif
471 471
472 F77_XFCN (dggbal, DGGBAL, 472 F77_XFCN (dggbal, DGGBAL,
473 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 473 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
474 nn, aa.fortran_vec (), nn, bb.fortran_vec (), 474 nn, aa.fortran_vec (), nn, bb.fortran_vec (),
475 nn, ilo, ihi, lscale.fortran_vec (), 475 nn, ilo, ihi, lscale.fortran_vec (),
476 rscale.fortran_vec (), work.fortran_vec (), info 476 rscale.fortran_vec (), work.fortran_vec (), info
477 F77_CHAR_ARG_LEN (1))); 477 F77_CHAR_ARG_LEN (1)));
478 } 478 }
479 479
480 // Since we just want the balancing matrices, we can use dggbal 480 // Since we just want the balancing matrices, we can use dggbal
481 // for both the real and complex cases; 481 // for both the real and complex cases;
482 // left first 482 // left first
483 483
484 if (compq == 'V') 484 if (compq == 'V')
485 { 485 {
486 F77_XFCN (dggbak, DGGBAK, 486 F77_XFCN (dggbak, DGGBAK,
487 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 487 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
488 F77_CONST_CHAR_ARG2 ("L", 1), 488 F77_CONST_CHAR_ARG2 ("L", 1),
489 nn, ilo, ihi, lscale.data (), rscale.data (), 489 nn, ilo, ihi, lscale.data (), rscale.data (),
490 nn, QQ.fortran_vec (), nn, info 490 nn, QQ.fortran_vec (), nn, info
491 F77_CHAR_ARG_LEN (1) 491 F77_CHAR_ARG_LEN (1)
492 F77_CHAR_ARG_LEN (1))); 492 F77_CHAR_ARG_LEN (1)));
493 493
494 #ifdef DEBUG 494 #ifdef DEBUG
495 if (compq == 'V') 495 if (compq == 'V')
496 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; 496 std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
497 #endif 497 #endif
498 } 498 }
499 499
500 // then right 500 // then right
501 if (compz == 'V') 501 if (compz == 'V')
502 { 502 {
503 F77_XFCN (dggbak, DGGBAK, 503 F77_XFCN (dggbak, DGGBAK,
504 (F77_CONST_CHAR_ARG2 (&bal_job, 1), 504 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
505 F77_CONST_CHAR_ARG2 ("R", 1), 505 F77_CONST_CHAR_ARG2 ("R", 1),
506 nn, ilo, ihi, lscale.data (), rscale.data (), 506 nn, ilo, ihi, lscale.data (), rscale.data (),
507 nn, ZZ.fortran_vec (), nn, info 507 nn, ZZ.fortran_vec (), nn, info
508 F77_CHAR_ARG_LEN (1) 508 F77_CHAR_ARG_LEN (1)
509 F77_CHAR_ARG_LEN (1))); 509 F77_CHAR_ARG_LEN (1)));
510 510
511 #ifdef DEBUG 511 #ifdef DEBUG
512 if (compz == 'V') 512 if (compz == 'V')
513 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; 513 std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
514 #endif 514 #endif
515 } 515 }
516 516
517 static char qz_job; 517 static char qz_job;
518 qz_job = (nargout < 2 ? 'E' : 'S'); 518 qz_job = (nargout < 2 ? 'E' : 'S');
519 519
520 if (complex_case) 520 if (complex_case)
521 { 521 {
522 // complex case 522 // complex case
523 if (args(0).is_real_type ()) 523 if (args(0).is_real_type ())
524 caa = ComplexMatrix (aa); 524 caa = ComplexMatrix (aa);
525 525
526 if (args(1).is_real_type ()) 526 if (args(1).is_real_type ())
527 cbb = ComplexMatrix (bb); 527 cbb = ComplexMatrix (bb);
528 528
529 if (compq == 'V') 529 if (compq == 'V')
530 CQ = ComplexMatrix (QQ); 530 CQ = ComplexMatrix (QQ);
531 531
532 if (compz == 'V') 532 if (compz == 'V')
533 CZ = ComplexMatrix (ZZ); 533 CZ = ComplexMatrix (ZZ);
534 534
535 error ("complex case not done yet"); 535 error ("complex case not done yet");
536 return retval; 536 return retval;
537 } 537 }
538 else // real matrices case 538 else // real matrices case
539 { 539 {
540 #ifdef DEBUG 540 #ifdef DEBUG
541 std::cout << "qz: peforming qr decomposition of bb" << std::endl; 541 std::cout << "qz: peforming qr decomposition of bb" << std::endl;
542 #endif 542 #endif
543 543
559 #ifdef DEBUG 559 #ifdef DEBUG
560 std::cout << "qz: updated aa " << std::endl; 560 std::cout << "qz: updated aa " << std::endl;
561 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl; 561 std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
562 562
563 if (compq == 'V') 563 if (compq == 'V')
564 std::cout << "QQ =" << QQ << std::endl; 564 std::cout << "QQ =" << QQ << std::endl;
565 #endif 565 #endif
566 566
567 if (compq == 'V') 567 if (compq == 'V')
568 QQ = QQ*bqr.Q (); 568 QQ = QQ*bqr.Q ();
569 569
570 #ifdef DEBUG 570 #ifdef DEBUG
571 std::cout << "qz: precursors done..." << std::endl; 571 std::cout << "qz: precursors done..." << std::endl;
572 #endif 572 #endif
573 573
575 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl; 575 std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl;
576 #endif 576 #endif
577 577
578 // reduce to generalized hessenberg form 578 // reduce to generalized hessenberg form
579 F77_XFCN (dgghrd, DGGHRD, 579 F77_XFCN (dgghrd, DGGHRD,
580 (F77_CONST_CHAR_ARG2 (&compq, 1), 580 (F77_CONST_CHAR_ARG2 (&compq, 1),
581 F77_CONST_CHAR_ARG2 (&compz, 1), 581 F77_CONST_CHAR_ARG2 (&compz, 1),
582 nn, ilo, ihi, aa.fortran_vec (), 582 nn, ilo, ihi, aa.fortran_vec (),
583 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn, 583 nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
584 ZZ.fortran_vec (), nn, info 584 ZZ.fortran_vec (), nn, info
585 F77_CHAR_ARG_LEN (1) 585 F77_CHAR_ARG_LEN (1)
586 F77_CHAR_ARG_LEN (1))); 586 F77_CHAR_ARG_LEN (1)));
587 587
588 // check if just computing generalized eigenvalues or if we're 588 // check if just computing generalized eigenvalues or if we're
589 // actually computing the decomposition 589 // actually computing the decomposition
590 590
591 // reduce to generalized Schur form 591 // reduce to generalized Schur form
592 F77_XFCN (dhgeqz, DHGEQZ, 592 F77_XFCN (dhgeqz, DHGEQZ,
593 (F77_CONST_CHAR_ARG2 (&qz_job, 1), 593 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
594 F77_CONST_CHAR_ARG2 (&compq, 1), 594 F77_CONST_CHAR_ARG2 (&compq, 1),
595 F77_CONST_CHAR_ARG2 (&compz, 1), 595 F77_CONST_CHAR_ARG2 (&compz, 1),
596 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (), 596 nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
597 nn, alphar.fortran_vec (), alphai.fortran_vec (), 597 nn, alphar.fortran_vec (), alphai.fortran_vec (),
598 betar.fortran_vec (), QQ.fortran_vec (), nn, 598 betar.fortran_vec (), QQ.fortran_vec (), nn,
599 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info 599 ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
600 F77_CHAR_ARG_LEN (1) 600 F77_CHAR_ARG_LEN (1)
601 F77_CHAR_ARG_LEN (1) 601 F77_CHAR_ARG_LEN (1)
602 F77_CHAR_ARG_LEN (1))); 602 F77_CHAR_ARG_LEN (1)));
603 } 603 }
604 604
605 // order the QZ decomposition? 605 // order the QZ decomposition?
606 if (! (ord_job == 'N' || ord_job == 'n')) 606 if (! (ord_job == 'N' || ord_job == 'n'))
607 { 607 {
608 if (complex_case) 608 if (complex_case)
609 { 609 {
610 // probably not needed, but better be safe 610 // probably not needed, but better be safe
611 error ("qz: cannot re-order complex qz decomposition."); 611 error ("qz: cannot re-order complex qz decomposition.");
612 return retval; 612 return retval;
613 } 613 }
614 else 614 else
615 { 615 {
616 #ifdef DEBUG_SORT 616 #ifdef DEBUG_SORT
617 std::cout << "qz: ordering eigenvalues: ord_job = " 617 std::cout << "qz: ordering eigenvalues: ord_job = "
618 << ord_job << std::endl; 618 << ord_job << std::endl;
619 #endif 619 #endif
620 620
621 // declared static to avoid vfork/long jump compiler complaints 621 // declared static to avoid vfork/long jump compiler complaints
622 static sort_function sort_test; 622 static sort_function sort_test;
623 sort_test = 0; 623 sort_test = 0;
624 624
625 switch (ord_job) 625 switch (ord_job)
626 { 626 {
627 case 'S': 627 case 'S':
628 case 's': 628 case 's':
629 sort_test = &fin; 629 sort_test = &fin;
630 break; 630 break;
631 631
632 case 'B': 632 case 'B':
633 case 'b': 633 case 'b':
634 sort_test = &fout; 634 sort_test = &fout;
635 break; 635 break;
636 636
637 case '+': 637 case '+':
638 sort_test = &fcrhp; 638 sort_test = &fcrhp;
639 break; 639 break;
640 640
641 case '-': 641 case '-':
642 sort_test = &folhp; 642 sort_test = &folhp;
643 break; 643 break;
644 644
645 default: 645 default:
646 // invalid order option (should never happen, since we 646 // invalid order option (should never happen, since we
647 // checked the options at the top). 647 // checked the options at the top).
648 panic_impossible (); 648 panic_impossible ();
649 break; 649 break;
650 } 650 }
651 651
652 octave_idx_type ndim, fail; 652 octave_idx_type ndim, fail;
653 double inf_norm; 653 double inf_norm;
654 654
655 F77_XFCN (xdlange, XDLANGE, 655 F77_XFCN (xdlange, XDLANGE,
656 (F77_CONST_CHAR_ARG2 ("I", 1), 656 (F77_CONST_CHAR_ARG2 ("I", 1),
657 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm 657 nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
658 F77_CHAR_ARG_LEN (1))); 658 F77_CHAR_ARG_LEN (1)));
659 659
660 double eps = DBL_EPSILON*inf_norm*nn; 660 double eps = DBL_EPSILON*inf_norm*nn;
661 661
662 #ifdef DEBUG_SORT 662 #ifdef DEBUG_SORT
663 std::cout << "qz: calling dsubsp: aa=" << std::endl; 663 std::cout << "qz: calling dsubsp: aa=" << std::endl;
664 octave_print_internal (std::cout, aa, 0); 664 octave_print_internal (std::cout, aa, 0);
665 std::cout << std::endl << "bb=" << std::endl; 665 std::cout << std::endl << "bb=" << std::endl;
666 octave_print_internal (std::cout, bb, 0); 666 octave_print_internal (std::cout, bb, 0);
667 if (compz == 'V') 667 if (compz == 'V')
668 { 668 {
669 std::cout << std::endl << "ZZ=" << std::endl; 669 std::cout << std::endl << "ZZ=" << std::endl;
670 octave_print_internal (std::cout, ZZ, 0); 670 octave_print_internal (std::cout, ZZ, 0);
671 } 671 }
672 std::cout << std::endl; 672 std::cout << std::endl;
673 std::cout << "alphar = " << std::endl; 673 std::cout << "alphar = " << std::endl;
674 octave_print_internal (std::cout, (Matrix) alphar, 0); 674 octave_print_internal (std::cout, (Matrix) alphar, 0);
675 std::cout << std::endl << "alphai = " << std::endl; 675 std::cout << std::endl << "alphai = " << std::endl;
676 octave_print_internal (std::cout, (Matrix) alphai, 0); 676 octave_print_internal (std::cout, (Matrix) alphai, 0);
677 std::cout << std::endl << "beta = " << std::endl; 677 std::cout << std::endl << "beta = " << std::endl;
678 octave_print_internal (std::cout, (Matrix) betar, 0); 678 octave_print_internal (std::cout, (Matrix) betar, 0);
679 std::cout << std::endl; 679 std::cout << std::endl;
680 #endif 680 #endif
681 681
682 Array<octave_idx_type> ind (nn); 682 Array<octave_idx_type> ind (nn);
683 683
684 F77_XFCN (dsubsp, DSUBSP, 684 F77_XFCN (dsubsp, DSUBSP,
685 (nn, nn, aa.fortran_vec (), bb.fortran_vec (), 685 (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
686 ZZ.fortran_vec (), sort_test, eps, ndim, fail, 686 ZZ.fortran_vec (), sort_test, eps, ndim, fail,
687 ind.fortran_vec ())); 687 ind.fortran_vec ()));
688 688
689 #ifdef DEBUG 689 #ifdef DEBUG
690 std::cout << "qz: back from dsubsp: aa=" << std::endl; 690 std::cout << "qz: back from dsubsp: aa=" << std::endl;
691 octave_print_internal (std::cout, aa, 0); 691 octave_print_internal (std::cout, aa, 0);
692 std::cout << std::endl << "bb=" << std::endl; 692 std::cout << std::endl << "bb=" << std::endl;
693 octave_print_internal (std::cout, bb, 0); 693 octave_print_internal (std::cout, bb, 0);
694 if (compz == 'V') 694 if (compz == 'V')
695 { 695 {
696 std::cout << std::endl << "ZZ=" << std::endl; 696 std::cout << std::endl << "ZZ=" << std::endl;
697 octave_print_internal (std::cout, ZZ, 0); 697 octave_print_internal (std::cout, ZZ, 0);
698 } 698 }
699 std::cout << std::endl; 699 std::cout << std::endl;
700 #endif 700 #endif
701 701
702 // manually update alphar, alphai, betar 702 // manually update alphar, alphai, betar
703 static int jj; 703 static int jj;
704 704
705 jj=0; 705 jj=0;
706 while (jj < nn) 706 while (jj < nn)
707 { 707 {
708 #ifdef DEBUG_EIG 708 #ifdef DEBUG_EIG
709 std::cout << "computing gen eig #" << jj << std::endl; 709 std::cout << "computing gen eig #" << jj << std::endl;
710 #endif 710 #endif
711 711
712 static int zcnt; // number of zeros in this block 712 static int zcnt; // number of zeros in this block
713 713
714 if (jj == (nn-1)) 714 if (jj == (nn-1))
715 zcnt = 1; 715 zcnt = 1;
716 else if (aa(jj+1,jj) == 0) 716 else if (aa(jj+1,jj) == 0)
717 zcnt = 1; 717 zcnt = 1;
718 else zcnt = 2; 718 else zcnt = 2;
719 719
720 if (zcnt == 1) // real zero 720 if (zcnt == 1) // real zero
721 { 721 {
722 #ifdef DEBUG_EIG 722 #ifdef DEBUG_EIG
723 std::cout << " single gen eig:" << std::endl; 723 std::cout << " single gen eig:" << std::endl;
724 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl; 724 std::cout << " alphar(" << jj << ") = " << aa(jj,jj) << std::endl;
725 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl; 725 std::cout << " betar( " << jj << ") = " << bb(jj,jj) << std::endl;
726 std::cout << " alphai(" << jj << ") = 0" << std::endl; 726 std::cout << " alphai(" << jj << ") = 0" << std::endl;
727 #endif 727 #endif
728 728
729 alphar(jj) = aa(jj,jj); 729 alphar(jj) = aa(jj,jj);
730 alphai(jj) = 0; 730 alphai(jj) = 0;
731 betar(jj) = bb(jj,jj); 731 betar(jj) = bb(jj,jj);
732 } 732 }
733 else 733 else
734 { 734 {
735 // complex conjugate pair 735 // complex conjugate pair
736 #ifdef DEBUG_EIG 736 #ifdef DEBUG_EIG
737 std::cout << "qz: calling dlag2:" << std::endl; 737 std::cout << "qz: calling dlag2:" << std::endl;
738 std::cout << "safmin=" 738 std::cout << "safmin="
739 << setiosflags (std::ios::scientific) << safmin << std::endl; 739 << setiosflags (std::ios::scientific) << safmin << std::endl;
740 740
741 for (int idr = jj; idr <= jj+1; idr++) 741 for (int idr = jj; idr <= jj+1; idr++)
742 { 742 {
743 for (int idc = jj; idc <= jj+1; idc++) 743 for (int idc = jj; idc <= jj+1; idc++)
744 { 744 {
745 std::cout << "aa(" << idr << "," << idc << ")=" 745 std::cout << "aa(" << idr << "," << idc << ")="
746 << aa(idr,idc) << std::endl; 746 << aa(idr,idc) << std::endl;
747 std::cout << "bb(" << idr << "," << idc << ")=" 747 std::cout << "bb(" << idr << "," << idc << ")="
748 << bb(idr,idc) << std::endl; 748 << bb(idr,idc) << std::endl;
749 } 749 }
750 } 750 }
751 #endif 751 #endif
752 752
753 // FIXME -- probably should be using 753 // FIXME -- probably should be using
754 // fortran_vec instead of &aa(jj,jj) here. 754 // fortran_vec instead of &aa(jj,jj) here.
755 755
756 double scale1, scale2, wr1, wr2, wi; 756 double scale1, scale2, wr1, wr2, wi;
757 const double *aa_ptr = aa.data () + jj*nn+jj; 757 const double *aa_ptr = aa.data () + jj*nn+jj;
758 const double *bb_ptr = bb.data () + jj*nn+jj; 758 const double *bb_ptr = bb.data () + jj*nn+jj;
759 F77_XFCN (dlag2, DLAG2, 759 F77_XFCN (dlag2, DLAG2,
760 (aa_ptr, nn, bb_ptr, nn, safmin, 760 (aa_ptr, nn, bb_ptr, nn, safmin,
761 scale1, scale2, wr1, wr2, wi)); 761 scale1, scale2, wr1, wr2, wi));
762 762
763 #ifdef DEBUG_EIG 763 #ifdef DEBUG_EIG
764 std::cout << "dlag2 returns: scale1=" << scale1 764 std::cout << "dlag2 returns: scale1=" << scale1
765 << "\tscale2=" << scale2 << std::endl 765 << "\tscale2=" << scale2 << std::endl
766 << "\twr1=" << wr1 << "\twr2=" << wr2 766 << "\twr1=" << wr1 << "\twr2=" << wr2
767 << "\twi=" << wi << std::endl; 767 << "\twi=" << wi << std::endl;
768 #endif 768 #endif
769 769
770 // just to be safe, check if it's a real pair 770 // just to be safe, check if it's a real pair
771 if (wi == 0) 771 if (wi == 0)
772 { 772 {
773 alphar(jj) = wr1; 773 alphar(jj) = wr1;
774 alphai(jj) = 0; 774 alphai(jj) = 0;
775 betar(jj) = scale1; 775 betar(jj) = scale1;
776 alphar(jj+1) = wr2; 776 alphar(jj+1) = wr2;
777 alphai(jj+1) = 0; 777 alphai(jj+1) = 0;
778 betar(jj+1) = scale2; 778 betar(jj+1) = scale2;
779 } 779 }
780 else 780 else
781 { 781 {
782 alphar(jj) = alphar(jj+1)=wr1; 782 alphar(jj) = alphar(jj+1)=wr1;
783 alphai(jj) = -(alphai(jj+1) = wi); 783 alphai(jj) = -(alphai(jj+1) = wi);
784 betar(jj) = betar(jj+1) = scale1; 784 betar(jj) = betar(jj+1) = scale1;
785 } 785 }
786 } 786 }
787 787
788 // advance past this block 788 // advance past this block
789 jj += zcnt; 789 jj += zcnt;
790 } 790 }
791 791
792 #ifdef DEBUG_SORT 792 #ifdef DEBUG_SORT
793 std::cout << "qz: back from dsubsp: aa=" << std::endl; 793 std::cout << "qz: back from dsubsp: aa=" << std::endl;
794 octave_print_internal (std::cout, aa, 0); 794 octave_print_internal (std::cout, aa, 0);
795 std::cout << std::endl << "bb=" << std::endl; 795 std::cout << std::endl << "bb=" << std::endl;
796 octave_print_internal (std::cout, bb, 0); 796 octave_print_internal (std::cout, bb, 0);
797 797
798 if (compz == 'V') 798 if (compz == 'V')
799 { 799 {
800 std::cout << std::endl << "ZZ=" << std::endl; 800 std::cout << std::endl << "ZZ=" << std::endl;
801 octave_print_internal (std::cout, ZZ, 0); 801 octave_print_internal (std::cout, ZZ, 0);
802 } 802 }
803 std::cout << std::endl << "qz: ndim=" << ndim << std::endl 803 std::cout << std::endl << "qz: ndim=" << ndim << std::endl
804 << "fail=" << fail << std::endl; 804 << "fail=" << fail << std::endl;
805 std::cout << "alphar = " << std::endl; 805 std::cout << "alphar = " << std::endl;
806 octave_print_internal (std::cout, (Matrix) alphar, 0); 806 octave_print_internal (std::cout, (Matrix) alphar, 0);
807 std::cout << std::endl << "alphai = " << std::endl; 807 std::cout << std::endl << "alphai = " << std::endl;
808 octave_print_internal (std::cout, (Matrix) alphai, 0); 808 octave_print_internal (std::cout, (Matrix) alphai, 0);
809 std::cout << std::endl << "beta = " << std::endl; 809 std::cout << std::endl << "beta = " << std::endl;
810 octave_print_internal (std::cout, (Matrix) betar, 0); 810 octave_print_internal (std::cout, (Matrix) betar, 0);
811 std::cout << std::endl; 811 std::cout << std::endl;
812 #endif 812 #endif
813 } 813 }
814 } 814 }
815 815
816 // compute generalized eigenvalues? 816 // compute generalized eigenvalues?
817 ComplexColumnVector gev; 817 ComplexColumnVector gev;
818 818
819 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) 819 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
820 { 820 {
821 if (complex_case) 821 if (complex_case)
822 { 822 {
823 error ("complex case not yet implemented"); 823 error ("complex case not yet implemented");
824 return retval; 824 return retval;
825 } 825 }
826 else 826 else
827 { 827 {
828 #ifdef DEBUG 828 #ifdef DEBUG
829 std::cout << "qz: computing generalized eigenvalues" << std::endl; 829 std::cout << "qz: computing generalized eigenvalues" << std::endl;
830 #endif 830 #endif
831 831
832 // return finite generalized eigenvalues 832 // return finite generalized eigenvalues
833 int cnt = 0; 833 int cnt = 0;
834 834
835 for (int ii = 0; ii < nn; ii++) 835 for (int ii = 0; ii < nn; ii++)
836 if (betar(ii) != 0) 836 if (betar(ii) != 0)
837 cnt++; 837 cnt++;
838 838
839 ComplexColumnVector tmp(cnt); 839 ComplexColumnVector tmp(cnt);
840 840
841 cnt = 0; 841 cnt = 0;
842 for (int ii = 0; ii < nn; ii++) 842 for (int ii = 0; ii < nn; ii++)
843 if (betar(ii) != 0) 843 if (betar(ii) != 0)
844 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii); 844 tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii);
845 gev = tmp; 845 gev = tmp;
846 } 846 }
847 } 847 }
848 848
849 // right, left eigenvector matrices 849 // right, left eigenvector matrices
850 if (nargout >= 5) 850 if (nargout >= 5)
851 { 851 {
852 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute? 852 char side = (nargout == 5 ? 'R' : 'B'); // which side to compute?
853 char howmny = 'B'; // compute all of them and backtransform 853 char howmny = 'B'; // compute all of them and backtransform
854 octave_idx_type *select = 0; // dummy pointer; select is not used. 854 octave_idx_type *select = 0; // dummy pointer; select is not used.
855 octave_idx_type m; 855 octave_idx_type m;
856 856
857 if (complex_case) 857 if (complex_case)
858 { 858 {
859 error ("complex type not yet implemented"); 859 error ("complex type not yet implemented");
860 return retval; 860 return retval;
861 } 861 }
862 else 862 else
863 { 863 {
864 #ifdef DEBUG 864 #ifdef DEBUG
865 std::cout << "qz: computing generalized eigenvectors" << std::endl; 865 std::cout << "qz: computing generalized eigenvectors" << std::endl;
866 #endif 866 #endif
867 867
868 VL = QQ; 868 VL = QQ;
869 VR = ZZ; 869 VR = ZZ;
870 870
871 F77_XFCN (dtgevc, DTGEVC, 871 F77_XFCN (dtgevc, DTGEVC,
872 (F77_CONST_CHAR_ARG2 (&side, 1), 872 (F77_CONST_CHAR_ARG2 (&side, 1),
873 F77_CONST_CHAR_ARG2 (&howmny, 1), 873 F77_CONST_CHAR_ARG2 (&howmny, 1),
874 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (), 874 select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
875 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn, 875 nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
876 m, work.fortran_vec (), info 876 m, work.fortran_vec (), info
877 F77_CHAR_ARG_LEN (1) 877 F77_CHAR_ARG_LEN (1)
878 F77_CHAR_ARG_LEN (1))); 878 F77_CHAR_ARG_LEN (1)));
879 879
880 // now construct the complex form of VV, WW 880 // now construct the complex form of VV, WW
881 int jj = 0; 881 int jj = 0;
882 882
883 while (jj < nn) 883 while (jj < nn)
884 { 884 {
885 OCTAVE_QUIT; 885 OCTAVE_QUIT;
886 886
887 // see if real or complex eigenvalue 887 // see if real or complex eigenvalue
888 int cinc = 2; // column increment; assume complex eigenvalue 888 int cinc = 2; // column increment; assume complex eigenvalue
889 889
890 if (jj == (nn-1)) 890 if (jj == (nn-1))
891 cinc = 1; // single column 891 cinc = 1; // single column
892 else if (aa(jj+1,jj) == 0) 892 else if (aa(jj+1,jj) == 0)
893 cinc = 1; 893 cinc = 1;
894 894
895 // now copy the eigenvector (s) to CVR, CVL 895 // now copy the eigenvector (s) to CVR, CVL
896 if (cinc == 1) 896 if (cinc == 1)
897 { 897 {
898 for (int ii = 0; ii < nn; ii++) 898 for (int ii = 0; ii < nn; ii++)
899 CVR(ii,jj) = VR(ii,jj); 899 CVR(ii,jj) = VR(ii,jj);
900 900
901 if (side == 'B') 901 if (side == 'B')
902 for (int ii = 0; ii < nn; ii++) 902 for (int ii = 0; ii < nn; ii++)
903 CVL(ii,jj) = VL(ii,jj); 903 CVL(ii,jj) = VL(ii,jj);
904 } 904 }
905 else 905 else
906 { 906 {
907 // double column; complex vector 907 // double column; complex vector
908 908
909 for (int ii = 0; ii < nn; ii++) 909 for (int ii = 0; ii < nn; ii++)
910 { 910 {
911 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); 911 CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
912 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); 912 CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
913 } 913 }
914 914
915 if (side == 'B') 915 if (side == 'B')
916 for (int ii = 0; ii < nn; ii++) 916 for (int ii = 0; ii < nn; ii++)
917 { 917 {
918 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); 918 CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
919 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); 919 CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
920 } 920 }
921 } 921 }
922 922
923 // advance to next eigenvectors (if any) 923 // advance to next eigenvectors (if any)
924 jj += cinc; 924 jj += cinc;
925 } 925 }
926 } 926 }
927 } 927 }
928 928
929 switch (nargout) 929 switch (nargout)
930 { 930 {
931 case 7: 931 case 7:
932 retval(6) = gev; 932 retval(6) = gev;
933 933
934 case 6: // return eigenvectors 934 case 6: // return eigenvectors
935 retval(5) = CVL; 935 retval(5) = CVL;
936 936
937 case 5: // return eigenvectors 937 case 5: // return eigenvectors
938 retval(4) = CVR; 938 retval(4) = CVR;
939 939
940 case 4: 940 case 4:
941 if (nargin == 3) 941 if (nargin == 3)
942 { 942 {
943 #ifdef DEBUG 943 #ifdef DEBUG
944 std::cout << "qz: sort: retval(3) = gev = " << std::endl; 944 std::cout << "qz: sort: retval(3) = gev = " << std::endl;
945 octave_print_internal (std::cout, gev); 945 octave_print_internal (std::cout, gev);
946 std::cout << std::endl; 946 std::cout << std::endl;
947 #endif 947 #endif
948 retval(3) = gev; 948 retval(3) = gev;
949 } 949 }
950 else 950 else
951 retval(3) = ZZ; 951 retval(3) = ZZ;
952 952
953 case 3: 953 case 3:
954 if (nargin == 3) 954 if (nargin == 3)
955 retval(2) = ZZ; 955 retval(2) = ZZ;
956 else 956 else
957 retval(2) = QQ; 957 retval(2) = QQ;
958 958
959 case 2: 959 case 2:
960 #ifdef DEBUG 960 #ifdef DEBUG
961 std::cout << "qz: retval (1) = bb = " << std::endl; 961 std::cout << "qz: retval (1) = bb = " << std::endl;
962 octave_print_internal (std::cout, bb, 0); 962 octave_print_internal (std::cout, bb, 0);