% Cyclic Block Penta-diagonal System Solver
% History: 20.jan.2008 MB Created
function X = CBPS( A, B, C, D, E, F, par)
    alp = 1;
    bet = -1;
    gam = 1;
    del = -1;
    if (nargin == 7) 
        if (length(par) == 4)
            alp = par(1);
            bet = par(2);
            gam = par(3);
            del = par(4);
        end
    end
% Local arrays    
    [m,m1,n]=size(A);
    T = zeros(m,m);
    S = eye(2*m,2*m);
    P = zeros(2*m,1);  
    B1(:,:) = B(:,:,1);
% Initialization    
    C(:,:,1) = C(:,:,1) - (bet/alp)*D(:,:,n);
    D(:,:,1) = D(:,:,1) - (bet/alp)*E(:,:,n);
    B(:,:,2) = B(:,:,2) - (del/gam)*E(:,:,n-1);    
    D(:,:,n-1) = D(:,:,n-1) - (gam/del)*A(:,:,2);    
    B(:,:,n) = B(:,:,n) - (alp/bet)*A(:,:,1);
    C(:,:,n) = C(:,:,n) - (alp/bet)*B(:,:,1);
% Factorization    
    C(:,:,1) = inv(C(:,:,1));
    D(:,:,1) = C(:,:,1)*D(:,:,1);
    E(:,:,1) = C(:,:,1)*E(:,:,1);
    C(:,:,2) = inv(C(:,:,2) - B(:,:,2)*D(:,:,1));
    D(:,:,2) = C(:,:,2)*(D(:,:,2) - B(:,:,2)*E(:,:,1));
    E(:,:,2) = C(:,:,2)*E(:,:,2);
    for k =3:n
        B(:,:,k) = B(:,:,k) - A(:,:,k)*D(:,:,k-2);
        C(:,:,k) = inv(C(:,:,k) - B(:,:,k)*D(:,:,k-1) - A(:,:,k)*E(:,:,k-2));
        if (k < n)
            D(:,:,k) = C(:,:,k)*(D(:,:,k)-B(:,:,k)*E(:,:,k-1));
        end
        if (k < n - 1)
            E(:,:,k) = C(:,:,k)*E(:,:,k);
        end
    end
% Intermediate solution
    F(:,1) = C(:,:,1)*F(:,1);
    C(:,:,1) = C(:,:,1)/alp;
    B(:,:,1) = 0;
    T = C(:,:,2);
    F(:,2)   = T*(F(:,2) - B(:,:,2)*F(:,1));
    C(:,:,2) = -T*(B(:,:,2)*C(:,:,1));
    B(:,:,2) = T/gam;
    for k = 3:n
        T = C(:,:,k);
        F(:,k) = T*(F(:,k)-B(:,:,k)*F(:,k-1)-A(:,:,k)*F(:,k-2));
        C(:,:,k) = -T*(B(:,:,k)*C(:,:,k-1)+A(:,:,k)*C(:,:,k-2));
        B(:,:,k) = -T*(B(:,:,k)*B(:,:,k-1)+A(:,:,k)*B(:,:,k-2));
        if (k == n)
            C(:,:,k) = C(:,:,k)+T/bet;
        end
        if (k == n-1)
            B(:,:,k) = B(:,:,k)+T/del;
        end
    end
% Back substitution
    F(:,n-1) = F(:,n-1)-D(:,:,n-1)*F(:,n);
    C(:,:,n-1) = C(:,:,n-1)-D(:,:,n-1)*C(:,:,n);
    B(:,:,n-1) = B(:,:,n-1)-D(:,:,n-1)*B(:,:,n);
    for k = n-2:-1:1
        F(:,k) = F(:,k) - D(:,:,k)*F(:,k+1) - E(:,:,k)*F(:,k+2);
        C(:,:,k) = C(:,:,k) - D(:,:,k)*C(:,:,k+1) - E(:,:,k)*C(:,:,k+2);
        B(:,:,k) = B(:,:,k) - D(:,:,k)*B(:,:,k+1) - E(:,:,k)*B(:,:,k+2);
    end
% Solution of auxiliary system 
    m1 = m+1;
    S(1:m,1:m) = S(1:m,1:m) + alp*(A(:,:,1)*C(:,:,n-1)+B1(:,:)*C(:,:,n))+bet*(D(:,:,n)*C(:,:,1)+E(:,:,n)*C(:,:,2));
    S(1:m,m1:2*m) = alp*(A(:,:,1)*B(:,:,n-1)+B1(:,:)*B(:,:,n))+bet*(D(:,:,n)*B(:,:,1)+E(:,:,n)*B(:,:,2));
    S(m1:2*m,1:m) = gam*A(:,:,2)*C(:,:,n) + del*E(:,:,n-1)*C(:,:,1);
    S(m1:2*m,m1:2*m) = S(m1:2*m,m1:2*m) + gam*A(:,:,2)*B(:,:,n) + del*E(:,:,n-1)*B(:,:,1);
    P(1:m) = alp*(A(:,:,1)*F(:,n-1)+B1(:,:)*F(:,n))+bet*(D(:,:,n)*F(:,1)+E(:,:,n)*F(:,2));
    P(m1:2*m) = gam*A(:,:,2)*F(:,n) + del*E(:,:,n-1)*F(:,1);
    P = S\P;
% Final solution
    for k = 1:n
        X(:,k) = F(:,k) - C(:,:,k)*P(1:m) - B(:,:,k)*P(m1:2*m);
    end
return
      
    
    