/* -*- MaTX -*-
 *
 *  NAME
 *      funm() - Matrix functions
 *
 *  SYNOPSIS
 *      F = funm(A, fun)
 *         Matrix F;
 *         Matrix A;
 *         Complex fun();
 *
 *  DESCRIPTION
 *      funm(A, fun) evaluates the matrix function specified by fun().
 *
 *  EXAMPLE
 *      funm(A, exp_f) calculates the matrix exponential, and funm(A, log_f)
 *      the matrix logarithm, and funm(A, sqrt_f) the matrix square root.
 *
 *        Func Complex exp_f(a) {
 *          Complex a;
 *        {
 *          return exp(a);
 *        }
 *
 *        Func Complex log_f(a) {
 *          Complex a;
 *        {
 *          return log(a);
 *        }
 *
 *        Func Complex sqrt_f(a) {
 *          Complex a;
 *        {
 *          return sqrt(a);
 *        }
 *
 *  REFERENCES
 *      [1] Parlett's method. Golub and VanLoan (1983)
 *
 *  SEE ALSO
 *      exp, log, and square
 */

Func Matrix funm(A, fun, tol, ...)
	Matrix A;
	Complex fun();
	Real tol;
{
	Integer i, j, k, n;
	Complex d, s;
	Matrix F, Qr, Tr;
	CoMatrix Q, T, tmp;

	error(nargchk(2, 3, nargs, "funm"));
	
	n = Rows(A);
	{Qr, Tr} = schur(A);
	{Q, T} = rsf2csf(Qr, Tr);
	F = (Z(n),:);
	for (i = 1; i <= n; i++) {
		F(i,i) = fun(T(i,i));
	}

	if (nargs < 3) {
		tol = EPS * norm(T, 1);
	}

	for (k = 1; k <= n-1; k++) {
		for (i = 1; i <= n-k; i++) {
			j = i + k;
			s = T(i,j) * (F(j,j) - F(i,i));
			if (k > 1) {
				tmp = T(i,i+1:j-1)*F(i+1:j-1,j) - F(i,i+1:j-1)*T(i+1:j-1,j);
				s = s + tmp(1);
			}
			d = T(j,j) - T(i,i);

			if (abs(d) <= tol) {
				warning("funm(): Result may be inaccurate.\n");
				d = (tol,0);
			}
			F(i,j) = s / d;
		}
	}

	F = Q * F * Q#;
	return F;
}