changeset 8:a587712dcf5f draft default tip

funm_atom.m: rename fun_atom to funm_atom * funm_atom.m: rename fun_atom to funm_atom.
author Antonio Pino Robles <data.script93@gmail.com>
date Fri, 29 May 2015 09:48:36 +0200
parents e0817ccb2834
children
files funm_files/fun_atom.m funm_files/funm_atom.m
diffstat 2 files changed, 78 insertions(+), 78 deletions(-) [+]
line wrap: on
line diff
--- a/funm_files/fun_atom.m	Tue May 26 18:14:54 2015 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-function [F,n_terms] = fun_atom(T,fun,tol,prnt)
-%FUNM_ATOM  Function of triangular matrix with nearly constant diagonal.
-%           [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function
-%           FUN at the upper triangular matrix T, where T has nearly constant
-%           diagonal.  A Taylor series is used.
-%           FUN(X,K) must return the K'th derivative of
-%           the function represented by FUN evaluated at the vector X.
-%           TOL is a convergence tolerance for the Taylor series,
-%           defaulting to EPS.
-%           If PRNT ~= 0 trace information is printed.
-%           N_TERMS is the number of terms taken in the Taylor series.
-%           N_TERMS  = -1 signals lack of convergence.
-
-if nargin < 3 | isempty(tol), tol = eps; end
-if nargin < 4, prnt = 0; end
-
-if isequal(fun,@fun_log)   % LOG is special case.
-   [F,n_terms]  = logm_isst(T,prnt);
-   return
-end
-
-itmax = 500;
-
-n = length(T);
-if n == 1, F = feval(fun,T,0); n_terms = 1; return, end
-
-lambda = sum(diag(T))/n;
-F = eye(n)*feval(fun,lambda,0);
-f_deriv_max = zeros(itmax+n-1,1);
-N = T - lambda*eye(n);
-mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf );
-
-P = N;
-max_d = 1;
-
-for k = 1:itmax
-
-    f = feval(fun,lambda,k);
-    F_old = F;
-    F = F + P*f;
-    rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf));
-    if prnt
-        fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k));
-        fprintf('  N^k/k! = %7.1e', norm(P,inf));
-        fprintf('  rel_d = %5.0e',rel_diff);
-        fprintf('  abs_d = %5.0e',norm(F - F_old,inf));
-    end
-    P = P*N/(k+1);
-
-    if rel_diff <= tol
-
-      % Approximate the maximum of derivatives in convex set containing
-      % eigenvalues by maximum of derivatives at eigenvalues.
-      for j = max_d:k+n-1
-          f_deriv_max(j) = norm(feval(fun,diag(T),j),inf);
-      end
-      max_d = k+n;
-      omega = 0;
-      for j = 0:n-1
-          omega = max(omega,f_deriv_max(k+j)/factorial(j));
-      end
-
-      trunc = norm(P,inf)*mu*omega;  % norm(F) moved to RHS to avoid / 0.
-      if prnt
-          fprintf('  [trunc,test] = [%5.0e %5.0e]', ...
-                   trunc, tol*norm(F,inf))
-      end
-      if prnt == 5, trunc = 0; end % Force simple stopping test.
-      if trunc <= tol*norm(F,inf)
-         n_terms = k;
-         if prnt, fprintf('\n'), end, return
-      end
-    end
-
-    if prnt, fprintf('\n'), end
-
-end
-n_terms = -1;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/funm_files/funm_atom.m	Fri May 29 09:48:36 2015 +0200
@@ -0,0 +1,78 @@
+function [F,n_terms] = fun_atom(T,fun,tol,prnt)
+%FUNM_ATOM  Function of triangular matrix with nearly constant diagonal.
+%           [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function
+%           FUN at the upper triangular matrix T, where T has nearly constant
+%           diagonal.  A Taylor series is used.
+%           FUN(X,K) must return the K'th derivative of
+%           the function represented by FUN evaluated at the vector X.
+%           TOL is a convergence tolerance for the Taylor series,
+%           defaulting to EPS.
+%           If PRNT ~= 0 trace information is printed.
+%           N_TERMS is the number of terms taken in the Taylor series.
+%           N_TERMS  = -1 signals lack of convergence.
+
+if nargin < 3 | isempty(tol), tol = eps; end
+if nargin < 4, prnt = 0; end
+
+if isequal(fun,@fun_log)   % LOG is special case.
+   [F,n_terms]  = logm_isst(T,prnt);
+   return
+end
+
+itmax = 500;
+
+n = length(T);
+if n == 1, F = feval(fun,T,0); n_terms = 1; return, end
+
+lambda = sum(diag(T))/n;
+F = eye(n)*feval(fun,lambda,0);
+f_deriv_max = zeros(itmax+n-1,1);
+N = T - lambda*eye(n);
+mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf );
+
+P = N;
+max_d = 1;
+
+for k = 1:itmax
+
+    f = feval(fun,lambda,k);
+    F_old = F;
+    F = F + P*f;
+    rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf));
+    if prnt
+        fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k));
+        fprintf('  N^k/k! = %7.1e', norm(P,inf));
+        fprintf('  rel_d = %5.0e',rel_diff);
+        fprintf('  abs_d = %5.0e',norm(F - F_old,inf));
+    end
+    P = P*N/(k+1);
+
+    if rel_diff <= tol
+
+      % Approximate the maximum of derivatives in convex set containing
+      % eigenvalues by maximum of derivatives at eigenvalues.
+      for j = max_d:k+n-1
+          f_deriv_max(j) = norm(feval(fun,diag(T),j),inf);
+      end
+      max_d = k+n;
+      omega = 0;
+      for j = 0:n-1
+          omega = max(omega,f_deriv_max(k+j)/factorial(j));
+      end
+
+      trunc = norm(P,inf)*mu*omega;  % norm(F) moved to RHS to avoid / 0.
+      if prnt
+          fprintf('  [trunc,test] = [%5.0e %5.0e]', ...
+                   trunc, tol*norm(F,inf))
+      end
+      if prnt == 5, trunc = 0; end % Force simple stopping test.
+      if trunc <= tol*norm(F,inf)
+         n_terms = k;
+         if prnt, fprintf('\n'), end, return
+      end
+    end
+
+    if prnt, fprintf('\n'), end
+
+end
+n_terms = -1;