diff src/DLD-FUNCTIONS/balance.cc @ 3185:9580887dd160

[project @ 1998-09-26 02:45:55 by jwe]
author jwe
date Sat, 26 Sep 1998 02:45:59 +0000
parents 3f6a813eb09e
children eba59b8c64dc
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/balance.cc	Fri Sep 25 04:35:09 1998 +0000
+++ b/src/DLD-FUNCTIONS/balance.cc	Sat Sep 26 02:45:59 1998 +0000
@@ -42,24 +42,27 @@
 
 extern "C"
 {
-  int F77_FCN( dggbal, DGGBAL) (const char* JOB,  const int& N,
-        double* A, const int& LDA, double* B, const int& LDB,
-        int& ILO, int & IHI, double* LSCALE,
-        double* RSCALE, double* WORK, int& INFO, long );
+  int F77_FCN (dggbal, DGGBAL) (const char* JOB, const int& N,
+				double* A, const int& LDA, double* B,
+				const int& LDB, int& ILO, int& IHI,
+				double* LSCALE, double* RSCALE,
+				double* WORK, int& INFO, long);
 
-  int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE,
-        const int& N, const int& ILO, const int& IHI,
-        double* LSCALE, double* RSCALE, int& M,
-        double* V, const int& LDV, int& INFO, long, long);
+  int F77_FCN (dggbak, DGGBAK) (const char* JOB, const char* SIDE,
+				const int& N, const int& ILO,
+				const int& IHI, double* LSCALE,
+				double* RSCALE, int& M,	double* V,
+				const int& LDV, int& INFO, long, long);
 
