function AM= invert_akaike(rPM, BM)% function AM= invert_akaike(rPM, BM)%   This file is part of the TFPM toolbox v0.5 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Implements the Akaike algorithm (SIAMAM 24(2), mar73,% pp234) for solving the block-Toeplitz system %% btoep(rPM)*AP = BP; M...block size, P...super size,% size(rPM)= [M, M(2P-1)].%% Dimensions:PP= size(rPM);M= PP(1);P= (PP(2)/M+1)/2;% Permute rPM and transpose BM!!rrPM= [];for p= 1:2*P-1   rrPM= [rPM(:, (p-1)*M+1:p*M).' rrPM];end;rPM= rrPM;BM= BM.';rhs_dim= size(BM, 1);% Dimension of right-hand-side block-vector%% The output tensorAM= zeros(P*M, rhs_dim, P);%% Here is the algorithm for the system [A1..AP] T= [B1..BP]!!% Init (p= 0):qinv= rPM(:, (P-1)*M+1:P*M);sinv= qinv;am  = BM(:, 1:M)*inv(qinv);AM(:, :, 1)= [am.'; zeros(M*(P-1), rhs_dim)];if(P==1)   return;end;Xp= -rPM(:, (P-2)*M+1:(P-1)*M) * inv(qinv);%X1Yp= -rPM(:, (P  )*M+1:(P+1)*M) * inv(qinv);%Y1ap= rPM(:, (P-2)*M+1:(P-1)*M);             %a1rp= rPM(:, (P  )*M+1:(P+1)*M);             %r1qinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P-2)*M+1:(P-1)*M)*inv(qinv)*rPM(:, (P  )*M+1:(P+1)*M);sinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P  )*M+1:(P+1)*M)*inv(qinv)*rPM(:, (P-2)*M+1:(P-1)*M);for p= 1:P-2   Ep= block_exchange(M, p);   qinv= qinv_new;   sinv= sinv_new;   App= zeros(rhs_dim, M);   for m= 1:p      App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);   end;   App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);   am_new= [];   for m= 1:p      am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];   end;   am= [am_new App];   AM(:, :, p+1)= [am.'; zeros(M*(P-1-p), rhs_dim)];   Xp_new= [zeros(M) Xp] - (rPM(:, (P-2-p)*M+1:(P-1-p)*M) + Xp*   ap)*inv(sinv)*[eye(M) Yp];   Yp_new= [Yp zeros(M)] - (rPM(:, (P  +p)*M+1:(P+1+p)*M) + Yp*Ep*rp)*inv(qinv)*[Xp eye(M)];      ap= [ap; rPM(:, (P-2-p)*M+1:(P-1-p)*M)];   rp= [rp; rPM(:, (P  +p)*M+1:(P+1+p)*M)];      qinv_new= qinv - Xp_new(:, 1:M)          *Yp_new(:, p*M+1:(p+1)*M) * qinv;   sinv_new= sinv - Yp_new(:, p*M+1:(p+1)*M)*Xp_new(:, 1:M)           * sinv;   Xp= Xp_new;   Yp= Yp_new;end;p= P-1;qinv= qinv_new;sinv= sinv_new;App= zeros(rhs_dim, M);for m= 1:p   App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);end;App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);am_new= [];for m= 1:p   am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];end;am= [am_new App];AM(:, :, P)= am.';