# HG changeset patch # User Antonio Pino Robles # Date 1432656894 -7200 # Node ID e0817ccb28345e66c03430b949898afbeea8a12a # Parent adcd5fc9b9a2421da5fc77a32483020e0be82820 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave. added funm_files/fun_atom.m added sqrtm2.m added toolbox/gecp.m toolbox/see.m: comment wrong call to subplot toolbox/tmtdemo.m: add a cast to double, as eig does not admit bool matrix input removed funm_files/funm_atom.m diff -r adcd5fc9b9a2 -r e0817ccb2834 funm_files/fun_atom.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/funm_files/fun_atom.m Tue May 26 18:14:54 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; diff -r adcd5fc9b9a2 -r e0817ccb2834 funm_files/funm_atom.m --- a/funm_files/funm_atom.m Sat May 16 11:39:36 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; diff -r adcd5fc9b9a2 -r e0817ccb2834 sqrtm2.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sqrtm2.m Tue May 26 18:14:54 2015 +0200 @@ -0,0 +1,121 @@ +function [X, alpha, condest] = sqrtm2(A) +%SQRTM2 Matrix square root. +% X = SQRTM2(A) is a square root of the matrix A (A = X*X). +% X is the unique square root for which every eigenvalue has +% nonnegative real part (the principal square root). +% If A is real with a negative real eigenvalue then a complex +% result will be produced. +% If A is singular then A may not have a square root. +% [X, ALPHA, CONDEST] = SQRTM2(A) returns a stability factor ALPHA and +% an estimate CONDEST for the matrix square root condition number of X. +% The residual NORM(A-X^2,'fro')/NORM(A,'fro') is bounded by +% (n+1)*ALPHA*EPS and the Frobenius norm relative error in X is +% bounded approximately by N*ALPHA*CONDEST*EPS, where N = MAX(SIZE(A)). + +% References: +% N. J. Higham, Computing real square roots of a real +% matrix, Linear Algebra and Appl., 88/89 (1987), pp. 405-430. +% A. Bjorck and S. Hammarling, A Schur method for the square root of a +% matrix, Linear Algebra and Appl., 52/53 (1983), pp. 127-140. + +n = max(size(A)); +[Q, T] = schur(A); % T is real/complex according to A. +[Q, T] = rsf2csf(Q, T); % T is now complex Schur form. + +% Compute upper triangular square root R of T, a column at a time. + +nzeig = length(find(diag(T)==0)); +if nzeig, warning('Matrix is singular and may not have a square root.'), end + +R = zeros(n); +for j=1:n + R(j,j) = sqrt(T(j,j)); + for i=j-1:-1:1 + s = 0; + for k=i+1:j-1 + s = s + R(i,k)*R(k,j); + end + if R(i,i) + R(j,j) ~= 0 + R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j)); + end + end +end + +if nargout > 1, alpha = norm(R,'fro')^2 / norm(T,'fro'); end + +X = Q*R*Q'; + +if nargout > 2 + + if nzeig + condest = inf; + else + + % Power method to get condition number estimate. + warns = warning; + warning('off'); + + tol = 1e-2; + x = ones(n^2,1); % Starting vector. + cnt = 1; + e = 1; + e0 = 0; + while abs(e-e0) > tol*e & cnt <= 6 + x = x/norm(x); + x0 = x; + e0 = e; + Sx = solve(R, x); + x = solve(R, Sx, 'T'); + e = sqrt(real(x0'*x)); % sqrt of Rayleigh quotient. + fprintf('cnt = %2.0f, e = %9.4e\n', cnt, e) + cnt = cnt+1; + end + + condest = e*norm(A,'fro')/norm(X,'fro'); + warning(warns); + + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function x = solve(R, b, tran) +%SOLVE Solves Kronecker system. +% x = SOLVE(R, b, TRAN) solves +% A*x = b if TRAN = '', +% A'*x = b if TRAN = 'T', +% where A = KRON(EYE,R) + KRON(TRANSPOSE(R),EYE). +% Default: TRAN = ''. + +if nargin < 3, tran = ''; end + +n = max(size(R)); +x = zeros(n^2,1); + +I = eye(n); + +if isempty(tran) + + % Forward substitution. + for i = 1:n + temp = b(n*(i-1)+1:n*i); + for j = 1:i-1 + temp = temp - R(j,i)*x(n*(j-1)+1:n*j); + end + x(n*(i-1)+1:n*i) = (R + R(i,i)*I) \ temp; + end + +elseif strcmp(tran,'T') + + % Back substitution. + for i = n:-1:1 + temp = b(n*(i-1)+1:n*i); + for j = i+1:n + temp = temp - conj(R(i,j))*x(n*(j-1)+1:n*j); + end + x(n*(i-1)+1:n*i) = (R' + conj(R(i,i))*I) \ temp; + end + +end + +return diff -r adcd5fc9b9a2 -r e0817ccb2834 toolbox/gecp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/gecp.m Tue May 26 18:14:54 2015 +0200 @@ -0,0 +1,56 @@ + +function [L, U, P, Q, rho] = gecp(A) +%GECP Gaussian elimination with complete pivoting. +% [L, U, P, Q, RHO] = GECP(A) computes the factorization P*A*Q = L*U, +% where L is unit lower triangular, U is upper triangular, +% and P and Q are permutation matrices. RHO is the growth factor. +% By itself, GECP(A) returns the final reduced matrix from the +% elimination containing both L and U. + +% from the Test Matrix Toolbox version 2.0 +% https://github.com/carandraug/testmatrix/blob/master/Test%20Matrix%20Toolbox/gecp.m + +[n, n] = size(A); +pp = 1:n; qq = 1:n; +if nargout == 5 + maxA = norm(A(:), inf); + rho = maxA; +end + +for k = 1:n-1 + + % Find largest element in remaining square submatrix. + % Note: when tie for max, no guarantee which element is chosen. + [colmaxima, rowindices] = max( abs(A(k:n, k:n)) ); + [biggest, colindex] = max(colmaxima); + row = rowindices(colindex)+k-1; col = colindex+k-1; + + % Permute largest element into pivot position. + A( [k, row], : ) = A( [row, k], : ); + A( :, [k, col] ) = A( :, [col, k] ); + pp( [k, row] ) = pp( [row, k] ); qq( [k, col] ) = qq( [col, k] ); + + if A(k,k) == 0 + break + end + + A(k+1:n,k) = A(k+1:n,k)/A(k,k); % Multipliers. + + % Elimination + i = k+1:n; + A(i,i) = A(i,i) - A(i,k) * A(k,i); + if nargout == 5, rho = max( rho, max(max(abs(A(i,i)))) ); end + +end + +if nargout <= 1 + L = A; + return +end + +L = tril(A,-1) + eye(n); +U = triu(A); + +if nargout >= 3, P = eye(n); P = P(pp,:); end +if nargout >= 4, Q = eye(n); Q = Q(:,qq); end +if nargout == 5, rho = rho/maxA; end diff -r adcd5fc9b9a2 -r e0817ccb2834 toolbox/see.m --- a/toolbox/see.m Sat May 16 11:39:36 2015 +0200 +++ b/toolbox/see.m Tue May 26 18:14:54 2015 +0200 @@ -59,6 +59,6 @@ end text(0,0,'Matrix not square.') end - subplot; +% subplot; end diff -r adcd5fc9b9a2 -r e0817ccb2834 toolbox/tmtdemo.m --- a/toolbox/tmtdemo.m Sat May 16 11:39:36 2015 +0200 +++ b/toolbox/tmtdemo.m Tue May 26 18:14:54 2015 +0200 @@ -217,9 +217,9 @@ c = []; e = []; j = 1; for i=1:matrix(0) A = full(matrix(i, 10)); - if norm(skewpart(A),1) % If not Hermitian... + if norm(skewpart(A),1) c1 = cond(A); - eg = eig(A); + eg = eig(double(A)); e1 = max(abs(eg)) / min(abs(eg)); % Filter out extremely ill-conditioned matrices. if c1 <= 1e10, c(j) = c1; e(j) = e1; j = j + 1; end