Mercurial > matrix-functions
comparison mftoolbox/sqrtm_pulay.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 [X,k] = sqrtm_pulay(A,D) | |
| 2 %SQRTM_PULAY Matrix square root by Pulay iteration. | |
| 3 % [X,K] = SQRTM_PULAY(A,D) computes the principal square root of the | |
| 4 % matrix A using the Pulay iteration with diagonal matrix D | |
| 5 % (default: D = DIAG(DIAG(A))). D must have positive diagonal entries. | |
| 6 % K is the number of iterations. | |
| 7 % Note: this iteration is linearly converent and converges only when | |
| 8 % SQRTM(D) sufficiently well approximates SQRTM(A). | |
| 9 | |
| 10 if nargin < 2, D = diag(diag(A)); end | |
| 11 | |
| 12 if any(diag(D)<0) || ~isreal(D) | |
| 13 error('D must have positive, real diagonal.') | |
| 14 end | |
| 15 | |
| 16 n = length(A); | |
| 17 dhalf = sqrt(diag(D)); | |
| 18 Dhalf = diag(dhalf); | |
| 19 B = zeros(n); | |
| 20 maxit = 50; | |
| 21 | |
| 22 tol = mft_tolerance(A); | |
| 23 | |
| 24 for k = 1:maxit | |
| 25 | |
| 26 Bold = B; | |
| 27 B = (A - D - Bold^2) ./ (dhalf(:,ones(1,n)) + dhalf(:,ones(1,n))'); | |
| 28 X = Dhalf + B; | |
| 29 | |
| 30 reldiff = norm(B - Bold,inf)/norm(X,inf); | |
| 31 | |
| 32 if reldiff <= tol, return, end | |
| 33 | |
| 34 end | |
| 35 error('Not converged after %2.0f iterations', maxit) |