# HG changeset patch
# User jwe
# Date 745022000 0
# Node ID f3215b07c171d45a223c4e448dee3dfd5a2a5577
# Parent  e399beacf758d6ac5cdfddb315eafc79de273436
[project @ 1993-08-10 22:33:20 by jwe]
Initial revision

diff -r e399beacf758 -r f3215b07c171 src/qzval.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/qzval.cc	Tue Aug 10 22:33:20 1993 +0000
@@ -0,0 +1,165 @@
+// tc-qzval.cc                                           -*- C++ -*-
+/*
+
+Copyright (C) 1993 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+// Written by A. S. Hodel <scotte@eng.auburn.edu>
+
+#ifdef __GNUG__
+#pragma implementation
+#endif
+
+#include "Matrix.h"
+
+#include "tree-const.h"
+#include "user-prefs.h"
+#include "gripes.h"
+#include "error.h"
+
+int F77_FCN (qzhes) (const int*, const int*, double*, double*, const
+		     long*, double*);
+ 
+int F77_FCN (qzit) (const int*, const int*, double*, double*, const
+		    double*, const long*, double*, int*);
+ 
+int F77_FCN (qzval) (const int*, const int*, double*, double*,
+		     double*, double*, double*, const long*, double*);
+
+// XXX FIXME XXX
+extern int empty_arg (tree_constant&);
+extern tree_constant* empty_tree (int , char*);
+extern ComplexMatrix ComplexMatrixLoad (tree_constant&);
+
+#ifdef WITH_DLD
+tree_constant *
+builtin_qzvalue_2 (tree_constant *args, int nargin, int nargout)
+{
+  return qzvalue (args, nargin, nargout);
+}
+#endif
+
+tree_constant *
+qzvalue (tree_constant *args, int nargin, int nargout)
+{
+  tree_constant *retval = NULL_TREE_CONST;
+
+  tree_constant arga = args[1].make_numeric ();
+  tree_constant argb = args[2].make_numeric();
+
+  if (empty_arg (arga) || empty_arg (argb))
+    retval = empty_tree (nargout, "qzvalue");
+  else
+    {
+
+// Arguments are not empty, so check for correct dimensions.
+
+      int a_rows = arga.rows();
+      int a_cols = arga.columns();
+      int b_rows = argb.rows();
+      int b_cols = argb.columns();
+  
+      if ((a_rows != a_cols) || (b_rows != b_cols))
+	{
+	  gripe_square_matrix_required ("qzvalue: first two parameters:");
+	  return retval;
+	}
+      else if (a_rows != b_rows)
+	{
+	  gripe_nonconformant ();
+	  return retval;
+	}
+  
+// Dimensions look o.k., let's solve the problem.
+
+      retval = new tree_constant[nargout+1];
+
+      if (arga.is_complex_type () || argb.is_complex_type ())
+	error ("qzvalue: cannot yet do complex matrix arguments\n");
+      else
+	{
+
+// Do everything in real arithmetic.
+
+	  Matrix jnk (a_rows, a_rows, 0.0);
+
+	  ColumnVector alfr (a_rows);
+	  ColumnVector alfi (a_rows);
+	  ColumnVector beta (a_rows);
+
+	  long matz = 0;
+	  int info;
+
+// XXX FIXME ??? XXX
+	  double eps = DBL_EPSILON;
+
+	  Matrix ca = arga.to_matrix ();
+	  Matrix cb = argb.to_matrix ();
+
+// Use EISPACK qz functions.
+
+	  F77_FCN (qzhes) (&a_rows, &a_rows, ca.fortran_vec (),
+			   cb.fortran_vec (), &matz, jnk.fortran_vec ());
+
+	  F77_FCN (qzit) (&a_rows, &a_rows, ca.fortran_vec (),
+			  cb.fortran_vec (), &eps, &matz,
+			  jnk.fortran_vec (), &info);  
+
+	  if (info)
+	    error ("qzvalue: trouble in qzit, info = %d", info);
+
+	  F77_FCN (qzval) (&a_rows, &a_rows, ca.fortran_vec (),
+			   cb.fortran_vec (), alfr.fortran_vec (),
+			   alfi.fortran_vec (), beta.fortran_vec (),
+			   &matz, jnk.fortran_vec ());
+
+// Count and extract finite generalized eigenvalues.
+
+	  int i, cnt;
+	  Complex Im (0, 1);
+	  for (i = 0, cnt = 0; i < a_rows; i++)
+	    if (beta (i) != 0)
+	      cnt++;
+
+	  ComplexColumnVector cx (cnt, 0.0);
+
+	  for (i = 0; i < a_rows; i++)
+	    {
+	      if (beta (i) != 0)
+		{
+
+// Finite generalized eigenvalue.
+
+		  cnt--;
+		  cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i);
+		}
+	    }
+	  retval[0] = tree_constant (cx);
+	}
+    }
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/