changeset 1819:8b8498bf8ec5

[project @ 1996-01-31 11:29:17 by jwe]
author jwe
date Wed, 31 Jan 1996 11:31:00 +0000
parents 947837af62ac
children c6fdd23c0e79
files liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/dMatrix.cc liboctave/dMatrix.h src/expm.cc src/givens.cc src/syl.cc
diffstat 7 files changed, 463 insertions(+), 423 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -35,7 +35,9 @@
 
 #include <sys/types.h> // XXX FIXME XXX
 
+#include "CmplxAEPBAL.h"
 #include "CmplxDET.h"
+#include "CmplxSCHUR.h"
 #include "CmplxSVD.h"
 #include "f77-uscore.h"
 #include "lo-error.h"
@@ -78,6 +80,20 @@
   int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);
 
   int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);
+
+  int F77_FCN (zlartg, ZLARTG) (const Complex&, const Complex&,
+				double&, Complex&, Complex&);
+
+  int F77_FCN (ztrsyl, ZTRSYL) (const char*, const char*, const int&,
+				const int&, const int&,
+				const Complex*, const int&,
+				const Complex*, const int&, 
+				const Complex*, const int&, double&,
+				int&, long, long);
+
+  double F77_FCN (zlange, ZLANGE) (const char*, const int&,
+				   const int&, const Complex*,
+				   const int&, double*); 
 }
 
 // Complex Matrix class
@@ -1385,6 +1401,124 @@
   return retval;
 }
 
+// Constants for matrix exponential calculation.
+
+static double padec [] =
+{
+  5.0000000000000000e-1,
+  1.1666666666666667e-1,
+  1.6666666666666667e-2,
+  1.6025641025641026e-3,
+  1.0683760683760684e-4,
+  4.8562548562548563e-6,
+  1.3875013875013875e-7,
+  1.9270852604185938e-9,
+};
+
+ComplexMatrix
+ComplexMatrix::expm (void) const
+{
+  ComplexMatrix retval;
+
+  ComplexMatrix m = *this;
+
+  int nc = columns ();
+
+  // trace shift value
+  Complex trshift = 0.0;
+
+  // Preconditioning step 1: trace normalization.
+
+  for (int i = 0; i < nc; i++)
+    trshift += m.elem (i, i);
+
+  trshift /= nc;
+
+  for (int i = 0; i < nc; i++)
+    m.elem (i, i) -= trshift;
+
+  // Preconditioning step 2: eigenvalue balancing.
+
+  ComplexAEPBALANCE mbal (m, "B");
+  m = mbal.balanced_matrix ();
+  ComplexMatrix d = mbal.balancing_matrix ();
+
+  // Preconditioning step 3: scaling.
+
+  ColumnVector work (nc);
+  double inf_norm
+    = F77_FCN (zlange, ZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
+				work.fortran_vec ());
+
+  int sqpow = (int) (inf_norm > 0.0
+		     ? (1.0 + log (inf_norm) / log (2.0))
+		     : 0.0);
+
+  // Check whether we need to square at all.
+
+  if (sqpow < 0)
+    sqpow = 0;
+
+  if (sqpow > 0)
+    {
+      double scale_factor = 1.0;
+      for (int i = 0; i < sqpow; i++)
+	scale_factor *= 2.0;
+
+      m = m / scale_factor;
+    }
+
+  // npp, dpp: pade' approx polynomial matrices.
+
+  ComplexMatrix npp (nc, nc, 0.0);
+  ComplexMatrix dpp = npp;
+
+  // Now powers a^8 ... a^1.
+
+  int minus_one_j = -1;
+  for (int j = 7; j >= 0; j--)
+    {
+      npp = m * npp + m * padec[j];
+      dpp = m * dpp + m * (minus_one_j * padec[j]);
+      minus_one_j *= -1;
+    }
+
+  // Zero power.
+
+  dpp = -dpp;
+  for (int j = 0; j < nc; j++)
+    {
+      npp.elem (j, j) += 1.0;
+      dpp.elem (j, j) += 1.0;
+    }
+
+  // Compute pade approximation = inverse (dpp) * npp.
+
+  retval = dpp.solve (npp);
+	
+  // Reverse preconditioning step 3: repeated squaring.
+
+  while (sqpow)
+    {
+      retval = retval * retval;
+      sqpow--;
+    }
+
+  // Reverse preconditioning step 2: inverse balancing.
+  // XXX FIXME XXX -- should probably do this with Lapack calls
+  // instead of a complete matrix inversion.
+
+  retval = retval.transpose ();
+  d = d.transpose ();
+  retval = retval * d;
+  retval = d.solve (retval);
+  retval = retval.transpose ();
+
+  // Reverse preconditioning step 1: fix trace normalization.
+
+  return retval * exp (trshift);
+}
+
 // column vector by row vector -> matrix operations
 
 ComplexMatrix
