Mercurial > matrix-functions

annotate sqrtm2.m @ 8:a587712dcf5f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
funm_atom.m: rename fun_atom to funm_atom * funm_atom.m: rename fun_atom to funm_atom.
author Antonio Pino Robles <data.script93@gmail.com>
date Fri, 29 May 2015 09:48:36 +0200
parents e0817ccb2834
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
1 function [X, alpha, condest] = sqrtm2(A)
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
2 %SQRTM2 Matrix square root.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
3 % X = SQRTM2(A) is a square root of the matrix A (A = X*X).
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
4 % X is the unique square root for which every eigenvalue has
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
5 % nonnegative real part (the principal square root).
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
6 % If A is real with a negative real eigenvalue then a complex
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
7 % result will be produced.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
8 % If A is singular then A may not have a square root.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
9 % [X, ALPHA, CONDEST] = SQRTM2(A) returns a stability factor ALPHA and
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
10 % an estimate CONDEST for the matrix square root condition number of X.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
11 % The residual NORM(A-X^2,'fro')/NORM(A,'fro') is bounded by
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
12 % (n+1)*ALPHA*EPS and the Frobenius norm relative error in X is
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
13 % bounded approximately by N*ALPHA*CONDEST*EPS, where N = MAX(SIZE(A)).
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
14
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
15 % References:
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
16 % N. J. Higham, Computing real square roots of a real
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
17 % matrix, Linear Algebra and Appl., 88/89 (1987), pp. 405-430.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
18 % A. Bjorck and S. Hammarling, A Schur method for the square root of a
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
19 % matrix, Linear Algebra and Appl., 52/53 (1983), pp. 127-140.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
20
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
21 n = max(size(A));
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
22 [Q, T] = schur(A); % T is real/complex according to A.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
23 [Q, T] = rsf2csf(Q, T); % T is now complex Schur form.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
24
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
25 % Compute upper triangular square root R of T, a column at a time.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
26
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
27 nzeig = length(find(diag(T)==0));
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
28 if nzeig, warning('Matrix is singular and may not have a square root.'), end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
29
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
30 R = zeros(n);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
31 for j=1:n
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
32 R(j,j) = sqrt(T(j,j));
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
33 for i=j-1:-1:1
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
34 s = 0;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
35 for k=i+1:j-1
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
36 s = s + R(i,k)*R(k,j);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
37 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
38 if R(i,i) + R(j,j) ~= 0
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
39 R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j));
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
40 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
41 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
42 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
43
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
44 if nargout > 1, alpha = norm(R,'fro')^2 / norm(T,'fro'); end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
45
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
46 X = Q*R*Q';
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
47
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
48 if nargout > 2
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
49
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
50 if nzeig
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
51 condest = inf;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
52 else
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
53
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
54 % Power method to get condition number estimate.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
55 warns = warning;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
56 warning('off');
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
57
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
58 tol = 1e-2;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
59 x = ones(n^2,1); % Starting vector.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
60 cnt = 1;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
61 e = 1;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
62 e0 = 0;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
63 while abs(e-e0) > tol*e & cnt <= 6
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
64 x = x/norm(x);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
65 x0 = x;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
66 e0 = e;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
67 Sx = solve(R, x);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
68 x = solve(R, Sx, 'T');
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
69 e = sqrt(real(x0'*x)); % sqrt of Rayleigh quotient.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
70 fprintf('cnt = %2.0f, e = %9.4e\n', cnt, e)
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
71 cnt = cnt+1;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
72 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
73
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
74 condest = e*norm(A,'fro')/norm(X,'fro');
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
75 warning(warns);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
76
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
77 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
78
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
79 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
80
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
82 function x = solve(R, b, tran)
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
83 %SOLVE Solves Kronecker system.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
84 % x = SOLVE(R, b, TRAN) solves
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
85 % A*x = b if TRAN = '',
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
86 % A'*x = b if TRAN = 'T',
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
87 % where A = KRON(EYE,R) + KRON(TRANSPOSE(R),EYE).
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
88 % Default: TRAN = ''.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
89
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
90 if nargin < 3, tran = ''; end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
91
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
92 n = max(size(R));
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
93 x = zeros(n^2,1);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
94
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
95 I = eye(n);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
96
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
97 if isempty(tran)
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
98
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
99 % Forward substitution.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
100 for i = 1:n
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
101 temp = b(n*(i-1)+1:n*i);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
102 for j = 1:i-1
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
103 temp = temp - R(j,i)*x(n*(j-1)+1:n*j);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
104 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
105 x(n*(i-1)+1:n*i) = (R + R(i,i)*I) \ temp;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
106 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
107
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
108 elseif strcmp(tran,'T')
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
109
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
110 % Back substitution.
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
111 for i = n:-1:1
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
112 temp = b(n*(i-1)+1:n*i);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
113 for j = i+1:n
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
114 temp = temp - conj(R(i,j))*x(n*(j-1)+1:n*j);
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
115 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
116 x(n*(i-1)+1:n*i) = (R' + conj(R(i,i))*I) \ temp;
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
117 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
118
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
119 end
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
120
e0817ccb2834 Add square root matrix function file, rename atom at funm_files, modify old toolbox to run it inside GNU Octave.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
121 return