Mercurial > matrix-functions
comparison mftoolbox/funm_condest1.m @ 0:8f23314345f4 draft
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 |