@@ -3372,6 +3506,70 @@
   return is;
 }
 
+ComplexMatrix
+Givens (const Complex& x, const Complex& y)
+{
+  double cc;
+  Complex cs, temp_r;
+ 
+  F77_FCN (zlartg, ZLARTG) (x, y, cc, cs, temp_r);
+
+  ComplexMatrix g (2, 2);
+
+  g.elem (0, 0) = cc;
+  g.elem (1, 1) = cc;
+  g.elem (0, 1) = cs;
+  g.elem (1, 0) = -conj (cs);
+
+  return g;
+}
+
+ComplexMatrix
+Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
+	   const ComplexMatrix& c)
+{
+  ComplexMatrix retval;
+
+  // XXX FIXME XXX -- need to check that a, b, and c are all the same
+  // size.
+
+  // Compute Schur decompositions
+
+  ComplexSCHUR as (a, "U");
+  ComplexSCHUR bs (b, "U");
+  
+  // Transform c to new coordinates.
+
+  ComplexMatrix ua = as.unitary_matrix ();
+  ComplexMatrix sch_a = as.schur_matrix ();
+
+  ComplexMatrix ub = bs.unitary_matrix ();
+  ComplexMatrix sch_b = bs.schur_matrix ();
+  
+  ComplexMatrix cx = ua.hermitian () * c * ub;
+
+  // Solve the sylvester equation, back-transform, and return the
+  // solution.
+
+  int a_nr = a.rows ();
+  int b_nr = b.rows ();
+
+  double scale;
+  int info;
+  
+  F77_FCN (ztrsyl, ZTRSYL) ("N", "N", 1, a_nr, b_nr,
+			    sch_a.fortran_vec (), a_nr,
+			    sch_b.fortran_vec (), b_nr,
+			    cx.fortran_vec (), a_nr, scale,
+			    info, 1L, 1L);
+
+  // XXX FIXME XXX -- check info?
+
+  retval = -ua * cx * ub.hermitian ();
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -161,6 +161,8 @@
   ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
 			       int& rank) const;
 
+  ComplexMatrix expm (void) const;
+
   // column vector by row vector -> matrix operations
 
   friend ComplexMatrix operator * (const ColumnVector& a,
@@ -348,6 +350,11 @@
   ComplexMatrix (Complex *d, int r, int c) : MArray2<Complex> (d, r, c) { }
 };
 
