## Copyright (C) 2010 Marco Caliari ## ## 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 3 of the License, 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, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{U}, @var{T}] =} rsf2csf (@var{UR}, @var{TR}) ## Converts a real, upper quasi-triangular Schur form @var{TR} to a complex, ## upper triangular Schur form @var{T}. ## ## Note that the following relations hold: ## ## @code{@var{UR} * @var{TR} * @var{UR}' = @var{U} * @var{T} * @var{U}'} and @code{@var{U}' * ## @var{U}} is identity. ## ## Note also that U and T are not unique. ## ## @end deftypefn function [U, T] = rsf2csf (U, T) n = length (T); if (n > 1) G = speye (n); blocks = find(diag(T,-1))'+1; for m = blocks s = sqrt(T(m,m-1)/(T(m-1,m)-T(m,m-1))); c = sqrt(T(m-1,m)/(T(m-1,m)-T(m,m-1))); G (m-1:m, m-1:m) = [c, s;-s', c]; endfor T = G' * T * G; ## Replace almost-zeros on the subdiagonal with zeros T (sub2ind ([n, n], blocks, blocks - 1)) = 0; U = U * G; endif endfunction %!test %! A = [1, 1, 1, 2; 1, 2, 1, 1; 1, 1, 3, 1; -2, 1, 1, 1]; %! [u, t] = schur (A); %! [U, T] = rsf2csf (u, t); %! assert (norm (u * t * u' - U * T * U'), 0, 1e-12) %! assert (norm (A - U * T * U'), 0, 1e-12) %!test %! A = rand (10); %! [u, t] = schur (A); %! [U, T] = rsf2csf (u, t); %! assert (norm (tril (T, -1)), 0) %! assert (norm (U * U'), 1, 1e-14)