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 ZZ= [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 YY= [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 DiagAinv= 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 Restfor 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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%