Mercurial > matrix-functions
comparison mftoolbox/sqrtm_dbp.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,M,k] = sqrtm_dbp(A,scale) | |
| 2 %SQRTM_DBP Matrix square root by product form of Denman-Beavers iteration. | |
| 3 % [X,M,k] = SQRTM_DBP(A,SCAL) computes the principal square root X | |
| 4 % of the matrix A using the product form of the Denman-Beavers | |
| 5 % iteration. The matrix M tends to EYE. | |
| 6 % SCAL specifies the scaling: | |
| 7 % SCAL == 0, no scaling. | |
| 8 % SCAL == 1, determinantal scaling (default). | |
| 9 % k is the number of iterations. | |
| 10 | |
| 11 n = length(A); | |
| 12 if nargin < 2, scale = 1; end | |
| 13 | |
| 14 tol = mft_tolerance(A); | |
| 15 X = A; | |
| 16 M = A; | |
| 17 maxit = 25; | |
| 18 | |
| 19 for k = 1:maxit | |
| 20 | |
| 21 if scale == 1 | |
| 22 g = (abs(det(M)))^(-1/(2*n)); | |
| 23 X = g*X; M = g^2*M; | |
| 24 end | |
| 25 | |
| 26 Xold = X; invM = inv(M); | |
| 27 | |
| 28 X = X*(eye(n) + invM)/2; | |
| 29 M = 0.5*(eye(n) + (M + invM)/2); | |
| 30 | |
| 31 Mres = norm(M - eye(n),'fro'); | |
| 32 | |
| 33 reldiff = norm(X - Xold,'fro')/norm(X,'fro'); | |
| 34 if reldiff < 1e-2, scale = 0; end % Switch to no scaling. | |
| 35 | |
| 36 if Mres <= tol, return; end | |
| 37 | |
| 38 end | |
| 39 error('Not converged after %2.0f iterations', maxit) |