## Copyright (C) 2000 P.R. Nienhuis
##
## This program 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.
##
## This program 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 this program; see the file COPYING.  If not, write to the
## Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## funm:  Matrix equivalent of function 'name'
##
## Usage:    B = funm(A, name)
##  where    A = square non-singular matrix, provisionally
##               real-valued
##           B = square result matrix
##        name = string, name of function to apply to A.
##        args = any arguments to pass to function 'name'
##               The function must accept a vector and apply
##               element-wise to that vector.
##
## Example:  To compute sqrtm(A), you could use funm(A, 'sqrt')

## NOTE: the following comments are withheld until they can be verified
##
## If you have a network of coupled systems, where for the individual
## systems a solution exists in terms of scalar variables, in many
## cases the network might be solved using the same form of the
## solution but with substituting the matrix equivalent of the function
## applied to the scalar variables.
## The approach is to do an eigen-analysis of the network to decouple
## the systems, apply the scalar functions to the eigenvalues,
## and then recombine the systems into a network.

## Author: P.R. Nienhuis, 106130.1515@compuserve.com
## Additions by P. Kienzle, .........

function B = funm(A, name, ...)

  if (nargin < 2)
    usage ("B = funm (A, 'f' [, args])");
  endif

  dfi = do_fortran_indexing;
  unwind_protect
    do_fortran_indexing = 1;
    n = size (A, 1);
    idx = 1 : n+1 : n*n;
    [V, D] = eig (A);
    D (idx) = feval (name, D (idx), all_va_args);
    B = V * D * inv (V);
  unwind_protect_cleanup
    do_fortran_indexing = dfi;
  end_unwind_protect

endfunction