www.pudn.com > tfarma10.rar > invert_blocktri.m, change:2004-11-18,size:3047b

function Ainv= invert_blocktri(A, B, C) % function Ainv= invert_blocktri(A, B, C) % This file is part of the TFPM toolbox v1.0 (c) % michael.jachan@tuwien.ac.at and underlies the GPL. % % Huang/McColl, J. Phys. A: Math. Gen. 30 (1997) 7919--7933 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;tfpm; M= 10; N= 10; A= []; B= []; C= []; %for n= 1:N % A= [A -n*eye(M)]; % B= [B eye(M)]; % C= [C n*eye(M)]; %% C= [C n*zeros(M)]; %end; for n= 1:N A= [A randn(M)]; B= [B randn(M)]; C= [C randn(M)]; end; A(:, 1:M)= zeros(M); C(:, (N-1)*M+1:end)= zeros(M); X= zeros(M*N); for i= 1:N for j= 1:N if(i-1==j) X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= A(:, (i-1)*M+1:(i)*M); end; if(i==j) X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= B(:, (j-1)*M+1:(j)*M); end; if(j-1==i) X((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= C(:, (j-2)*M+1:(j-1)*M); end; end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M= size(A, 1); N= size(A, 2)/M; % Compute Z Z= [eye(M) inv(C(:, 1:M)) * B(:, 1:M)]; for i= 2:N-1 Z= [Z inv(C(:, (i-1)*M+1:(i)*M)) * ( B(:, (i-1)*M+1:(i)*M)*Z(:, (i-1)*M+1:(i)*M) - A(:, (i-1)*M+1:(i)*M)*Z(:, (i-2)*M+1:(i-1)*M) )]; end; % Compute Y Y= [inv(A(:, (N-1)*M+1:end)) * B(:, (N-1)*M+1:end) eye(M)]; for j= N-1:-1:2 Y= [inv(A(:, (j-1)*M+1:(j)*M)) * ( B(:, (j-1)*M+1:(j)*M)*Y(:, 1:M) - C(:, (j-1)*M+1:(j)*M)*Y(:, M+1:2*M) ) Y]; end; % Main Block Diag Ainv= zeros(M*N); Ainv(1:M, 1:M)= inv( B(:, 1:M) - C(:, 1:M)*Y(:, 1*M+1:2*M)*inv(Y(:, 1:M)) ); for j= 2:N-1 Ainv((j-1)*M+1:j*M, (j-1)*M+1:j*M)= inv( B(:, (j-1)*M+1:(j)*M) - A(:, (j-1)*M+1:(j)*M)*Z(:, (j-2)*M+1:(j-1)*M)*inv(Z(:, (j-1)*M+1:(j)*M)) - C(:, (j-1)*M+1:(j)*M)*Y(:, (j)*M+1:(j+1)*M)*inv(Y(:, (j-1)*M+1:(j)*M))); end; Ainv((N-1)*M+1:end, (N-1)*M+1:end)= inv( B(:, (N-1)*M+1:end) - A(:, (N-1)*M+1:end)*Z(:, (N-2)*M+1:(N-1)*M)*inv(Z(:, (N-1)*M+1:N*M)) ); % The Rest for i= 1:N for j= 1:N if(j<i) Ainv((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= ... -Y(:, (i-1)*M+1:(i)*M)*inv(Y(:, (i-2)*M+1:(i-1)*M))*Ainv((i-2)*M+1:(i-1)*M, (j-1)*M+1:(j)*M); end; end; end; for i= N:-1:1 for j= N:-1:1 if(i<j) Ainv((i-1)*M+1:(i)*M, (j-1)*M+1:(j)*M)= ... -Z(:, (i-1)*M+1:(i)*M)*inv(Z(:, (i )*M+1:(i+1)*M))*Ainv((i )*M+1:(i+1)*M, (j-1)*M+1:(j)*M); end; end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% norm(Ainv-inv(X)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%