Mercurial > octave-nkf

changeset 10307:4e4270ab70d6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
support complex case in qz
author Jyh-miin Lin <jyhmiinlin@gmail.com>
date Thu, 11 Feb 2010 08:13:49 +0100
parents 7b5e8527441e
children 0d928dd9eeb8
files src/ChangeLog src/DLD-FUNCTIONS/qz.cc
diffstat 2 files changed, 232 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog	Thu Feb 11 07:52:00 2010 +0100
+++ b/src/ChangeLog	Thu Feb 11 08:13:49 2010 +0100
@@ -1,3 +1,9 @@
+2010-01-30  Jyh-Miin Lin <jyhmiin@gmail.com>
+	* DLD-FUNCTIONS/qz.cc (Fqz):
+	  Add support for complex case (no reordering yet).
+	  Place dggbak and zggbak later.
+	  Change Q to Q' to meet MATLAB's results.
+
 2010-02-11  Jaroslav Hajek  <highegg@gmail.com>
 
 	* symtab.cc: Reverse the effect of 2ceae0b40515.
--- a/src/DLD-FUNCTIONS/qz.cc	Thu Feb 11 07:52:00 2010 +0100
+++ b/src/DLD-FUNCTIONS/qz.cc	Thu Feb 11 08:13:49 2010 +0100
@@ -39,6 +39,7 @@
 #include <iomanip>
 
 #include "CmplxQRP.h"
+#include "CmplxQR.h"
 #include "dbleQR.h"
 #include "f77-fcn.h"
 #include "lo-math.h"
@@ -72,6 +73,14 @@
                              double* WORK, octave_idx_type& INFO
                              F77_CHAR_ARG_LEN_DECL);
 
+F77_RET_T
+  F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type& N, Complex* A, const octave_idx_type& LDA,
+                             Complex* B, const octave_idx_type& LDB, octave_idx_type& ILO,
+                             octave_idx_type& IHI, double* LSCALE, double* RSCALE,
+                             double* WORK, octave_idx_type& INFO
+                             F77_CHAR_ARG_LEN_DECL);	
+
   F77_RET_T
   F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
                              F77_CONST_CHAR_ARG_DECL,
@@ -82,6 +91,16 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
 
+F77_RET_T
+  F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type& N, const octave_idx_type& ILO,
+                             const octave_idx_type& IHI, const double* LSCALE,
+                             const double* RSCALE, octave_idx_type& M, Complex* V,
+                             const octave_idx_type& LDV, octave_idx_type& INFO
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   F77_RET_T
   F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL,
                              F77_CONST_CHAR_ARG_DECL,
@@ -94,6 +113,18 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
 
+ F77_RET_T
+  F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type& N, const octave_idx_type& ILO,
+                             const octave_idx_type& IHI, Complex* A,
+                             const octave_idx_type& LDA, Complex* B,
+                             const octave_idx_type& LDB, Complex* Q,
+                             const octave_idx_type& LDQ, Complex* Z,
+                             const octave_idx_type& LDZ, octave_idx_type& INFO
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   F77_RET_T
   F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL,
                              F77_CONST_CHAR_ARG_DECL,
@@ -109,6 +140,21 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
 
+F77_RET_T
+  F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type& N, const octave_idx_type& ILO, const octave_idx_type& IHI,
+                             Complex* A, const octave_idx_type& LDA, Complex* B,
+                             const octave_idx_type& LDB, Complex* ALPHA, Complex* BETA, Complex* CQ, const octave_idx_type& LDQ, 
+			     Complex* CZ, const octave_idx_type& LDZ, 
+			     Complex* WORK,
+                             const octave_idx_type& LWORK, double* RWORK,
+			     octave_idx_type& INFO
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   F77_RET_T
   F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA, const double* B,
                            const octave_idx_type& LDB, const double& SAFMIN,
@@ -138,6 +184,18 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
 
+F77_RET_T
+  F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             octave_idx_type* SELECT, const octave_idx_type& N,const Complex* A,
+                             const octave_idx_type& LDA,const Complex* B,
+                             const octave_idx_type& LDB, Complex* xVL,
+                             const octave_idx_type& LDVL, Complex* xVR,
+                             const octave_idx_type& LDVR, const octave_idx_type& MM,
+                             octave_idx_type& M, Complex* CWORK, double* RWORK, octave_idx_type& INFO
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   F77_RET_T
   F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL,
                                double& retval
@@ -237,7 +295,7 @@
 \n\
     A*V = B*V*diag(lambda)\n\
     W'*A = diag(lambda)*W'*B\n\
-    AA = Q'*A*Z, BB = Q'*B*Z\n\
+    AA = Q*A*Z, BB = Q*B*Z\n\
 \n\
 @end group\n\
 @end example\n\
@@ -426,9 +484,9 @@
   // Both matrices loaded, now let's check what kind of arithmetic:
   //declared static to avoid compiler warnings about long jumps, vforks.
 
-  static int complex_case
+  int complex_case
     = (args(0).is_complex_type () || args(1).is_complex_type ());
