############################################################################
#
#	File:     matrix.icn
#
#	Subject:  Procedures for matrix manipulation
#
#	Authors:  Stephen B. Wampler and Ralph E. Griswold
#
#	Date:     December 2, 2000
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#
#  This file contains procedures for matrix manipulation.
#
############################################################################
#
#  Links: lu
#
############################################################################

link lu

procedure matrix_width(M)

   return *M[1]

end

procedure matrix_height(M)

   return *M

end

procedure write_matrix(file, M, x, s)
   local r, c, row, col
   /s := " "
   r := matrix_height(M)
   c := matrix_width(M)

   if /x then {				# packed, no punctuation
      every row := 1 to r do {
         every col := 1 to c do {
            writes(file, M[row][col], s)
            }
         write(file)
         }
      }
   else {
      every row := 1 to r do {
         writes(file, "[")
         every col := 1 to c do {
            writes(file, M[row][col], ", ")
            }
         write(file, "]")
         }
      }

end

procedure copy_matrix(M)
   local M1, n, i

   n := *M

   M1 := list(n)

   every i := 1 to n do
      M1[i] := copy(M[i])

   return M1

end

procedure create_matrix(n, m, x)
   local M

   M := list(n)
   every !M := list(m, x)

   return M

end

procedure identity_matrix(n, m)
   local r, c, M

   M := create_matrix(n, m, 0)

   every r := 1 to n do {
      every c := 1 to m do {
         if r = c then M[r][c] := 1
         }
      }

   return M

end

procedure add_matrix(M1, M2)
   local M3, r, c, n, m

   if ((n := *M1) ~= *M2) | ((m := *M1[1]) ~= *M2[1]) then
      stop("*** incorrect matrix sizes")

   M3 := create_matrix(n, m)

   every r := 1 to n do
      every c := 1 to m do
         M3[r][c] := M1[r][c] + M2[r][c]

   return M3

end

procedure mult_matrix(M1, M2)
   local M3, r, c, n, k

   if (n := *M1[1]) ~= *M2 then stop("*** incorrect matrix sizes")

   M3 := create_matrix(*M1,*M2[1])
   every r := 1 to *M1 do {
      every c := 1 to *M2[1] do {
         M3[r][c] := 0
         every k := 1 to n do {
             M3[r][c] +:= M1[r][k] * M2[k][c]
             }
         }
      }

   return M3

end

procedure invert_matrix(M)
   local M1, Y, I, d, i, n, B, j

   n := *M
   if n ~= *M[1] then stop("*** matrix not square")

   M1 := copy_matrix(M)
   Y := identity_matrix(n, n)
   I := list(n, 0)		# index vector

#  First perform LH decomposition on M1 (which changes it and produces
#  an index vector, I.

   d := lu_decomp(M1, I) | stop("*** singular matrix")

   every j := 1 to n do {
      B := list(n)		# work on columns
      every i := 1 to n do
         B[i] := Y[i][j]
      lu_back_sub(M1, I, B)	# does  not change M1 or I
      every i := 1 to n do	# put column in result
         Y[i][j] := B[i]
      }

   return Y

end

procedure determinant(M)
   local M1, I, result, i, n

   n := *M
   if n ~= *M[1] then stop("*** matrix not square")

   M1 := copy_matrix(M)
   I := list(n, 0)		# not used but required by lu_decomp()

   result := lu_decomp(M1, I) | stop("*** singular matrix")

   every i := 1 to n do		# determinant is product of diagonal
      result *:= M1[i][i]	# elements of the decomposed matrix

   return result

end