From octave-maintainers-request at bevo dot che dot wisc dot edu  Thu Dec  2 12:32:19 1999
Subject: polyvalm.m
From: "Ross A. Lippert" <ripper at cs dot sandia dot gov>
To: "octave-maintainers at bevo dot che dot wisc dot edu"   <octave-maintainers@bevo.che.wisc.edu>
Date: Thu, 02 Dec 1999 11:31:46 -0700

This is a multi-part message in MIME format.
--------------D6135A5C2ADB95D32687BC91
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit


Here is a replacement for polyvalm.m which will give correct
answers for matrices which are non-diagonalizable.

The current version of polyvalm.m assumes the matrix has a 
trivial eigenstructure.



-r

P.S. the file was produced by consulting polyvalm.m and polyval.m
and nothing else.
--------------D6135A5C2ADB95D32687BC91
Content-Type: text/plain; charset=us-ascii;
 name="polyvalm.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="polyvalm.m"

## usage: polyvalm (c, x)
##
## Evaluate a polynomial in the matrix sense.
##
## In octave, a polynomial is represented by it's coefficients (arranged
## in descending order). For example a vector c of length n+1 corresponds
## to the following nth order polynomial
##
##   p(x) = c(1) x^n + ... + c(n) x + c(n+1).
##
## polyvalm(c,X) will evaluate the polynomial in the matrix sense, i.e. matrix
## multiplication is used instead of element by element multiplication as is
## used in polyval.
##
## X must be a square matrix.
##
## SEE ALSO: polyval, poly, roots, conv, deconv, residue, filter,
##           polyderiv, polyinteg

## Author: Tony Richardson <arichard at stark dot cc dot oh dot us>
## Created: June 1994
## Adapted-By: jwe
## Changed-By: Ross Lippert <ripper at cs,sandia dot gov>

function y = polyvalm (c, x)

  if (nargin != 2)
    usage ("polyvalm (c, x)");
  endif

  if (! (is_vector (c) || isempty (c)))
    error ("polyvalm: first argument must be a vector.");
  endif

  if (! is_square (x))
    error("polyvalm: second argument must be a square matrix.");
  endif

  if (isempty (c))
    y = [];
    return;
  endif

  n = length(c);
  I = eye(rows(x),columns(x));
  y = c(1) * I
  for index = 2:n,
     y = c(index)*I + x*y;
  endfor

  if (is_symmetric (x))
    y = (y+y')/2;
  endif

endfunction

--------------D6135A5C2ADB95D32687BC91--