-
+  // static complex_case causing random segfault, so it is removed
   if (nargin == 3 && complex_case)
     {
       error ("qz: cannot re-order complex qz decomposition.");
@@ -438,7 +496,7 @@
   // first, declare variables used in both the real and complex case
   Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
   RowVector alphar(nn), alphai(nn), betar(nn);
-
+  ComplexRowVector xalpha(nn), xbeta(nn);
   ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
   octave_idx_type ilo, ihi, info;
   char compq = (nargout >= 3 ? 'V' : 'N');
@@ -455,12 +513,32 @@
 
   // always perform permutation balancing
   const char bal_job = 'P';
-  RowVector lscale(nn), rscale(nn), work(6*nn);
+  RowVector lscale(nn), rscale(nn), work(6*nn), rwork(nn);
 
   if (complex_case)
     {
-      error ("Complex case not implemented yet");
-      return retval;
+#ifdef DEBUG
+      if (compq == 'V')
+	std::cout << "qz: performing balancing; CQ=" << std::endl << CQ << std::endl;
+#endif
+      if (args(0).is_real_type ())
+	caa = ComplexMatrix (aa);
+
+      if (args(1).is_real_type ())
+	cbb = ComplexMatrix (bb);
+
+      if (compq == 'V')
+	CQ = ComplexMatrix (QQ);
+
+      if (compz == 'V')
+	CZ = ComplexMatrix (ZZ);
+
+          F77_XFCN (zggbal, ZGGBAL,
+		(F77_CONST_CHAR_ARG2 (&bal_job, 1),
+		 nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
+		 nn, ilo, ihi, lscale.fortran_vec (),
+		 rscale.fortran_vec (), work.fortran_vec (), info
+		 F77_CHAR_ARG_LEN (1)));
     }
   else
     {
@@ -481,7 +559,7 @@
   // for both the real and complex cases;
   // left first
 
-  if (compq == 'V')
+ /* if (compq == 'V')
     {
       F77_XFCN (dggbak, DGGBAK,
                 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
@@ -512,7 +590,7 @@
       if (compz == 'V')
         std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
 #endif
-    }
+    }   */
 
   static char qz_job;
   qz_job = (nargout < 2 ? 'E' : 'S');   
@@ -520,20 +598,59 @@
   if (complex_case)
     {
       // complex case
-      if (args(0).is_real_type ())
-        caa = ComplexMatrix (aa);
-
-      if (args(1).is_real_type ())
-        cbb = ComplexMatrix (bb);
+      ComplexQR cbqr (cbb); // declare cbqr as the QR decomposition of cbb
+      cbb = cbqr.R ();  // the R matrix of QR decomposition for cbb
+      caa = (cbqr.Q ().hermitian ())*caa; // (Q*)caa for following work
+      //if (compq == 'V')
+      CQ = CQ*cbqr.Q ();
+      F77_XFCN (zgghrd, ZGGHRD,
+                (F77_CONST_CHAR_ARG2 (&compq, 1),
+                 F77_CONST_CHAR_ARG2 (&compz, 1),
+                 nn, ilo, ihi, caa.fortran_vec (),
+                 nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn,
+                 CZ.fortran_vec (), nn, info
+                 F77_CHAR_ARG_LEN (1)
+                 F77_CHAR_ARG_LEN (1)));
+      ComplexRowVector cwork(1*nn); 
+      F77_XFCN (zhgeqz, ZHGEQZ,
+                (F77_CONST_CHAR_ARG2 (&qz_job, 1),
+                 F77_CONST_CHAR_ARG2 (&compq, 1),
+                 F77_CONST_CHAR_ARG2 (&compz, 1),
+                 nn, ilo, ihi, 
+                 caa.fortran_vec (), nn, 
+                 cbb.fortran_vec (),nn, 
+                 xalpha.fortran_vec (), xbeta.fortran_vec (), 
+                 CQ.fortran_vec (), nn,
+                 CZ.fortran_vec (), nn, 
+                 cwork.fortran_vec (), nn, rwork.fortran_vec (), info
+                 F77_CHAR_ARG_LEN (1)
+                 F77_CHAR_ARG_LEN (1)
+                 F77_CHAR_ARG_LEN (1))); 
 
-      if (compq == 'V')
-        CQ = ComplexMatrix (QQ);
+      if (compq == 'V') 
+        {
+          // Left eigenvector
+          F77_XFCN (zggbak, ZGGBAK,
+                    (F77_CONST_CHAR_ARG2 (&bal_job, 1),
+                     F77_CONST_CHAR_ARG2 ("L", 1),
+                     nn, ilo, ihi, lscale.data (), rscale.data (),
+                     nn, CQ.fortran_vec (), nn, info
+                     F77_CHAR_ARG_LEN (1)
+                     F77_CHAR_ARG_LEN (1)));
+        }        
 
+      // then right
       if (compz == 'V')
-        CZ = ComplexMatrix (ZZ);
+        {
+          F77_XFCN (zggbak, ZGGBAK,
+                    (F77_CONST_CHAR_ARG2 (&bal_job, 1),
+                     F77_CONST_CHAR_ARG2 ("R", 1),
+                     nn, ilo, ihi, lscale.data (), rscale.data (),
+                     nn, CZ.fortran_vec (), nn, info
+                     F77_CHAR_ARG_LEN (1)
+                     F77_CHAR_ARG_LEN (1)));
+        } 
 
-      error ("complex case not done yet");
-      return retval;
     }
   else          // real matrices case
     {
@@ -600,6 +717,39 @@
                  F77_CHAR_ARG_LEN (1)
                  F77_CHAR_ARG_LEN (1)
                  F77_CHAR_ARG_LEN (1)));
+      if (compq == 'V')
+       {
+         F77_XFCN (dggbak, DGGBAK,
+		(F77_CONST_CHAR_ARG2 (&bal_job, 1),
+		 F77_CONST_CHAR_ARG2 ("L", 1),
+		 nn, ilo, ihi, lscale.data (), rscale.data (),
+		 nn, QQ.fortran_vec (), nn, info
+		 F77_CHAR_ARG_LEN (1)
+		 F77_CHAR_ARG_LEN (1)));
+
+#ifdef DEBUG
+         if (compq == 'V')
+           std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
+#endif
+        }
+
+  // then right
+      if (compz == 'V')
+        {
+           F77_XFCN (dggbak, DGGBAK,
+		(F77_CONST_CHAR_ARG2 (&bal_job, 1),
+		 F77_CONST_CHAR_ARG2 ("R", 1),
+		 nn, ilo, ihi, lscale.data (), rscale.data (),
+		 nn, ZZ.fortran_vec (), nn, info
+		 F77_CHAR_ARG_LEN (1)
+		 F77_CHAR_ARG_LEN (1)));
+
+#ifdef DEBUG
+           if (compz == 'V')
+             std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
+#endif
+        }
+
     }
 
   // order the QZ decomposition?
