Mercurial > matrix-functions

comparison mftoolbox/funm_condest1.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [c,est] = funm_condest1(A,fun,fun_frechet,flag1,varargin)
2 %FUNM_CONDEST1 Estimate of 1-norm condition number of matrix function.
3 % C = FUNM_CONDEST1(A,FUN,FUN_FRECHET,FLAG) produces an estimate of
4 % the 1-norm relative condition number of function FUN at the matrix A.
5 % FUN and FUN_FRECHET are function handles:
6 % - FUN(A) evaluates the function at the matrix A.
7 % - If FLAG == 0 (default)
8 % FUN_FRECHET(B,E) evaluates the Frechet derivative at B
9 % in the direction E;
10 % if FLAG == 1
11 % - FUN_FRECHET('notransp',E) evaluates the
12 % Frechet derivative at A in the direction E.
13 % - FUN_FRECHET('transp',E) evaluates the
14 % Frechet derivative at A' in the direction E.
15 % If FUN_FRECHET is empty then the Frechet derivative is approximated
16 % by finite differences. More reliable results are obtained when
17 % FUN_FRECHET is supplied.
18 % MATLAB'S NORMEST1 (block 1-norm power method) is used, with a random
19 % starting matrix, so the approximation can be different each time.
20 % C = FUNM_CONDEST1(A,FUN,FUN_FRECHET,FLAG,P1,P2,...) passes extra inputs
21 % P1,P2,... to FUN and FUN_FRECHET.
22 % [C,EST] = FUNM_CONDEST1(A,...) also returns an estimate EST of the
23 % 1-norm of the Frechet derivative.
24 % Note: this function makes an assumption on the adjoint of the
25 % Frechet derivative that, for f having a power series expansion,
26 % is equivalent to the series having real coefficients.
27
28 if nargin < 3 || isempty(fun_frechet), fte_diff = 1; else fte_diff = 0; end
29 if nargin < 4 || isempty(flag1), flag1 = 0; end
30
31 n = length(A);
32 funA = feval(fun,A,varargin{:});
33 if fte_diff, d = sqrt( eps*norm(funA,1) ); end
34
35 factor = norm(A,1)/norm(funA,1);
36
37 [est,v,w,iter] = normest1(@afun);
38 c = est*factor;
39
40 %%%%%%%%%%%%%%%%%%%%%%%%% Nested function.
41 function Z = afun(flag,X)
42 %AFUN Function to evaluate matrix products needed by NORMEST1.
43
44 if isequal(flag,'dim')
45 Z = n^2;
46 elseif isequal(flag,'real')
47 Z = isreal(A);
48 else
49
50 [p,q] = size(X);
51 if p ~= n^2, error('Dimension mismatch'), end
52 E = zeros(n);
53 Z = zeros(n^2,q);
54 for j=1:q
55
56 E(:) = X(:,j);
57
58 if isequal(flag,'notransp')
59
60 if fte_diff
61 Y = (feval(fun,A+d*E/norm(E,1),varargin{:}) - funA)/d;
62 else
63 if flag1
64 Y = feval(fun_frechet,'notransp',E,varargin{:});
65 else
66 Y = feval(fun_frechet,A,E,varargin{:});
67 end
68 end
69
70 elseif isequal(flag,'transp')
71
72 if fte_diff
73 Y = (feval(fun,A'+d*E/norm(E,1),varargin{:}) - funA')/d;
74 else
75 if flag1
76 Y = feval(fun_frechet,'transp',E,varargin{:});
77 else
78 Y = feval(fun_frechet,A',E,varargin{:});
79 end
80 end
81
82 end
83
84 Z(:,j) = Y(:);
85 end
86
87 end
88
89 end
90
91 end