+ComplexMatrix Givens (const Complex&, const Complex&);
+
+ComplexMatrix Sylvester (const ComplexMatrix&, const ComplexMatrix&,
+			 const ComplexMatrix&);
+
 #endif
 
 /*
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -37,7 +37,9 @@
 
 #include <sys/types.h>  // XXX FIXME XXX
 
+#include "dbleAEPBAL.h"
 #include "dbleDET.h"
+#include "dbleSCHUR.h"
 #include "dbleSVD.h"
 #include "f77-uscore.h"
 #include "lo-error.h"
@@ -81,6 +83,19 @@
   int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);
 
   int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);
+
+  int F77_FCN (dlartg, DLARTG) (const double&, const double&, double&,
+				double&, double&);
+
+  int F77_FCN (dtrsyl, DTRSYL) (const char*, const char*, const int&,
+				const int&, const int&, const double*,
+				const int&, const double*, const int&,
+				const double*, const int&, double&,
+				int&, long, long);
+
+  double F77_FCN (dlange, DLANGE) (const char*, const int&,
+				   const int&, const double*,
+				   const int&, double*); 
 }
 
 // Matrix class.
@@ -1158,6 +1173,122 @@
   return tmp.lssolve (b, info, rank);
 }
 
+// Constants for matrix exponential calculation.
+
+static double padec [] =
+{
+  5.0000000000000000e-1,
+  1.1666666666666667e-1,
+  1.6666666666666667e-2,
+  1.6025641025641026e-3,
+  1.0683760683760684e-4,
+  4.8562548562548563e-6,
+  1.3875013875013875e-7,
+  1.9270852604185938e-9,
+};
+
+Matrix
+Matrix::expm (void) const
+{
+  Matrix retval;
+
+  Matrix m = *this;
+
+  int nc = columns ();
+
+  // trace shift value
+  double trshift = 0;
+
+  // Preconditioning step 1: trace normalization.
+
+  for (int i = 0; i < nc; i++)
+    trshift += m.elem (i, i);
+
+  trshift /= nc;
+
+  for (int i = 0; i < nc; i++)
+    m.elem (i, i) -= trshift;
+
+  // Preconditioning step 2: balancing.
+
+  AEPBALANCE mbal (m, "B");
+  m = mbal.balanced_matrix ();
+  Matrix d = mbal.balancing_matrix ();
+
+  // Preconditioning step 3: scaling.
+
+  ColumnVector work(nc);
+  double inf_norm
+    = F77_FCN (dlange, DLANGE) ("I", nc, nc, m.fortran_vec (),nc,
+				work.fortran_vec ());
+
+  int sqpow = (int) (inf_norm > 0.0
+		     ? (1.0 + log (inf_norm) / log (2.0))
+		     : 0.0);
+
+  // Check whether we need to square at all.
+
+  if (sqpow < 0)
+    sqpow = 0;
+
+  if (sqpow > 0)
+    {
+      double scale_factor = 1.0;
+      for (int i = 0; i < sqpow; i++)
+	scale_factor *= 2.0;
+
+      m = m / scale_factor;
+    }
+
+  // npp, dpp: pade' approx polynomial matrices.
+
+  Matrix npp (nc, nc, 0.0);
+  Matrix dpp = npp;
+
+  // Now powers a^8 ... a^1.
+
+  int minus_one_j = -1;
+  for (int j = 7; j >= 0; j--)
+    {
+      npp = m * npp + m * padec[j];
+      dpp = m * dpp + m * (minus_one_j * padec[j]);
+      minus_one_j *= -1;
+    }
+
+  // Zero power.
+
+  dpp = -dpp;
+  for(int j = 0; j < nc; j++)
+    {
+      npp.elem (j, j) += 1.0;
+      dpp.elem (j, j) += 1.0;
+    }
+
+  // Compute pade approximation = inverse (dpp) * npp.
+
+  retval = dpp.solve (npp);
+
+  // Reverse preconditioning step 3: repeated squaring.
+
+  while (sqpow)
+    {
+      retval = retval * retval;
+      sqpow--;
+    }
+
+  // Reverse preconditioning step 2: inverse balancing.
+
+  retval = retval.transpose();
+  d = d.transpose ();
+  retval = retval * d;
+  retval = d.solve (retval);
+  retval = retval.transpose ();
+
+  // Reverse preconditioning step 1: fix trace normalization.
+
+  return retval * exp (trshift);
+}
+
 Matrix&
 Matrix::operator += (const Matrix& a)
 {
@@ -2357,6 +2488,69 @@
   return count;
 }
 
+Matrix
+Givens (double x, double y)
+{
+  double cc, s, temp_r;
+
+  F77_FCN (dlartg, DLARTG) (x, y, cc, s, temp_r);
+
+  Matrix g (2, 2);
+
+  g.elem (0, 0) = cc;
+  g.elem (1, 1) = cc;
+  g.elem (0, 1) = s;
+  g.elem (1, 0) = -s;
+
+  return g;
+}
+
+Matrix
+Sylvester (const Matrix& a, const Matrix& b, const Matrix& c)
+{
+  Matrix retval;
+
+  // XXX FIXME XXX -- need to check that a, b, and c are all the same
+  // size.
+
+  // Compute Schur decompositions.
+
+  SCHUR as (a, "U");
+  SCHUR bs (b, "U");
+  
+  // Transform c to new coordinates.
+
+  Matrix ua = as.unitary_matrix ();
+  Matrix sch_a = as.schur_matrix ();
+
+  Matrix ub = bs.unitary_matrix ();
+  Matrix sch_b = bs.schur_matrix ();
+  
+  Matrix cx = ua.transpose () * c * ub;
+  
+  // Solve the sylvester equation, back-transform, and return the
+  // solution.
+
+  int a_nr = a.rows ();
+  int b_nr = b.rows ();
+
+  double scale;
+  int info;
+
+  F77_FCN (dtrsyl, DTRSYL) ("N", "N", 1, a_nr, b_nr,
+			    sch_a.fortran_vec (), a_nr, 
+			    sch_b.fortran_vec (), b_nr,
+			    cx.fortran_vec (), a_nr, scale,
+			    info, 1L, 1L);
+
+
+  // XXX FIXME XXX -- check info?
+  
+  retval = -ua*cx*ub.transpose ();
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -157,6 +157,8 @@
   ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
 			       int& rank) const;
 
+  Matrix expm (void) const;
+
   Matrix& operator += (const Matrix& a);
   Matrix& operator -= (const Matrix& a);
 
@@ -240,6 +242,10 @@
   Matrix (double *d, int r, int c) : MArray2<double> (d, r, c) { }
 };
 
+Matrix Givens (double, double);
+
+Matrix Sylvester (const Matrix&, const Matrix&, const Matrix&);
+
 #endif
 
 /*
--- a/src/expm.cc
+++ b/src/expm.cc
@@ -27,30 +27,13 @@
 #include <config.h>
 #endif
 
-#include "CmplxAEPBAL.h"
-#include "dbleAEPBAL.h"
-#include "f77-uscore.h"
-#include "oct-math.h"
-
 #include "defun-dld.h"
 #include "error.h"
 #include "gripes.h"
 #include "help.h"
 #include "oct-obj.h"
-#include "user-prefs.h"
 #include "utils.h"
 
-extern "C"
-{
-  double F77_FCN (dlange, DLANGE) (const char*, const int&,
-				   const int&, const double*,
-				   const int&, double*); 
-
-  double F77_FCN (zlange, ZLANGE) (const char*, const int&,
-				   const int&, const Complex*,
-				   const int&, double*); 
-}
-
 DEFUN_DLD_BUILTIN ("expm", Fexpm, Sexpm, FSexpm, 10,
   "expm (X): matrix exponential, e^A")
 {
@@ -66,20 +49,6 @@
 
   tree_constant arg = args(0);
 
-  // Constants for matrix exponential calculation.
-
-  static double padec [] =
-    {
-      5.0000000000000000e-1,
-      1.1666666666666667e-1,
-      1.6666666666666667e-2,
-      1.6025641025641026e-3,
-      1.0683760683760684e-4,
-      4.8562548562548563e-6,
-      1.3875013875013875e-7,
-      1.9270852604185938e-9,
-    };
-
   int nr = arg.rows ();
   int nc = arg.columns ();
 
@@ -96,109 +65,14 @@
       return retval;
     }
 
-  char* balance_job = "B";	// variables for balancing
-
   if (arg.is_real_type ())
     {
-      // Compute the exponential.
-
       Matrix m = arg.matrix_value ();
 
       if (error_state)
 	return retval;
-
-      double trshift = 0;		// trace shift value
-
-      // Preconditioning step 1: trace normalization.
-
-      for (int i = 0; i < nc; i++)
-	trshift += m.elem (i, i);
-
-      trshift /= nc;
-
-      for (int i = 0; i < nc; i++)
-	m.elem (i, i) -= trshift;
-
-      // Preconditioning step 2: balancing.
-
-      AEPBALANCE mbal (m, balance_job);
-      m = mbal.balanced_matrix ();
-      Matrix d = mbal.balancing_matrix ();
-
-      // Preconditioning step 3: scaling.
-
-      ColumnVector work(nc);
-      double inf_norm
-	= F77_FCN (dlange, DLANGE) ("I", nc, nc, m.fortran_vec (),nc,
-				    work.fortran_vec ());
-
-      int sqpow = (int) (inf_norm > 0.0
-			 ? (1.0 + log (inf_norm) / log (2.0))
-			 : 0.0);
-
-      // Check whether we need to square at all.
-
-      if (sqpow < 0)
-	sqpow = 0;
-
-      if (sqpow > 0)
-	{
-	  double scale_factor = 1.0;
-	  for (int i = 0; i < sqpow; i++)
-	    scale_factor *= 2.0;
-
-	  m = m / scale_factor;
-	}
-
-      // npp, dpp: pade' approx polynomial matrices.
-
-      Matrix npp (nc, nc, 0.0);
-      Matrix dpp = npp;
-
-      // Now powers a^8 ... a^1.
-
-      int minus_one_j = -1;
-      for (int j = 7; j >= 0; j--)
-	{
-	  npp = m * npp + m * padec[j];
-	  dpp = m * dpp + m * (minus_one_j * padec[j]);
-	  minus_one_j *= -1;
-	}
-
-      // Zero power.
-
-      dpp = -dpp;
-      for(int j = 0; j < nc; j++)
-	{
-	  npp.elem (j, j) += 1.0;
-	  dpp.elem (j, j) += 1.0;
-	}
-
-      // Compute pade approximation = inverse (dpp) * npp.
-
-      Matrix result = dpp.solve (npp);
-
-      // Reverse preconditioning step 3: repeated squaring.
-
-      while (sqpow)
-	{
-	  result = result * result;
-	  sqpow--;
-	}
-
-      // Reverse preconditioning step 2: inverse balancing.
-
-      result = result.transpose();
-      d = d.transpose ();
-      result = result * d;
-      result = d.solve (result);
-      result = result.transpose ();
-
-      // Reverse preconditioning step 1: fix trace normalization.
-
-      result = result * exp (trshift);
-
-      retval = result;
+      else
+	retval = m.expm ();
     }
   else if (arg.is_complex_type ())
     {
@@ -206,101 +80,8 @@
 
       if (error_state)
 	return retval;
-
-      Complex trshift = 0.0;		// trace shift value
-
-      // Preconditioning step 1: trace normalization.
-
-      for (int i = 0; i < nc; i++)
-	trshift += m.elem (i, i);
-
-      trshift /= nc;
-
-      for (int i = 0; i < nc; i++)
-	m.elem (i, i) -= trshift;
-
-      // Preconditioning step 2: eigenvalue balancing.
-
-      ComplexAEPBALANCE mbal (m, balance_job);
-      m = mbal.balanced_matrix ();
-      ComplexMatrix d = mbal.balancing_matrix ();
-
-      // Preconditioning step 3: scaling.
-
-      ColumnVector work (nc);
-      double inf_norm
-	= F77_FCN (zlange, ZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
-				    work.fortran_vec ());
-
-      int sqpow = (int) (inf_norm > 0.0
-			 ? (1.0 + log (inf_norm) / log (2.0))
-			 : 0.0);
-
-      // Check whether we need to square at all.
-
-      if (sqpow < 0)
-	sqpow = 0;
-
-      if (sqpow > 0)
-	{
-	  double scale_factor = 1.0;
-	  for (int i = 0; i < sqpow; i++)
-	    scale_factor *= 2.0;
-
-	  m = m / scale_factor;
-	}
-
-      // npp, dpp: pade' approx polynomial matrices.
-
-      ComplexMatrix npp (nc, nc, 0.0);
-      ComplexMatrix dpp = npp;
-
-      // Now powers a^8 ... a^1.
-
-      int minus_one_j = -1;
-      for (int j = 7; j >= 0; j--)
-	{
-	  npp = m * npp + m * padec[j];
-	  dpp = m * dpp + m * (minus_one_j * padec[j]);
-	  minus_one_j *= -1;
-	}
-
-      // Zero power.
-
-      dpp = -dpp;
-      for (int j = 0; j < nc; j++)
-	{
-	  npp.elem (j, j) += 1.0;
-	  dpp.elem (j, j) += 1.0;
-	}
-
-      // Compute pade approximation = inverse (dpp) * npp.
-
-      ComplexMatrix result = dpp.solve (npp);
-	
-      // Reverse preconditioning step 3: repeated squaring.
-
-      while (sqpow)
-	{
-	  result = result * result;
-	  sqpow--;
-	}
-
-      // Reverse preconditioning step 2: inverse balancing.
-      // XXX FIXME XXX -- should probably do this with Lapack calls
-      // instead of a complete matrix inversion.
-
-      result = result.transpose ();
-      d = d.transpose ();
-      result = result * d;
-      result = d.solve (result);
-      result = result.transpose ();
-
-      // Reverse preconditioning step 1: fix trace normalization.
-
-      result = result * exp (trshift);
-
-      retval = result;
+      else
+	retval = m.expm ();
     }
   else
     {
--- a/src/givens.cc
+++ b/src/givens.cc
@@ -21,29 +21,16 @@
 
 */
 