@@ -820,8 +970,19 @@
     {
       if (complex_case)
         {
-          error ("complex case not yet implemented");
-          return retval;
+	  int cnt = 0;
+
+	  for (int ii = 0; ii < nn; ii++)
+	//    if (cbetar(ii) != 0)
+	      cnt++;
+
+	  ComplexColumnVector tmp(cnt);
+
+	  cnt = 0;
+	  for (int ii = 0; ii < nn; ii++)
+	  //  if (cbetar(ii) != 0)
+	      tmp(cnt++) = xalpha(ii)/xbeta(ii);
+	  gev = tmp;
         }
       else
         {
@@ -852,12 +1013,24 @@
       char side = (nargout == 5 ? 'R' : 'B');   // which side to compute?
       char howmny = 'B';  // compute all of them and backtransform
       octave_idx_type *select = 0; // dummy pointer; select is not used.
-      octave_idx_type m;
 
       if (complex_case)
         {
-          error ("complex type not yet implemented");
-          return retval;
+          octave_idx_type m;
+	  CVL=CQ;
+	  CVR=CZ;
+	  ComplexRowVector cwork2(2*nn);
+	  RowVector rwork2(8*nn);
+
+	  //octave_idx_type n=nn;
+	  F77_XFCN (ztgevc, ZTGEVC,
+		    (F77_CONST_CHAR_ARG2 (&side, 1),
+		     F77_CONST_CHAR_ARG2 (&howmny, 1),
+		     select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
+		     nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn,
+		     m, cwork2.fortran_vec (), rwork2.fortran_vec (), info
+		     F77_CHAR_ARG_LEN (1)
+		     F77_CHAR_ARG_LEN (1)));	 
         }
       else
         {
@@ -867,7 +1040,7 @@
 
           VL = QQ;
           VR = ZZ;
-
+	  octave_idx_type m;
           F77_XFCN (dtgevc, DTGEVC,
                     (F77_CONST_CHAR_ARG2 (&side, 1),
                      F77_CONST_CHAR_ARG2 (&howmny, 1),
@@ -948,15 +1121,36 @@
           retval(3) = gev;
         }
       else
-        retval(3) = ZZ;
+      {if (complex_case)
+	  retval(3) = CZ;
+	else
+	  retval(3) = ZZ;
+      }
 
     case 3:
       if (nargin == 3)
-        retval(2) = ZZ;
+	retval(2) = CZ;
       else
-        retval(2) = QQ;
-
+      {if (complex_case)
+	retval(2) = CQ.hermitian(); // compabible with MATLAB output
+       else
+	 retval(2) = QQ.transpose();
+      }
     case 2:
+      {if (complex_case)
+      {
+#ifdef DEBUG
+      std::cout << "qz: retval (1) = cbb = " << std::endl;
+      octave_print_internal (std::cout, cbb, 0);
+      std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
+      octave_print_internal (std::cout, caa, 0);
+      std::cout << std::endl;
+#endif	
+      retval(1) = cbb;
+      retval(0) = caa;
+      }
+      else
+      { // real case
 #ifdef DEBUG
       std::cout << "qz: retval (1) = bb = " << std::endl;
       octave_print_internal (std::cout, bb, 0);
@@ -966,7 +1160,9 @@
 #endif
       retval(1) = bb;
       retval(0) = aa;
-      break;
+      }
+      break;}      
+
 
     case 1:
     case 0: