% csc148, Spring 1996
% Solution program for Assignment 3
% ---------------------------------------------------------------------------


% "factor.t"
%
% This file contains code for symbolically factoring a matrix as
% outlined in the Assignment 4 handout.
%
% A matrix is stored as an array of rows. Each row is represented
% by a list of row items. Each row item represents a non-zero value
% in the matrix.
%
% Note: Both a first and last pointer is kept for each row. This
% is because the procedure Cholesky always adds to the end of a row.
% Having a last pointer takes advantage of this fact.

type RowItem :
    record
        col : int
        next : ^RowItem
    end record

type Row :
    record
        first, last : ^RowItem
    end record



% The length on an output line.

const lineLength : int := 80



% setValue(M, row, col)
%     - Adds an item (row,col) to matrix M. Duplicates are NOT added.

procedure setValue (var M : array 1 .. * of Row, row, col : int)

    const n := upper (M)

    assert 1 <= row and row <= n and 1 <= col and col <= n

    % Add the row item at column 'col'.

    var RI : ^RowItem
    new RI
    assert RI not= nil
    RI -> col := col
    RI -> next := nil

    const first := M (row).first
    const last := M (row).last

    % Special cases:

    if first = nil then % First element ever.

        M (row).first := RI
        M (row).last := RI

    elsif col < first -> col then % New first element

        M (row).first := RI
        RI -> next := first

    elsif col > last -> col then % New last element

        last -> next := RI
        M (row).last := RI

    else % General case.

        % Look for insertion point

        var scan : ^RowItem := first
        loop
            exit when scan -> next -> col >= col
            scan := scan -> next
        end loop

        % Only insert if not already present.

        if scan -> next -> col not= col then
            RI -> next := scan -> next
            scan -> next := RI
        else
            free RI
        end if
    end if

end setValue



% obtainRow(M, row)
%     - This function returns a row of the matrix M.

function obtainRow (M : array 1 .. * of Row, row : int) : Row

    const n := upper (M)

    assert 1 <= row and row <= n

    result M (row)
end obtainRow



% obtainRowItem(R, col)
%     - This function returns the row item at column 'col' if it
% exists. Otherwise nil is returned.

function obtainRowItem (R : Row, col : int) : ^RowItem

    if R.first = nil or col > R.last -> col then
        result nil
    end if

    % Look for requested row item.

    % Special case of looking for last item.

    if R.last not= nil and R.last -> col = col then
        result R.last
    end if

    % General case of searching from the beginning of the list.

    var scan : ^RowItem := R.first
    loop
        exit when scan = nil or scan -> col > col
        if scan -> col = col then
            result scan
        end if
        scan := scan -> next
    end loop

    % Requested item wasn't in the row.

    result nil

end obtainRowItem



% value(M, i, j)
%     - This function returns the (boolean) value of M(i,j).

function value (M : array 1 .. * of Row, i, j : int) : boolean

    const n := upper (M)

    assert 1 <= i and i <= n and
        1 <= j and j <= n

    % Value is true iff it is present

    result obtainRowItem (obtainRow (M, i), j) not= nil

end value



% initMatrix(M)
%     - This procedure initializes a matrix to all zero. It is
% assumed that any allocated storage has already been freed.

procedure initMatrix (var M : array 1 .. * of Row)

    const n := upper (M)

    for row : 1 .. n
        M (row).first := nil
        M (row).last := nil
    end for

end initMatrix



% freeRowItems(R)
%     - This procedure frees any storage used by a list of row items.
% R is set to nil.

procedure freeRowItems (var R : ^RowItem)

    loop
        exit when R = nil

        var tmp : ^RowItem := R

        R := R -> next
        free tmp

    end loop

end freeRowItems



% freeMatrix(M)
%     - This procedure frees any storage used by a matrix. It
% is assumed that M has already been initialized.

procedure freeMatrix (var M : array 1 .. * of Row)

    const n := upper (M)

    for row : 1 .. n

        freeRowItems (M (row).first)
        M (row).last := nil

    end for

end freeMatrix



% copyMatrix(M1, M2)
%     - This procedures copies a matrix M1 into M2. This procedure
% initializes M2 before copying.

procedure copyMatrix (M1 : array 1 .. * of Row, var M2 : array 1 .. * of Row)

    const n := upper (M1)

    assert n = upper (M2)

    initMatrix (M2)

    for row : 1 .. n
        var scan : ^RowItem := obtainRow (M1, row).first
        loop
            exit when scan = nil
            setValue (M2, row, scan -> col)
            scan := scan -> next
        end loop
    end for

end copyMatrix



% permuteSymmetrixMatrix(M1, M2, p)
%     - This procedure permutes a matrix M1 according to a permutation 'p'
% and stores the result in M2. M2 is initialized before the permutation.
% M1 is assumed to be symmetric and only the lower left triangle of M2
% is actually stored.