-  int F77_FCN( zggbal, ZGGBAL) ( const char* JOB,  const int& N,
-        Complex* A, const int& LDA, Complex* B, const int& LDB,
-        int& ILO, int & IHI, double* LSCALE,
-        double* RSCALE, double* WORK, int& INFO, long );
+  int F77_FCN (zggbal, ZGGBAL) (const char* JOB, const int& N,
+				Complex* A, const int& LDA, Complex* B,
+				const int& LDB, int& ILO, int& IHI,
+				double* LSCALE, double* RSCALE,
+				double* WORK, int& INFO, long);
 }
 
 DEFUN_DLD (balance, args, nargout,
-  "AA = balance (A [, OPT]) or [[DD,] AA] =  balance (A [, OPT])\n\
+  "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\
 \n\
 generalized eigenvalue problem:\n\
 \n\
@@ -77,7 +80,6 @@
 \n\
 [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)")
 {
-
   octave_value_list retval;
 
   int nargin = args.length ();
@@ -89,7 +91,7 @@
     }
 
   // determine if it's AEP or GEP
-  int AEPcase = (nargin == 1 ? 1 : args(1).is_string() );
+  int AEPcase = nargin == 1 ? 1 : args(1).is_string ();
   string bal_job;
 
   // problem dimension
@@ -99,181 +101,224 @@
 
   if (arg_is_empty < 0)
     return retval;
+
   if (arg_is_empty > 0)
     return octave_value_list (2, Matrix ());
 
   if (nn != args(0).columns())
-  {
-    gripe_square_matrix_required ("balance");
-    return retval;
-  }
+    {
+      gripe_square_matrix_required ("balance");
+      return retval;
+    }
 
   // Extract argument 1 parameter for both AEP and GEP.
   Matrix aa;
   ComplexMatrix caa;
-  if (args(0).is_complex_type ()) caa = args(0).complex_matrix_value ();
-  else                            aa = args(0).matrix_value ();
-  if (error_state)                return retval;
+
+  if (args(0).is_complex_type ())
+    caa = args(0).complex_matrix_value ();
+  else
+    aa = args(0).matrix_value ();
+
+  if (error_state)
+    return retval;
 
   // Treat AEP/GEP cases.
-  if(AEPcase)
-  {   
-    // Algebraic eigenvalue problem.
-    if(nargin == 1)
-      bal_job = "B";
-    else if(args(1).is_string())
-      bal_job = args(1).string_value();
-    // the next line should never execute, but better safe than sorry.
-    else error("balance: AEP argument 2 must be a string");
+  if (AEPcase)
+    {  
+      // Algebraic eigenvalue problem.
 
-    // balance the AEP
-    if (args(0).is_complex_type ())
-    {
-      ComplexAEPBALANCE result (caa, bal_job);
-
-      if (nargout == 0 || nargout == 1)
-        retval(0) = result.balanced_matrix ();
+      if (nargin == 1)
+	bal_job = "B";
+      else if (args(1).is_string ())
+	bal_job = args(1).string_value ();
       else
-      {
-        retval(1) = result.balanced_matrix ();
-        retval(0) = result.balancing_matrix ();
-      }
-    }
-    else
-    {
-      AEPBALANCE result (aa, bal_job);
+	{
+	  error ("balance: AEP argument 2 must be a string");
+	  return retval;
+	}
+
+      // balance the AEP
+      if (args(0).is_complex_type ())
+	{
+	  ComplexAEPBALANCE result (caa, bal_job);
 
-      if (nargout == 0 || nargout == 1)
-        retval(0) = result.balanced_matrix ();
+	  if (nargout == 0 || nargout == 1)
+	    retval(0) = result.balanced_matrix ();
+	  else
+	    {
+	      retval(1) = result.balanced_matrix ();
+	      retval(0) = result.balancing_matrix ();
+	    }
+	}
       else
-      {
-        retval(1) = result.balanced_matrix ();
-        retval(0) = result.balancing_matrix ();
-       }
+	{
+	  AEPBALANCE result (aa, bal_job);
+
+	  if (nargout == 0 || nargout == 1)
+	    retval(0) = result.balanced_matrix ();
+	  else
+	    {
+	      retval(1) = result.balanced_matrix ();
+	      retval(0) = result.balancing_matrix ();
+	    }
+	}
     }
-  }
-  //
-  // end of AEP case, now do GEP case
   else
-  {
-    // Generalized eigenvalue problem.
-    if(nargin == 2)
-      bal_job = "B";
-    else if(args(2).is_string())
-      bal_job = args(2).string_value();
-    else error("balance: GEP argument 3 must be a string");
-
-    if( (nn != args(1).columns()) || (nn != args(1).rows() ) )
-    {
-      gripe_nonconformant ();
-      return retval;
-    }
-    Matrix bb;
-    ComplexMatrix cbb;
-    if (args(1).is_complex_type ()) cbb = args(1).complex_matrix_value ();
-    else                            bb = args(1).matrix_value ();
-    if (error_state) return retval;
-
-    //
-    // Both matrices loaded, now let's check what kind of arithmetic:
-    // first, declare variables used in both the real and complex case
-    int ilo, ihi, info;
-    RowVector lscale(nn), rscale(nn), work(6*nn);
-    char job = bal_job[0];
-    static int complex_case = (args(0).is_complex_type() 
-                       || args(1).is_complex_type());
-
-    // now balance
-    if (complex_case)
-    {
-      if (args(0).is_real_type ()) caa = aa;
-      if (args(1).is_real_type ()) cbb = bb;
-  
-      F77_XFCN( zggbal, ZGGBAL, ( &job, nn, caa.fortran_vec(), nn,
-          cbb.fortran_vec(), nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), work.fortran_vec(), info, 1L));
-    }
-    else          // real matrices case
     {
-      F77_XFCN( dggbal, DGGBAL, (&job,  nn, aa.fortran_vec(),
-          nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), work.fortran_vec(), info , 1L));
-      
-      if(f77_exception_encountered)
-        (*current_liboctave_error_handler) 
-         ("unrecoverable error in balance GEP");
-    }
-      
-    // Since we just want the balancing matrices, we can use dggbal
-    // for both the real and complex cases;
-    Matrix Pl(nn,nn), Pr(nn,nn);
-    for(int ii=0; ii < nn ; ii++)
-      for( int jj=0; jj < nn ; jj++)
-        Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0);
+      // Generalized eigenvalue problem.
+      if (nargin == 2)
+	bal_job = "B";
+      else if (args(2).is_string ())
+	bal_job = args(2).string_value ();
+      else
+	{
+	  error ("balance: GEP argument 3 must be a string");
+	  return retval;
+	}
+
+      if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
+	{
+	  gripe_nonconformant ();
+	  return retval;
+	}
+
+      Matrix bb;
+      ComplexMatrix cbb;
+
+      if (args(1).is_complex_type ())
+	cbb = args(1).complex_matrix_value ();
+      else
+	bb = args(1).matrix_value ();
+
+      if (error_state)
+	return retval;
+
+      // Both matrices loaded, now let's check what kind of arithmetic:
+      // first, declare variables used in both the real and complex case
+
+      int ilo, ihi, info;
+      RowVector lscale(nn), rscale(nn), work(6*nn);
+      char job = bal_job[0];
+
+      static int complex_case
+	= (args(0).is_complex_type () || args(1).is_complex_type ());
+
+      // now balance
+      if (complex_case)
+	{
+	  if (args(0).is_real_type ())
+	    caa = aa;
+
+	  if (args(1).is_real_type ())
+	    cbb = bb;
   
-    // left first
-    F77_XFCN( dggbak, DGGBAK, (&job, "L",
-          nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), nn, Pl.fortran_vec(),
-          nn, info, 1L, 1L));
-      
-    if(f77_exception_encountered)
-      (*current_liboctave_error_handler) 
-        ("unrecoverable error in balance GEP(L)");
-      
-    // then right
-    F77_XFCN(dggbak, DGGBAK, (&job, "R",
-          nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), nn, Pr.fortran_vec(),
-          nn, info, 1L, 1L));
-    if(f77_exception_encountered)
-      (*current_liboctave_error_handler) 
-        ("unrecoverable error in balance GEP(R)");
+	  F77_XFCN (zggbal, ZGGBAL,
+		    (&job, nn, caa.fortran_vec(), nn,
+		     cbb.fortran_vec(), nn, ilo, ihi,
+		     lscale.fortran_vec(), rscale.fortran_vec(),
+		     work.fortran_vec(), info, 1L));
+
+	  if (f77_exception_encountered)
+	    {
+	      error ("unrecoverable error in balance GEP");
+	      return retval;
+	    }
+	}
+      else
+	{
+	  // real matrices case
 
-    switch (nargout)
-    {
-    case 0:
-    case 1:
-      warning ("balance: used GEP, should have two output arguments");
-      if(complex_case)
-        retval(0) = caa;
-      else
-        retval(0) = aa;
-      break;
+	  F77_XFCN (dggbal, DGGBAL,
+		    (&job,  nn, aa.fortran_vec(), nn, bb.fortran_vec(),
+		     nn, ilo, ihi, lscale.fortran_vec(),
+		     rscale.fortran_vec(), work.fortran_vec(), info, 1L));
+      
+	  if (f77_exception_encountered)
+	    {
+	      error ("unrecoverable error in balance GEP");
+	      return retval;
+	    }
+	}
+      
+      // Since we just want the balancing matrices, we can use dggbal
+      // for both the real and complex cases.
+
+      Matrix Pl(nn,nn), Pr(nn,nn);
+
+      for (int ii = 0; ii < nn; ii++)
+	for (int jj = 0; jj < nn; jj++)
+	  Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0);
+  
+      // left first
+      F77_XFCN (dggbak, DGGBAK,
+		(&job, "L", nn, ilo, ihi, lscale.fortran_vec(),
+		 rscale.fortran_vec(), nn, Pl.fortran_vec(),
+		 nn, info, 1L, 1L));
+      
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in balance GEP(L)");
+	  return retval;
+	}
+      
+      // then right
+      F77_XFCN (dggbak, DGGBAK,
+		(&job, "R", nn, ilo, ihi, lscale.fortran_vec(),
+		 rscale.fortran_vec(), nn, Pr.fortran_vec(),
+		 nn, info, 1L, 1L));
+
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in balance GEP(R)");
+	  return retval;
+	}
 
