1 files changed, 183 insertions, 0 deletions
diff --git a/ipl/procs/matrix.icn b/ipl/procs/matrix.icn
new file mode 100644
index 0000000..48007c4
--- /dev/null
+++ b/ipl/procs/matrix.icn
@@ -0,0 +1,183 @@
+############################################################################
+#
+#	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
+
+   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 produce of diagonal
+      result *:= M1[i][i]	# elements of the decomposed matrix
+
+   return result
+
+end