Mercurial > matrix-functions
comparison mftoolbox/sqrtm_triang_min_norm.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 R = sqrtm_triang_min_norm(T) | |
| 2 %SQRTM_TRIANG_MIN_NORM Estimated min norm square root of triangular matrix. | |
| 3 % R = SQRTM_TRIANG_MIN_NORM(T) computes a primary square root of the | |
| 4 % upper triangular matrix T and attempts to minimize its 1-norm. | |
| 5 | |
| 6 if ~isequal(T,triu(T)), error('T must be upper triangular'), end | |
| 7 | |
| 8 n = length(T); | |
| 9 rp = zeros(n,1); | |
| 10 rm = zeros(n,1); | |
| 11 | |
| 12 R = zeros(n); | |
| 13 for j=1:n | |
| 14 rp(j) = sqrt(T(j,j)); | |
| 15 rm(j) = -sqrt(T(j,j)); | |
| 16 for i=j-1:-1:1 | |
| 17 rp(i) = (T(i,j) - R(i,i+1:j-1)*rp(i+1:j-1))/(R(i,i) + rp(j)); | |
| 18 rm(i) = (T(i,j) - R(i,i+1:j-1)*rm(i+1:j-1))/(R(i,i) + rm(j)); | |
| 19 end | |
| 20 if norm(rp(1:j),1) <= norm(rm(1:j),1) | |
| 21 R(1:j,j) = rp(1:j); | |
| 22 else | |
| 23 R(1:j,j) = rm(1:j); | |
| 24 end | |
| 25 end |