procedure permuteSymmetricMatrix (M1 : array 1 .. * of Row,
        var M2 : array 1 .. * of Row,
        p : array 1 .. * of int)

    const n := upper (p)
    assert n = upper (M1) and n = upper (M2)

    initMatrix (M2)

    for row : 1 .. n
        const sourceRow := p (row)
        var scan : ^RowItem := obtainRow (M1, sourceRow).first
        loop
            exit when scan = nil
            const sourceColumn := scan -> col

            % Invert the permutation

            for col : 1 .. n
                if sourceColumn = p (col) then
                    if (col <= row) then
                        setValue (M2, row, col)
                    else
                        setValue (M2, col, row)
                    end if
                    exit
                end if
            end for

            scan := scan -> next
        end loop
    end for
end permuteSymmetricMatrix



% transposeMatrix(M1, M2)
%     - This procedure transposes a matrix M1 and stores the result in M2.
% M2 is initialized before the transposition occurs.

procedure transposeMatrix (M1 : array 1 .. * of Row,
        var M2 : array 1 .. * of Row)

    const n := upper (M1)

    assert n = upper (M2)

    initMatrix (M2)

    for row : 1 .. n
        var scan : ^RowItem := obtainRow (M1, row).first

        loop
            exit when scan = nil

            setValue (M2, scan -> col, row)

            scan := scan -> next
        end loop
    end for
end transposeMatrix



% removeUpperTriangle(M)
%     - This procedure frees the memory used by the upper triangular
% part of the matrix. It is assumed that the matrix has a non-zero
% diaganol.

procedure removeUpperTriangle (M : array 1 .. * of Row)

    const n := upper (M)

    for i : 1 .. n

        var diagElement : ^RowItem := obtainRowItem (obtainRow (M, i), i)

        assert diagElement not= nil

        freeRowItems (diagElement -> next)

    end for

end removeUpperTriangle



% putRow(M, row)
%     - This procedure outputs the sparsity pattern of row 'row' of M.
% No newlines are outputted.

procedure putRow (M : array 1 .. * of Row, row : int)

    const n := upper (M)

    assert 1 <= row and row <= n

    var numPrinted := 0

    var scan : ^RowItem := obtainRow (M, row).first
    loop
        exit when scan = nil
        put repeat (" ", scan -> col - numPrinted - 1) ..
        put "*" ..
        numPrinted := scan -> col

        scan := scan -> next
    end loop
    put repeat (" ", n - numPrinted) ..

end putRow



% putSymmetricRow(M, row)
%     - This procedure outputs the sparsity pattern of row 'row' of M.
% No newlines are outputted. M is assumed to be symmetric so only
% the lower left triangle is accessed.

procedure putSymmetricRow (M : array 1 .. * of Row, row : int)

    const n := upper (M)

    assert 1 <= row and row <= n

    var numPrinted := 0

    % First print columns <= row.

    var scan : ^RowItem := obtainRow (M, row).first
    loop
        exit when scan = nil or scan -> col > row

        put repeat (" ", scan -> col - numPrinted - 1) ..
        put "*" ..
        numPrinted := scan -> col

        scan := scan -> next
    end loop
    put repeat (" ", row - numPrinted) ..

    % Now print remaining columns

    for col : row + 1 .. n
        if value (M, col, row) then
            put "*" ..
        else
            put " " ..
        end if
    end for

end putSymmetricRow



% putMatrix(M)
%     - This procedure outputs the sparsity pattern of M.

procedure putMatrix (M : array 1 .. * of Row)

    const n := upper (M)

    for i : 1 .. n
        putRow (M, i)
        put ""
    end for

end putMatrix



% putOutputMatrices(M, L, Lt)
%     - This procedure outputs a matrix M, its factor L and L transpose
% as per the example in Figure 2 of the handout. M is assumed to be
% symmetric and only the lower left triangle is accessed.

procedure putOutputMatrices (M, L, Lt : array 1 .. * of Row)

    const indentAmount := 6
    const firstSeparationAmount := 8
    const secondSeparationAmount := 4

    const n := upper (M)

    assert n = upper (L) and n = upper (Lt)

    for i : 1 .. n
        put repeat (" ", indentAmount) ..
        putSymmetricRow (M, i)
        put repeat (" ", firstSeparationAmount) ..
        putRow (L, i)
        put repeat (" ", secondSeparationAmount) ..
        putRow (Lt, i)
        put ""
    end for

end putOutputMatrices



% putBanner(message)
%     - This procedure outputs a message centered on the line and
% surrounded by dashes. Two spaces are placed before and after the
% message. The output line is two less than maximum. This procedure
% uses the global constant 'lineLength'.

procedure putBanner (message : string)

    const messageLength := length (message)
    const dashLength := (lineLength - messageLength - 6) div 2

    put repeat ("-", dashLength) ..
    put "  ", message, "  " ..
    put repeat ("-", dashLength)

end putBanner



% putPermutation(P)
%     - This procedure outputs a permutation P.

procedure putPermutation (P : array 1 .. * of int)

    const n := upper (P)

    for i : 1 .. upper (P)
        put " ", P (i) ..
    end for
    put ""

end putPermutation



% doIntersect(rowA, rowB) : boolean
%     - This function returns true iff rowA and rowB both
% have items in some column. It is basically a merge (see
% mergesort, p. 304 in your notes).

