jacobi.m

jacobi algorithm matlab code

%%%%%%%%%%%%%%%%%%%%%JACOBI ALGORITHM WITH SORTING%%%%%%%%%%%%%%%%%%%%%%%%%

%%INPUT
%%randMatrix is the matrix we wish to diagonalize. It's only random for the
% first interation
%%k is our current iteration
%%A is the matrix previous to randMatrix 

%%OUTPUT
%%D is the more diagonal matrix
%%bk is the off(D)
%%k is the number of iterations

function [D, bk, k, originalMatrix] = jacobi(randMatrix, k, A)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%Creating the random symmetric matrix [START]%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (size(randMatrix)==[0,0])
        
        %randMatrix is our random 5 x 5 symmetric matrix and is initialized
        %to be completely filled with 0's
        randMatrix = zeros(5,5);
        
        %If the row is equal to the column then we are on the diagonal
        %and thus the values needn't be mirrored. If we are not on
        %the diagonal, then we need to mirror the random value on the
        %other side of the diagonal to create a symmetric matrix. When 
        %a random value is stored in randMatrix(row, column), it is also 
        %stored in (column, row), thus "mirroring" the value. 
        %
        %We will loop through the x portion of the randMatrix:
        %      | x 0 0 0 0 |
        %      | x x 0 0 0 |
        %      | x x x 0 0 |
        %      | x x x x 0 |
        %      | x x x x x |
        %and mirror the off-diagonal elements. The diagonal elements will a
        %sense be "mirrored" too, but mirrored with its own position. 
        %Ex: row=2, column=2, randNumber is some random value
        %    then randMatrix (row, col) = randNumber;
        %         randMatrix (col, row) = randNumber;
        %    is   randMatrix (2,2) = randNumber;
        %         randMatrix (2,2) = randNumber;
        
        for (row = 1:5)
            for (col = 1:row)
                  %randNumber is a random number on the interval [0,10]  
                  randNumber = floor(rand*10);
                  %change the value at randMatrix(row, col) to the random number
                  randMatrix (row, col) = randNumber;
                  %mirror the value at randMatrix(col, row)
                  randMatrix (col, row) = randNumber;
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Creating the random symmetric matrix randMatrix [END]%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

originalMatrix = randMatrix;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Finding the max off-diagonal entry and the 2 x 2 matrix%%%%%%%%%% 
%%%%%%%%%%containing that entry [START]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %B will our 2 x 2 matrix
    B = zeros(2,2);
    %maxRow will be the row of the max entry
    maxRow = 1;
    %maxCol will be the col of the max entry
    maxCol = 2;

    for (row = 1:5)
        for (col = 1:5)
            %If the row is not equal to the column, then we are looping
            %throught the off-diagonal elements. 
            if (row ~= col)
                %%If the current entry is greater than previous max entry then
                %%set maxRow equal to the current row and the maxCol to the
                %%current column. 
                if (abs(randMatrix(row, col))> abs(randMatrix(maxRow, maxCol)))
                    maxRow = row;
                    maxCol = col;
                end
            end
        end
    end
    %%
    %So we've found the max entry. Now we need to use it to find a 2 x 2
    %symmetric matrix that we can diagonalize. The matrix basically follows
    %this formula which should hold true any n x n symmetric matrix where 
    %n >=3:
    %B = | randMatrix(maxRow, maxRow)  randMatrix(maxRow, maxCol) |
    %    | randMatrix(maxRow, maxCol)  randMatrix(maxCol, maxCol) | 
    B(1,2) = randMatrix(maxRow, maxCol);
    B(2,1) = randMatrix(maxRow, maxCol);
    B(1,1) = randMatrix(maxRow, maxRow);
    B(2,2) = randMatrix(maxCol, maxCol);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%Finding the max off-diagonal entry and the 2 x 2 matrix%%%%%%%%%% 