-    case 2:
-      if(complex_case)
-      {
-        retval(1) = cbb;
-        retval(0) = caa;
-      }
-      else
-      {
-        retval(1) = bb;
-        retval(0) = aa;
-      }
-      break;
+      switch (nargout)
+	{
+	case 0:
+	case 1:
+	  warning ("balance: used GEP, should have two output arguments");
+	  if (complex_case)
+	    retval(0) = caa;
+	  else
+	    retval(0) = aa;
+	  break;
 
-    case 4:
-      if(complex_case)
-      {
-        retval(3) = cbb;
-        retval(2) = caa;
-      }
-      else
-      {
-        retval(3) = bb;
-        retval(2) = aa;
-      }
-      retval(1) = Pr;
-      retval(0) = Pl.transpose();  // so that aa_bal = cc*aa*dd, etc.
-      break;
+	case 2:
+	  if (complex_case)
+	    {
+	      retval(1) = cbb;
+	      retval(0) = caa;
+	    }
+	  else
+	    {
+	      retval(1) = bb;
+	      retval(0) = aa;
+	    }
+	  break;
 
-    default:
-      error ("balance: invalid number of output arguments");
-      break;
+	case 4:
+	  if (complex_case)
+	    {
+	      retval(3) = cbb;
+	      retval(2) = caa;
+	    }
+	  else
+	    {
+	      retval(3) = bb;
+	      retval(2) = aa;
+	    }
+	  retval(1) = Pr;
+	  retval(0) = Pl.transpose ();  // so that aa_bal = cc*aa*dd, etc.
+	  break;
+
+	default:
+	  error ("balance: invalid number of output arguments");
+	  break;
+	}
     }
-  }
+
   return retval;
 }