-// Written by A. S. Hodel <scotte@eng.auburn.edu>
+// Originally written by A. S. Hodel <scotte@eng.auburn.edu>
 
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
-#include "f77-uscore.h"
-
 #include "defun-dld.h"
 #include "error.h"
-#include "gripes.h"
 #include "help.h"
 #include "oct-obj.h"
-#include "user-prefs.h"
-
-extern "C"
-{
-  int F77_FCN (dlartg, DLARTG) (const double&, const double&, double&,
-				double&, double&);
-
-  int F77_FCN (zlartg, ZLARTG) (const Complex&, const Complex&,
-				double&, Complex&, Complex&);
-}
 
 DEFUN_DLD_BUILTIN ("givens", Fgivens, Sgivens, FSgivens, 11,
   "G = givens (X, Y)\n\
@@ -62,120 +49,67 @@
       print_usage ("givens");
       return retval;
     }
-
-  tree_constant arg_a = args(0);
-  tree_constant arg_b = args(1);
-
-  if (! arg_a.is_scalar_type () && arg_b.is_scalar_type ())
-    {
-      error("givens: requires two scalar arguments");
-      return retval;
-    }
-
-  Complex cx, cy;
-  double x, y;
-
-  if (arg_a.is_complex_type ())
-    {
-      cx = arg_a.complex_value ();
-
-      if (error_state)
-	return retval;
-    }
-  else 
-    {
-      x = arg_a.double_value ();
-
-      if (error_state)
-	return retval;
-
-      // Convert to complex just in case...
-
-      cx = x;
-    }
-
-  if (arg_b.is_complex_type ())
-    {
-      cy = arg_b.complex_value ();
-
-      if (error_state)
-	return retval;
-    }
   else
     {
-      y = arg_b.double_value ();
-
-      if (error_state)
-	return retval;
+      if (args(0).is_complex_type () || args(1).is_complex_type ())
+	{
+	  Complex cx = args(0).complex_value ();
+	  Complex cy = args(1).complex_value ();
 
-      // Convert to complex just in case...
-
-      cy = y;
-    }
-
-  // Now compute the rotation.
+	  if (! error_state)
+	    {
+	      ComplexMatrix result = Givens (cx, cy);
 
-  double cc;
-  if (arg_a.is_complex_type () || arg_b.is_complex_type ())
-    {
-      Complex cs, temp_r;
- 
-      F77_FCN (zlartg, ZLARTG) (cx, cy, cc, cs, temp_r);
-
-      switch (nargout)
-	{
-	case 0:
-	case 1:
-	  {
-	    ComplexMatrix g (2, 2);
-	    g.elem (0, 0) = cc;
-	    g.elem (1, 1) = cc;
-	    g.elem (0, 1) = cs;
-	    g.elem (1, 0) = -conj (cs);
-
-	    retval(0) = g;
-	  }
-	  break;
+	      if (! error_state)
+		{
+		  switch (nargout)
+		    {
+		    case 0:
+		    case 1:
+		      retval(0) = result;
+		      break;
    
-	case 2:
-	  retval(0) = cc;
-	  retval(1) = cs;
-	  break;
+		    case 2:
+		      retval(1) = result (0, 1);
+		      retval(0) = result (0, 0);
+		      break;
 
-	default:  
-	  error ("givens: invalid number of output arguments");
-	  break;
+		    default:
+		      error ("givens: invalid number of output arguments");
+		      break;
+		    }
+		}
+	    }
 	}
-    }
-  else
-    {
-      double s, temp_r;
-
-      F77_FCN (dlartg, DLARTG) (x, y, cc, s, temp_r);
-
-      switch (nargout)
+      else
 	{
-	case 0:
-	case 1:
-	  {
-	    Matrix g (2, 2);
-	    g.elem (0, 0) = cc;
-	    g.elem (1, 1) = cc;
-	    g.elem (0, 1) = s;
-	    g.elem (1, 0) = -s;
+	  double x = args(0).double_value ();
+	  double y = args(1).double_value ();
+
+	  if (! error_state)
+	    {
+	      Matrix result = Givens (x, y);
 
-	    retval(0) = g;
-	  }
-	  break;
+	      if (! error_state)
+		{
+		  switch (nargout)
+		    {
+		    case 0:
+		    case 1:
+		      retval(0) = result;
+		      break;
    
-	case 2:
-	  retval(0) = cc;
-	  retval(1) = s;
-	  break;
-   
-	default:
-	  error ("givens: invalid number of output arguments");
-	  break;
+		    case 2:
+		      retval(1) = result (0, 1);
+		      retval(0) = result (0, 0);
+		      break;
+
+		    default:
+		      error ("givens: invalid number of output arguments");
+		      break;
+		    }
+		}
+	    }
 	}
     }
 
--- a/src/syl.cc
+++ b/src/syl.cc
@@ -27,34 +27,13 @@
 #include <config.h>
 #endif
 
-#include "CmplxSCHUR.h"
-#include "dbleSCHUR.h"
-#include "f77-uscore.h"
-
 #include "defun-dld.h"
 #include "error.h"
 #include "gripes.h"
 #include "help.h"
 #include "oct-obj.h"
-#include "user-prefs.h"
 #include "utils.h"
 
-extern "C"
-{
-  int F77_FCN (dtrsyl, DTRSYL) (const char*, const char*, const int&,
-				const int&, const int&, const double*,
-				const int&, const double*, const int&,
-				const double*, const int&, double&,
-				int&, long, long);
- 
-  int F77_FCN (ztrsyl, ZTRSYL) (const char*, const char*, const int&,
-				const int&, const int&,
-				const Complex*, const int&,
-				const Complex*, const int&, 
-				const Complex*, const int&, double&,
-				int&, long, long);
-}
-
 DEFUN_DLD_BUILTIN ("syl", Fsyl, Ssyl, FSsyl, 11,
   "X = syl (A, B, C): solve the Sylvester equation A X + X B + C = 0")
 {
@@ -126,35 +105,7 @@
 	if (error_state)
 	  return retval;
 
-	// Compute Schur decompositions
-
-	ComplexSCHUR as (ca, "U");
-	ComplexSCHUR bs (cb, "U");
-  
-	// Transform cc to new coordinates.
-
-	ComplexMatrix ua = as.unitary_matrix ();
-	ComplexMatrix sch_a = as.schur_matrix ();
-	ComplexMatrix ub = bs.unitary_matrix ();
-	ComplexMatrix sch_b = bs.schur_matrix ();
-  
-	ComplexMatrix cx = ua.hermitian () * cc * ub;
-  
-	// Solve the sylvester equation, back-transform, and return
-	// the solution.
-  
-	double scale;
-	int info;
-  
-	F77_FCN (ztrsyl, ZTRSYL) ("N", "N", 1, a_nr, b_nr,
-				  sch_a.fortran_vec (), a_nr,
-				  sch_b.fortran_vec (), b_nr,
-				  cx.fortran_vec (), a_nr, scale,
-				  info, 1L, 1L);
-
-	cx = -ua * cx * ub.hermitian ();
-  
-	retval = cx;
+	retval = Sylvester (ca, cb, cc);
       }
     else
       {
@@ -175,38 +126,7 @@
 	if (error_state)
 	  return retval;
 
-	// Compute Schur decompositions.
-
-	SCHUR as (ca, "U");
-	SCHUR bs (cb, "U");
-  
-	// Transform cc to new coordinates.
-
-	Matrix ua = as.unitary_matrix ();
-	Matrix sch_a = as.schur_matrix ();
-	Matrix ub = bs.unitary_matrix ();
-	Matrix sch_b = bs.schur_matrix ();
-  
-	Matrix cx = ua.transpose () * cc * ub;
-  
-	// Solve the sylvester equation, back-transform, and return
-	// the solution.
-  
-	double scale;
-	int info;
-
-	F77_FCN (dtrsyl, DTRSYL) ("N", "N", 1, a_nr, b_nr,
-				  sch_a.fortran_vec (), a_nr, 
-				  sch_b.fortran_vec (), b_nr,
-				  cx.fortran_vec (), a_nr, scale,
-				  info, 1L, 1L);
-
-	if (info)
-	  error ("syl: trouble in dtrsyl info = %d", info);
-  
-	cx = -ua*cx*ub.transpose ();
-  
-	retval = cx;
+	retval = Sylvester (ca, cb, cc);
       }
 
   return retval;