%%%%%%%%%%containing that entry [END]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%Diagonalizing 2 x 2 symmetric matrix B [START]%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % a, b, d will be used the calculation of the eigenvalues where:
    % B = |a b |
    %     |b d |
    a = B(1,1);
    b = B(1,2);
    d = B(2,2);
    I = [1 0; 0 1];
    
    %In order to diagonalize the matrix, we set the det(B-(mu*I)) equal to 
    %0 and solve for mu:
    %    0 = det | a - mu  b    |
    %            | b       d-mu |
    %   0 = (a-mu)*(d-mu) - b^2
    %   0 = mu^2 - (d+a)mu + (ad-b^2)
    %   In this form, we can use the quadatric formula to solve for mu:
    %       mu+ = myPlus  = ((d+a) + sqrt((d-a)^2 + 4*b^2))/2
    %       mu- = muMinus =((d+a) - sqrt((d-a)^2 + 4*b^2))/2
    muPlus = ((d+a) + sqrt((d-a).^2 + 4.*b.^2))./2;
    muMinus = ((d+a) - sqrt((d-a).^2 + 4.*b.^2))./2;
    
    %Now we plug in the eigenvalues where u1 and u2 are the results of that
    %plugging in, though it's not really necessary that we calculate the
    %u2
    u1 = B - (muPlus*I);
    u2 = B - (muMinus*I);
    
    %u1R1 is the first column of our soon to be normalized matrix u. u1R1
    %is obtained from the first row of the u1 with the appropriate sign
    %flipping. Ex:  u1 = |a b|
    %                    |c d|
    %               u1R1 = |-b|
    %                      | a|
    %but technically we could skip the skip of actually looking for the
    %eigenvalues is the sum of the first row is equal to the sum of the
    %second. Ex:    | 2 1 | = 3
    %               | 1 2 | = 3
    %We can assume that |1| to be our u1R1.
    %                   |1|
    u1R1 = [-u1(1,2) u1(1,1)]';
    if (sum(B(1,:))==sum(B(2,:)))
        u1R1 = [1 1]';
    end
    
    %u2R1 is the second column of matrix u and is obtained from u1R1 again
    %with the appropriate flipping. 
    %           Ex: u1R1 = |a|
    %                      |b|
    %               u2R1 = |-b|
    %                      | a|
    u2R1 = [-u1R1(2) u1R1(1)]';
  
    %u1Length is the 1/||u1R1||. The length for u2R1 is the same as u1R1 so
    %there is no need to calculate it again.
    u1Length = 1/(sqrt(u1R1(1)^2 + u1R1(2)^2));
 
    %The first column of matrix u is u1R1 and the second is u2R1. u is
    %normalized because its components are both normalized.
     u = u1Length*[u1R1 u2R1];
    
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%Creating the Givens Matrix [START]%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%G will be our Givens matrix, first filled with 0's 
        G = zeros(5,5);
        
        %We place entries from the u into G that correspond to the entries
        %that were taken from the original matrix randMatrix
        G(maxRow, maxCol) = u(1,2);
        G(maxCol, maxRow) = u(2,1);
        G(maxRow, maxRow) = u(1,1);
        G(maxCol, maxCol) = u(2,2);
       
        %We loop through G placing 1's on the diagonal where the value
        %is 0.
        for (row=1:5)
            for (col=1:5)
                if (row==col)
                    if (G(row,col)==0)
                        G(row,col)=1;
                    end
                end
            end
        end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%Creating the Givens Matrix [END]%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %D is our final more-diagonal matrix that we return
    D = G'*randMatrix*G;
    
    %bk is the value of the OFF(D), where OFF(x) is the sum of the square of 
    %each off-diagonal elements of x 
    bk = 0;
    for (row = 2:5)
            for (col = 1:row-1)
                bk = bk + 2*(D(row, col))^2;
            end
    end
    
    %increment k
    k=k+1;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%Diagonalizing 2 x 2 symmetric matrix B [END]%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%