qr_method.m

MATLAB, computing matrix by QR

%%%%%%%%%%%%%%%%%%%%%%
%%   Main Function
function Matrix( )
    fprintf(' -------------------------------------------- \n');
    fid = fopen('matrix.txt', 'r');
    A = fscanf(fid, '%5f', [8 8]);
    B = fscanf(fid, '%5f', [1 8]);
    C = fscanf(fid, '%f',  [5 5]);
    A = C;
    %A = A';
    [Q, R] = GetQR( A );
    B = Q * R;
    fprintf('Compare matrix A with matrix B, which is equal to Q * R\n');
    check = Compare_Matrix(A, B);
    if( check == 1)
        fprintf('Matrix is equal to each other..\n');
        fprintf('Matrix Q is \n');
        Q
        fprintf('Matrix R is \n');
        R
        fprintf('\n');
    else
        fprintf('Matrix is not equal to each other..\n');
        fprintf('Something must be wrong!!!\n');
    end
    res = Singular_QR( A )
    res;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to Merge reflecting matrix to Unitary
function P = Merge_P( A, Q, v )
    dim_A = size(A);
    dim_v = size(v);
    if dim_A(1, 1) ~= dim_A(1, 2)
        H = 0;
        return
    end
    n_A = dim_A(1,1);
    n_v = dim_v(1,1);
    P   = eye(n_A);
    pos = n_A - n_v;
    
    for i = 1:n_v
        for j = 1:n_v
            P(pos+i, pos+j) = Q(i, j);
        end
    end
    
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to get the reflecting matrix
function [P, v] = SetReflecting( X, n ) %Fucntion is checked.
    NORM_X = norm(X);
    if X(1,1) == 0
        alpha  = NORM_X;
    else
        alpha = -sign(X(1,1))*NORM_X;
    end
    w = X;
    w(1, 1) = w(1, 1) - alpha;
    beta = 1/(NORM_X * (abs(X(1,1)) + NORM_X));
    P = eye(n) - beta * w * w';
    v = P*X;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to Merge reflecting matrix to Unitary
function P = Merge_P( A, Q, v )
    dim_A = size(A);
    dim_v = size(v);
    if dim_A(1, 1) ~= dim_A(1, 2)
        H = 0;
        return
    end
    n_A = dim_A(1,1);
    n_v = dim_v(1,1);
    P   = eye(n_A);
    pos = n_A - n_v;
    
    for i = 1:n_v
        for j = 1:n_v
            P(pos+i, pos+j) = Q(i, j);
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to Get the local Column of the matrix
function  v = Get_Column( A, col, beg )   
    dim = size(A);
    if dim(1, 1) ~= dim(1, 2)
        H = 0;
        return
    end
    N = dim(1,1);
    
    j = 1;
    for i = 1:N
        if i >= beg
            v(j, 1) = A(i,col);
            j = j + 1;
        end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to Build the Hessenberg Matrix 
%%  * PS: not only return the Hessenberg Matrix, but also return the Public
%%  Reflecting matrix.
% function  [H, pbc_P] = Build_Hess( A )
%     dim = size(A);
%     if dim(1, 1) ~= dim(1, 2)
%         H = 0;
%         return
%     end
%     N = dim(1,1);
%     pbc_P = eye(N);
%     
%     for i = 1:(N-2)
%         w = Get_Column(A, i, i+1)
%         dim = size(w);
%         dim_w = dim(1,1);
%         [P, v] = SetReflecting( w, dim_w )
%         P = Merge_P(A, P, v)
%         pbc_P = pbc_P * P
%         A = P * A * P
%     end
%     %H = A;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to Check two Matrix isequal
function check = Compare_Matrix(A, B)
    dim   = size(A);
    dim_A = dim(1,1);
    dim   = size(B);
    dim_B = dim(1,1);
    check = 0;
    TOL = 1e-10;
    
    if dim_A ~= dim_B
        fprintf('dim_A is not equal to dim_B...\n');
    else
        N = dim_A;
        for i = 1:N
            for j = 1:N
               if abs(A(i,j) - B(i, j)) > TOL
                    return
               end
           end
        end
    end
    check = 1;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to get the track of matrix
function sum = track( H )
    dim = size(H);
    N = dim(1,1);
    sum = 0;
    for i = 1:N
        sum = sum + H(i, i);
    end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Get the Power of the matrix
function res = pow_matrix( H, pow )
    dim = size(H);
    N = dim(1,1);
    res = eye(N);
    for i = 1:pow
        res = res * H;
    end
    
%%%%%%%%%%%%%%%%%
%%  Get coef
function P = get_coef( H )
    dim = size(H);
    N = dim(1,1);
    P(1:N) = 0;
    S(1:N) = 0;
          
    for i = 1:N
        M = pow_matrix(H, i);
        S(i) = track( M );
    end
    
    for i = 1:N
        P(i) = 0;
        sum = 0;
        for k = 1:(i-1)
            sum = sum - P(k)*S(i-k);
        end
        P(i) =( S(i) + sum )/i;
    end
    P = P * -1;
    
    Q(1:(N+1)) = 0;
    Q(1) = 1;
    for i = 2:(N+1)
        Q(i) = P(i-1);
    end
    P = Q;
    
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Set the polynormal
function   res = Singular_Polynormal( P, x )   
    dim = size(P);
    N = dim(1, 2);
    X(1:N) = 0;
    
    for i = 1:N
        X(i) = power(x, N-i);
    end
    res = P*X';
    
function   res = special(H, x)
    dim = size(H);
    N = dim(1,1);
    x
    for i = 1:N
        H(i,i) = H(i,i) - x;
    end
    res = det(H);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to get the Q, R matrix
function  [Q, R] = GetQR( A )
    dim = size(A);
    if dim(1, 1) ~= dim(1, 2)
        H = 0;
        return
    end
    N = dim(1,1);
    pbc_P = eye(N);
    
   for i = 1:(N)
        w = Get_Column(A, i, i);
        dim = size(w);
        dim_w = dim(1,1);
        [P, v] = SetReflecting( w, dim_w );
        P = Merge_P(A, P, v);
        % With this process, we can save 
        % the process to get the inverse of the matrix pbc_P
         pbc_P = pbc_P * P;
         A = P * A;
    end
    Q = pbc_P;
    R = A;
    
function res = Singular_QR( A )
    dim = size(A);
    if dim(1, 1) ~= dim(1, 2)
        H = 0;
        return
    end
    N = dim(1,1);
    
    Times = 100;
    for i = 1:Times
        [Q, R] = GetQR( A );
        A = R * Q
    end
    
    for i = 1:N
        res(i) = A(i,i);
    end