# HG changeset patch # User jwe # Date 745020770 0 # Node ID 2947a1ad8ca3066f9f35d99ed08956f833205d10 # Parent 7541831988d2405fe2c80224f26f3ccaaac547e1 [project @ 1993-08-10 22:12:50 by jwe] Initial revision diff -r 7541831988d2 -r 2947a1ad8ca3 src/syl.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/syl.cc Tue Aug 10 22:12:50 1993 +0000 @@ -0,0 +1,260 @@ +// tc-syl.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 + +#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 (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) (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); + +// Local function: check for empty matrix arguments. Probably should make +// this available elsewhere, since tc-xxx functions do this a lot. + +static inline int +empty_arg (tree_constant& arg) +{ + return (arg.rows () == 0 || arg.columns () == 0); +} + + +// Local function: construct return vector of empty matrices. Return +// empty matrices and/or gripe when appropriate. Probably should make +// this available elsewhere, since tc-xxx functions do this a lot. + +tree_constant * +empty_tree (int nargout, char* fcn_name) +{ + tree_constant *retval = NULL_TREE_CONST; + +// Got an empty argument, check if should gripe/return empty values. + + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) + { + if (flag < 0) + gripe_empty_arg (fcn_name, 0); + + Matrix m; + retval = new tree_constant [nargout+1]; + for (int i = 0; i < nargout; i++) + retval[i] = tree_constant (m); + } + else + gripe_empty_arg (fcn_name, 1); + + return retval; +} + +// Return value of tree_constant argument as ComplexMatrix. + +ComplexMatrix +ComplexMatrixLoad (tree_constant& arg) +{ + ComplexMatrix retval; + +// Set argument size for scalar (for later). + + switch (arg.const_type ()) + { + case tree_constant_rep::scalar_constant: + retval.resize (1, 1); + { + double real_val = arg.double_value (); + retval.elem (0, 0) = real_val; + } + break; + case tree_constant_rep::complex_scalar_constant: + retval.resize (1, 1); + retval.elem (0, 0) = arg.complex_value (); + break; + case tree_constant_rep::matrix_constant: + { + Matrix tmp = arg.matrix_value (); + retval = tmp; + } + break; + case tree_constant_rep::complex_matrix_constant: + retval = arg.complex_matrix_value (); + break; + default: + panic_impossible (); + break; + } + return retval; +} + +#ifdef WITH_DLD +tree_constant * +builtin_syl_2 (tree_constant *args, int nargin, int nargout) +{ + return svd (args, nargin, nargout); +} +#endif + +tree_constant * +syl (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 (); + tree_constant argc = args[3].make_numeric (); + + if (empty_arg (arga) || empty_arg (argb) || empty_arg (argc)) + retval = empty_tree (nargout, "syl"); + 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 (); + int c_rows = argc.rows (); + int c_cols = argc.columns (); + + if ((a_rows != a_cols) || (b_rows != b_cols)) + { + gripe_square_matrix_required ("syl: first two parameters:"); + return retval; + } + else if ((a_rows != c_rows) || (b_rows != c_cols)) + { + 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 () + || argc.is_complex_type ()) + { + +// Do everything in complex arithmetic; + + ComplexMatrix ca = ComplexMatrixLoad (arga); + ComplexMatrix cb = ComplexMatrixLoad (argb); + ComplexMatrix cc = ComplexMatrixLoad (argc); + +// 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 (); + + cx = ua.hermitian () * cc * ub; + +// Solve the sylvester equation, back-transform, and return the solution. + + double scale; + int info; + int one = 1; + + F77_FCN (ztrsyl) ("N", "N", &one, &a_rows, &b_rows, + sch_a.fortran_vec (), &a_rows, + sch_b.fortran_vec (), &b_rows, + cx.fortran_vec (), &a_rows, &scale, &info, + 1L, 1L); + + cx = -ua * cx * ub.hermitian (); + + retval[0] = tree_constant (cx); + } + else + { + +// Do everything in real arithmetic; + + Matrix ca = arga.to_matrix (); + Matrix cb = argb.to_matrix (); + Matrix cc = argc.to_matrix (); + +// 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 (); + + cx = ua.transpose () * cc * ub; + +// Solve the sylvester equation, back-transform, and return the solution. + + double scale; + int info; + int one = 1; + + F77_FCN (dtrsyl) ("N", "N", &one, &a_rows, &b_rows, + sch_a.fortran_vec (), &a_rows, + sch_b.fortran_vec (), &b_rows, + cx.fortran_vec (), &a_rows, &scale, &info, + 1L, 1L); + + if (info) + error ("syl: trouble in dtrsyl info = %d", info); + + cx = -ua*cx*ub.transpose (); + + retval[0] = tree_constant (cx); + } + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/