function doIntersect (rowA, rowB : Row) : boolean

    var scanA : ^RowItem := rowA.first
    var scanB : ^RowItem := rowB.first
    loop
        exit when scanA = nil or scanB = nil

        if scanA -> col = scanB -> col then
            result true
        elsif scanA -> col < scanB -> col then
            scanA := scanA -> next
        else % scanA->col > scanB->col then
            scanB := scanB -> next
        end if
    end loop

    result false

end doIntersect



% symbolicCholeskyRow(M, L)
%     - This procedure symbolically factors a matrix M and
% stores the lower triangular factor in L. This version
% is row-oriented in that it processes a row at a time. The
% handout version of Cholesky processes a column at a time.
%
% The matrix L is initialized before the factoring begins.

procedure symbolicCholeskyRow (M : array 1 .. * of Row,
        var L : array 1 .. * of Row)

    const n := upper (M)

    assert n = upper (L)

    initMatrix (L)

    for i : 1 .. n
        var scan : ^RowItem := obtainRow (M, i).first

        for j : 1 .. i - 1
            if scan not= nil and scan -> col = j then

                setValue (L, i, j)
                scan := scan -> next

            elsif doIntersect (L (i), L (j)) then

                setValue (L, i, j)

            end if
        end for
        setValue (L, i, i)
    end for

end symbolicCholeskyRow



% Column-wise version.

procedure symbolicCholeskyCol (M : array 1 .. * of Row,
        var L : array 1 .. * of Row)

    const n := upper (M)

    assert n = upper (L)

    initMatrix (L)

    for j : 1 .. n
        for i : j + 1 .. n
            if value (M, i, j) or doIntersect (L (i), L (j)) then
                setValue (L, i, j)
            end if
        end for
        setValue (L, j, j)
    end for

end symbolicCholeskyCol



% getRow(M, row)
%     - This procedure inputs row 'row' of M according to the format
% specified by the handout. M is assumed to be initialized.

procedure getRow (var M : array 1 .. * of Row, row : int)

    const n := upper (M)

    assert 1 <= row and row <= n

    var line : string

    assert not eof

    get line : *

    assert length (line) = n

    for col : 1 .. n
        if line (col) = "*" then

            setValue (M, row, col)

        elsif line (col) not= " " then

            assert false % Can't happen

        end if
    end for
end getRow



% testMatrix(M, choleskyCandidate, error)
%     - This procedure tests a sparsity pattern to see if it is a
% valid candidate for factoring. See the handout for details. M is
% assumed to be initialized.

procedure testMatrix (M : array 1 .. * of Row,
        var choleskyCandidate : boolean,
        var reasonForRejection : string)

    const n := upper (M)

    for i : 1 .. n

        % Diaganol elements must be present

        if not value (M, i, i) then
            reasonForRejection := "Pattern has a zero on the diaganol."
            choleskyCandidate := false
            return
        end if

        % Matrix must be symmetric

        for j : 1 .. i - 1
            if value (M, i, j) not= value (M, j, i) then
                reasonForRejection := "Pattern is not symmetric."
                choleskyCandidate := false
                return
            end if
        end for
    end for

    reasonForRejection := ""
    choleskyCandidate := true

end testMatrix



% testSymbolicFactoring
%    - This procedure tests the symbolic factorization procedure
% and is based on the code-structure at the bottom of page 4 of
% the handout.


procedure testSymbolicFactoring

    var patternNumber : int := 0

    loop
        exit when eof

        patternNumber += 1

        var n, p : int
        var canFactorMatrix : boolean
        var reasonWhyNot : string

        get n
        var S : array 1 .. n of Row

        initMatrix (S)

        % Read in sparsity pattern

        for row : 1 .. n
            getRow (S, row)
        end for

        putBanner ("sparsity pattern "
            + intstr (patternNumber)
            + " start")

        % Test if pattern can be factored.

        testMatrix (S, canFactorMatrix, reasonWhyNot)

        if not canFactorMatrix then

            put "Pattern cannot be factored."
            put reasonWhyNot

        else

            % We don't need upper triangular part of the matrix
            % after testing symmetry.

            removeUpperTriangle (S)

        end if

        % Read in number of permutations.

        get p

        for : 1 .. p
            var permutation : array 1 .. n of int

            for i : 1 .. n
                get permutation (i)

                % Check validity of permutation

                assert 1 <= i and i <= n

                for j : 1 .. i - 1
                    assert permutation (i) not= permutation (j)
                end for
            end for

            if canFactorMatrix then
                var Scopy, L, Lt : array 1 .. n of Row

                permuteSymmetricMatrix (S, Scopy, permutation)

                symbolicCholeskyRow (Scopy, L)

                transposeMatrix (L, Lt)

                putPermutation (permutation)
                putOutputMatrices (Scopy, L, Lt)
                put ""

                % Free allocated memory

                freeMatrix (Scopy)
                freeMatrix (L)
                freeMatrix (Lt)
            end if
        end for

        freeMatrix (S)

        putBanner ("sparsity pattern "
            + intstr (patternNumber)
            + " end")
    end loop

end testSymbolicFactoring

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main program

testSymbolicFactoring