Mercurial > hg > octave-max
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;