matrix.e
Test.
--
--
--
--
--Matrix.e file by A. P. Adamson, Jan. 1997, for use with EUPHORIA language,
--based in part on the matrix Ada programs in matrix directory of
--Walnut Creek Ada CDRom and its predecessors.
--Send error notes to me at euclid at horandata.net
--Additional routines can also be added. If you have one you like, send it
--to me or post it to the news group.
--A modest amount of checking has been done but there may still be errors.
--Use at your own risk. I make no claim for the correctness of the routines.
--
--I define row vector as: sequence {1,2,3,...} and
--I define column vector as: sequence {{1},{2},{3},...}
--
--The following functions are implemented in the order shown.
--
--global function makeRowVector(integer iIn)
--global function makeColumnVector(integer iIn)
--global function make2DMatrix(integer rows, integer cols)
--global function make3DMatrix(integer rows, integer cols, integer layers)
--global function transpose( sequence matIn)
--global function dotRowRowVec (sequence seqIn)
--global function dotRowColVec (sequence seqIn)
--global function vectorsCrossProduct (sequence seqIn)
--global function matrixProduct (sequence seqIn)
--global function matrixTimesRowVector(sequence seqIn)
--global function matrixTimesColumnVector(sequence seqIn)
--global function rowVectorTimesMatrix(sequence seqIn)
--global function columnVectorTimesMatrix(sequence seqIn)
--global function matrixTrace(sequence matIn)
--global function abs(atom numIn)
--global function matrixPower(sequence seqIn)
--global function rot2DVectorCW90 (sequence seqIn)
--global function rot2DVectorCCW90 (sequence seqIn)
--global function reverseXin2DVector(sequence seqIn)
--global function reverseYin2DVector(sequence seqIn)
--global function getValueACrossB(sequence seqIn)
--global function getVectorLength(sequence seqIn)
--global function getCosIncludedAngleBetweenAB(sequence seqIn)
--global function crout(sequence seqIn, sequence vecIn)
--global function simultaneousInvert(sequence seqIn, sequence vecIn)
--
global function makeRowVector(integer iIn)
return repeat(0,iIn)
end function
--
-- ******************************************************************--
--
global function makeColumnVector(integer iIn)
return repeat(repeat(0,1),iIn)
end function
--
-- ******************************************************************--
--
global function make2DMatrix(integer rows, integer cols)
return repeat(repeat(0,cols), rows)
end function
--
-- ******************************************************************--
--
global function make3DMatrix(integer rows, integer cols, integer layers)
return repeat(repeat(repeat(0,cols), rows),layers)
end function
--
-- ******************************************************************--
--
global function transpose( sequence matIn)
-- This function performs the tranpose of a vector or 2 dimensional
-- matrix matIn.
-- (A) {1,2,3,4} becomes a column vector {{1},{2},{3},{4}}
-- (B) {{1},{2},{3},{4}} becomes a column vector {1,2,3,4}
-- (C) {{1,2,3,4},{-1,3,-1,0}} becomes {{1,-1},{2,3},{3,-2},{4,0}}
-- etc
--
sequence matOut
if atom(matIn[1]) then --case A
matOut = make2DMatrix(length(matIn),1)
for i = 1 to length(matIn) by 1 do
matOut[i][1] = matIn[i]
end for
return matOut
elsif
length(matIn[1]) = 1 then --case B
matOut = makeRowVector(length(matIn))
for i = 1 to length(matIn) by 1 do
matOut[i] = matIn[i][1]
end for
return matOut
else
matOut = make2DMatrix(length(matIn[1]),length(matIn)) --case C
for i = 1 to length(matIn) by 1 do
for j = 1 to length(matIn[1]) by 1 do
matOut[j][i] = matIn[i][j]
end for
end for
end if
return matOut
end function
--
-- ******************************************************************--
--
global function dotRowRowVec (sequence seqIn)
--seqIn is of the form {{a, b, c}, {d, e, f}} but more than 3
--terms cam be used
-- This function performs the inner (dot) product of two vectors A
-- and B resulting in a floating point number.
-- Comparability of dimensions is checked.
--
atom flo
flo = 0.0
if (length(seqIn[1]) = length(seqIn[2]) and
atom(seqIn[1][1]) and
atom(seqIn[2][1]))
then
for i = 1 to length(seqIn[1]) do
flo = flo + seqIn[1][i] * seqIn[2][i]
end for
return flo
else
puts (1, "vector lengths unequal\n")
return 0
end if
end function
--
-- ******************************************************************--
--
global function dotRowColVec (sequence seqIn)
--seqIn is of the form {{a, b, c}, {{d}, {e}, {f}}} but more than 3
--terms cam be used
--
-- This function performs the inner (dot) product of two vectors A,
-- a row vwctor and B, a column vector resulting in a floating point
-- number.
-- Comparability of dimensions is checked.
--
atom flo
flo = 0.0
if (length(seqIn[1]) = length(seqIn[2]) and
atom(seqIn[1][1]) and
length(seqIn[2][1] = 1))
then
for i = 1 to length(seqIn[1]) do
flo = flo + seqIn[1][i] * seqIn[2][i][1]
end for
return flo
else
puts (1, "vector lengths unequal\n")
return 0
end if
end function
--
-- *******************************************************************
--
global function vectorsCrossProduct (sequence seqIn)
-- This function performs the cross product of two 3 element vectors
-- resulting in a 3 element vector.
-- Usage, v = vectorsCrossProduct(seqIn) where
-- seqIn is of the form {{{number},{number},{number}},
-- {{number},{number},{number}}}
-- Compatability of dimensions is checked.
--
sequence sOut
sOut = makeRowVector(3)
if (length(seqIn[1]) = 3 and
length(seqIn[2]) = 3 and
atom(seqIn[1][1]) and
atom(seqIn[2][1]) )
then
sOut[1] = seqIn[1][2] * seqIn[2][3] - seqIn[1][3] * seqIn[2][2]
sOut[2] = seqIn[1][3] * seqIn[2][1] - seqIn[1][1] * seqIn[2][3]
sOut[3] = seqIn[1][1] * seqIn[2][2] - seqIn[1][2] * seqIn[2][1]
--
end if
return sOut
end function
-- ******************************************************************--
--
-- functions "+, -, *, /" to term by term add, subtract, multiply, divide
-- 2 vectors (or 2 matricees, or a number and a matrix or vector) is
-- inherent in Euphoria. C = A + B
--
-- ******************************************************************--
--
global function matrixProduct (sequence seqIn)
--seqIn is of the form seqIn[1] = mat1In with mat1In having m rows
--and n columns
-- seqIn[2] = mat2In with mat2In having n rows
--and q columns, ie, the number of columns in the first matrix equals
--the number of rows in the second matrix.
--m, n, q must be > 1.
--The output matrix will have m rows and q columns.
-- Comparability of dimensions is checked.
--
atom sum
sequence vecOut
sum = 0
if (length(seqIn[1][1]) != length(seqIn[2])) then
puts(1,"Error in dimensions of input matricees")
return 1=0
else
vecOut = make2DMatrix(length(seqIn[1]),length(seqIn[2][2]) )
for i = 1 to length(seqIn[1]) by 1 do
for j =1 to length(seqIn[2][1]) by 1 do
sum = 0
for k = 1 to length(seqIn[1][1]) by 1 do
sum = sum + (seqIn[1][i][k] * seqIn[2][k][j])
end for
vecOut[i][j] = sum
sum = 0
end for
end for
end if
return vecOut
end function
--
-- **********************************************************************
--
global function matrixTimesRowVector(sequence seqIn)
--
-- This function performs the multiplication of a matrix (seqIn[1])
-- and a row vector (seqIn[2] resulting in a row vector.
-- seqIn format is:
-- {{matIn},{vecIn}} where matIn is a 2 dimensional matrix and
-- vecIn is a 1 dimensional vector of
-- the same order and form {1,2,3...}. vecOut form is the same as vecIn.
-- Comparability of dimensions is checked.
-- **********************************************************************
--
sequence vecOut
atom sum
vecOut = makeRowVector(length(seqIn[2]))
if length(seqIn[1][1]) != length(seqIn[2])
then
return 1=0
else
for i = 1 to length(seqIn[1]) do
sum = 0.0
for k = 1 to length(seqIn[2]) do
sum = sum + seqIn[1][i][k] * seqIn[2][k]
vecOut[i] = sum
end for
end for
end if
return vecOut
end function
--
-- **********************************************************************
global function matrixTimesColumnVector(sequence seqIn)
--
-- **********************************************************************
-- This function performs the multiplication of a matrix (seqIn[1])
-- and a column vector (seqIn[2] resulting in a column vector.
-- seqIn format is:
-- {{matIn},{vecIn}} where matIn is a 2 dimensional matrix and
-- vecIn is a 1 dimensional vector of
-- the same order and form {{1},{2},{3}...}. vecOut form is the
-- same as vecIn.
-- Comparability of dimensions is checked.
-- **********************************************************************
--
sequence vecOut
atom sum
vecOut = makeRowVector(length(seqIn[2]))
if length(seqIn[1][1]) != length(seqIn[2])
then
return 1=0
else
for i = 1 to length(seqIn[1]) do
sum = 0.0
for k = 1 to length(seqIn[2]) do
sum = sum + seqIn[1][i][k] * seqIn[2][k][1]
vecOut[i] = sum
end for
end for
end if
return vecOut
end function
--
-- **********************************************************************
--
global function rowVectorTimesMatrix(sequence seqIn)
--
-- **********************************************************************
-- This function performs the multiplication of a row vector
-- (seqIn[1] and a matrix (seqIn[2]) resulting in a row vector.
-- seqIn format is: {{vecIn},{matIn}} where matIn is a 2 dimensional
-- matrix and vecIn is a 1 dimensional vector of form {1,2,3} and of
-- the same order. vecOut form is the same as vecIn.
-- Comparability of dimensions is checked.
-- **********************************************************************
--
sequence vecOut
atom sum
--
vecOut = makeRowVector(length(seqIn[1]))
--
if length(seqIn[1]) != length(seqIn[2][1])
then
return 1=0
else
for i = 1 to length(seqIn[2]) do
sum = 0.0
for k = 1 to length(seqIn[1]) do
sum = sum + seqIn[2][k][i] * seqIn[1][i]
vecOut[i] = sum
end for
end for
end if
return vecOut
end function
--***************************************************************************
--
global function columnVectorTimesMatrix(sequence seqIn)
--
-- **********************************************************************
-- This function performs the multiplication of a (column) vector
-- (seqIn[1] and a matrix (seqIn[2]) resulting in a (column) vector.
-- seqIn format is: {{vecIn},{matIn}} where matIn is a 2 dimensional
-- matrix and vecIn is a 1 dimensional vector of the form {{1},{2},{3}}
-- and of the same order. vecOut form is the same as vecIn.
-- Comparability of dimensions is checked.
-- **********************************************************************
--
sequence vecOut
atom sum
--
vecOut = makeColumnVector(length(seqIn[1]))
--
if length(seqIn[1]) != length(seqIn[2][1])
then
return 1=0
else
for i = 1 to length(seqIn[2]) do
sum = 0.0
for k = 1 to length(seqIn[1]) do
sum = sum + seqIn[2][k][i] * seqIn[1][i][1]
vecOut[i][1] = sum
end for
end for
end if
return vecOut
end function
--***************************************************************************
--
global function matrixTrace(sequence matIn)
--matIn is a 2 dimensional square matrix, format {{1,2,3},{},...}
--This function returns the "trace", ie the sum of the major diagonal terms
--
atom sum
sum = 0
for i = 1 to length(matIn[1]) by 1 do
sum = sum + matIn[i][i]
end for
return sum
end function
--
--************************************************************************
--
global function abs(atom numIn)
if numIn < 0 then return -numIn else return numIn end if
end function
--
--************************************************************************
--
global function matrixPower(sequence seqIn)
sequence A, B, C, L, M
integer rows, cols, P, I_PIVOT, J_PIVOT
atom BIG_ENTRY, TEMP, EPSILON
A = seqIn[1]
rows = length(seqIn[1])
cols = length(seqIn[1][1])
P = seqIn[2]
B = make2DMatrix(rows,rows)
C = make2DMatrix(rows,rows)
L = makeRowVector(rows)
M = makeRowVector(rows)
-- *******************************************************************
-- This function performs the square matrix operation of " matrix A
-- raise to integer power P ". When P is negative , say P = -N ,
-- matrixPower(A, -N) = (inverse(A))**N , that is, the inverse of A
-- raise to power N . In this case, matrix A must be non-singular.
-- Exceptions will be raised if the matrix A is not a square matrix,
-- or if matrix A is singular.
-- The form of seqIn is {square matrix, power} where seqIn[1] is
-- the input matrix and power is an integer, + or -. If power = -1,
-- the returned matrix is the inverse of the matrix in.
--
-- I tested the inversion operation for matrix size capability
-- on a Pentium 120, 32MHZ, DOS6.22 and found:
-- 80 X 80 OK with emm386 active OK, 100 NG
-- 240 X 240 OK with emm386 inactive OK,85 seconds for inversion.
-- Locked up for 300 X 300
--
-- Bob Craig tested on Pentium 150, WIN 95, 32MHZ and found:
-- 500 X 500 worked OK, 707 seconds for the inversion step
-- The size capability is no doubt due to config.sys details.
--
-- *******************************************************************
--
if rows != cols then
puts(1,"Input matrix is not square")
return 1=0
end if
--
if P = 0 then
-- B = identity matrix
for I = 1 to rows do
for J = 1 to rows do
if I != J then
B[I][J] = 0.0
else
B[I][J] = 1.0
end if
end for
end for
return B
end if
--
B = A
--
if P > 0 then
--& B = A multiplied itself for P times
for c = 1 to P-1 by 1 do
B = matrixProduct({A, B})
end for
return B
end if
--
-- P is negative, find inverse first
-- initiate the row and column interchange information
--
for K = 1 to rows do
L[K] = K -- row interchage information
M[K] = K -- column interchange information
end for
--
-- major loop for inverse
--
for K = 1 to rows do
-- & search for row and column index I_PIVOT, J_PIVOT
-- & both in (K .. B'LAST(1) ) for maximum B(I,J)
-- & in absolute value :BIG_ENTRY
BIG_ENTRY = 0.0
-- check matrix singularity
for I = K to rows do
for J = K to rows do
if abs(B[I][J]) > abs(BIG_ENTRY) then
BIG_ENTRY = B[I][J]
I_PIVOT = I
J_PIVOT = J
end if
end for
end for
--
if K = 1 then
if BIG_ENTRY = 0.0 then
puts(1, "SINGULAR Matrix") return 1=0
else
EPSILON = rows * (abs(BIG_ENTRY)) * 0.000001
end if
else
if abs(BIG_ENTRY) < EPSILON then
puts(1, "SINGULAR Matrix")
end if
end if
--
-- interchange row and column
--
--& interchange K-th and I_PIVOT-th rows
if I_PIVOT != K then
for J = 1 to rows do
TEMP = B[I_PIVOT][J]
B[I_PIVOT][J] = B[K][J]
B[K][J] = TEMP
end for
L[K] = I_PIVOT
end if
--& interchange K-th and J_PIVOT-th columns
if J_PIVOT !=K then
for I = 1 to rows do
TEMP = B[I][J_PIVOT]
B[I][J_PIVOT] = B[I][K]
B[I][K] = TEMP
end for
M[K] = J_PIVOT
end if
--
--& divide K-th column by minus pivot (-BIG_ENTRY)
--
for I = 1 to rows do
if I != K then
B[I][K] = B[I][K] / (-BIG_ENTRY)
end if
end for
--
-- reduce matrix row by row
--
for I = 1 to rows do
if I != K then
for J = 1 to rows do
if J != K then
B[I][J] = B[I][J] + B[I][K] * B[K][J]
end if
end for
end if
end for
--
--& divide K-th row by pivot
--
for J = 1 to rows do
if J != K then
B[K][J] = B[K][J] / BIG_ENTRY
end if
end for
B[K][K] = 1.0 / BIG_ENTRY
--
end for -- end of major inverse loop
--
-- final column and row interchange to obtain
-- inverse of A, i.e. A**(-1)
--
for K = rows to 1 by -1 do
-- column interchage
J_PIVOT = L[K]
if J_PIVOT != K then
--& intechange B(I,J_PIVOT) and B(I,K) for each row I
for I = 1 to rows do
TEMP = B[I][J_PIVOT]
B[I][J_PIVOT] = B[I][K]
B[I][K] = TEMP
end for
end if
-- row interchage
I_PIVOT = M[K]
if I_PIVOT != K then
--& INTECHANGE B(I_PIVOT,J) and B(K,J) for each column J
for J = 1 to rows do
TEMP = B[I_PIVOT][J]
B[I_PIVOT][J] = B[K][J]
B[K][J] = TEMP
end for
end if
end for
--
-- inverse of A is obtained and stored in B
-- now ready to handle the negative power
--
-- & C = B**(-P)
if P = -1 then
return B
end if
--
C = B
--
for I = P+1 to -1 by -1 do
--C = C * B
C = matrixProduct({C, B})
end for
return C
end function
--****************************************************************************
--
global function rot2DVectorCW90 (sequence seqIn)
--**********************************************************************
******
-- This function rotates a 2D vector 90 degrees cw and returns the
-- rotated vector. seqIn format is {2,4}
--**********************************************************************
******
return {seqIn[2],-seqIn[1]}
end function
--
--**********************************************************************
******
global function rot2DVectorCCW90 (sequence seqIn)
--**********************************************************************
******
-- This function rotates a 2D vector 90 degrees ccw and returns the
-- rotated vector. seqIn format is {2,4}
--**********************************************************************
******
return {-seqIn[2],seqIn[1]}
end function
--**********************************************************************
******
--
global function reverseXin2DVector(sequence seqIn)
--************************************************************************
-- This function rotates a 2D vector 180 degrees about the Y axis.
-- seqIn is of the form {2,4}
--************************************************************************
return {-seqIn[1], seqIn[2]}
end function
--****************************************************************************
--
global function reverseYin2DVector(sequence seqIn)
--**********************************************************************
******
-- Rotates 2D vector 180 degrees about the X axis. seqIn form is
-- {2,4}
--**********************************************************************
******
return {seqIn[1], -seqIn[2]}
end function
--************************************************************************
--
global function getValueACrossB(sequence seqIn)
--A : VEC2T; B : VEC2T) return FLOAT is
--************************************************************************
--Gets magnitude of A cross B for 2 2D vectors.
--************************************************************************
sequence A, B
A = seqIn[1]
B = seqIn[2]
return A[1] * B[2] - A[2] * B[1]
end function
--
-- ******************************************************************--
global function getVectorLength(sequence seqIn)
--seqIn is of the form {2,3,4}, assumed to be x, y, z coordinates of a vector
--in x,y,z space with x, y, z mutually perpendicular ()more than 3 coords
--can be used
atom a
a = 0
--return sqrt(seqIn[1] * seqIn[1] + seqIn[2] * seqIn[2] + seqIn[3] * seqIn[3])
for i = 1 to length(seqIn) by 1 do
a = a + seqIn[i] * seqIn[i]
end for
return sqrt(a)
end function
--
--************************************************************************
--
global function getCosIncludedAngleBetweenAB(sequence seqIn)
--seqIn is of form {vecA , vecB}
--Gets cos(THETA) where THETA is angle between 2 vectors.
--************************************************************************
sequence A, B
atom lenA, lenB
A = seqIn[1]
B = seqIn[2]
lenA = getVectorLength(A)
lenB = getVectorLength(B)
if lenA = 0 or lenB = 0 then puts(1,"One vec has 0 length!")
return 1=0
else
return ((dotRowRowVec({A,B})) / ( getVectorLength(A) *
getVectorLength(B)))
end if
end function
--
--************************************************************************
--
--following function solves a simultaneous equation set via Crout's
--reduction. It should take half the time of a solution via inversion.
--Actually it takes 4 times as long????
--include matrix.e
--
global function crout(sequence seqIn, sequence vecIn)
--
--For a n by n simultaneous equation system such as:
-- x1 + 2*x2 - x4 = 1
-- x1 + x2 - x3 + x4 = 4
-- 2*x1 - x2 + 2*x3 = 6
-- x1 + x2 - x4 = -1
--
--Seqin is the square matrix 1 2 0 -1
-- 1 1 -1 1
-- 2 -1 2 0
-- 1 1 0 -1
--Vecin is the column vector {{1},{4},{6},{-1}}
--
--They will be manipulated by the crout reduction and forward and backward
--substitution to give the solution to the equations, a new column vector,
--which is returned.
--Reference, Matricees for Engineers, by Alan D. Kraus, page 118
--
--The code fragment is:
sequence U, L, Y, X
integer mord
atom suml, sumu
mord = length(seqIn)
U = make2DMatrix(mord,mord) --upper triangle
L = make2DMatrix(mord,mord) --lower triangle
Y = makeRowVector(mord) --Y = U * X
X = makeRowVector(mord) --X = answer vector
--code for finding U and L follows
for k = 1 to mord do
U[k][k] = U[k][k] + 1.0
for i = k to mord do
if k = 1 then
L[i][k] = seqIn[i][k]
for j = k+1 to mord do
if k != j then
U[k][j] = (1.0/L[k][k]) * (seqIn[k][j])
end if
end for
else
suml = 0.0
for m = 1 to k-1 do
suml = suml + L[i][m] * U[m][k]
L[i][k] = seqIn[i][k] - suml
end for
end if
for j = k to mord do
sumu = 0.0
for m = 1 to k-1 do
sumu = sumu + L[k][m] * U[m][j]
end for
U[k][j] = (1.0/L[k][k]) * (seqIn[k][j] - sumu)
end for
end for
end for
--U and L are set by above
--
--Code Fragment for substitutions, set Y
for i = 1 to mord do --value Y per 5.20, page 125 for g
if i = 1 then
Y[i] = (1.0 / L[i][i]) * vecIn[i][1]
else
suml = 0.0
for m = 1 to i - 1 do
suml = suml + L[i][m] * Y[m]
end for
Y[i] = (1.0 / L[i][i]) * (vecIn[i][1] - suml)
end if
end for
--
--Now use back subs to set X
for k = mord to 1 by -1 do --value X per 5.21, page 125
sumu = 0.0
for m = k+1 to mord do
sumu = sumu + U[k][m] * X[m]
end for
X[k] = Y[k] - sumu
end for
return X
end function
--
--************************************************************************
--
global function simultaneousInvert(sequence seqIn, sequence vecIn)
--
--For a n by n simultaneous equation system such as:
-- x1 + 2*x2 - x4 = 1
-- x1 + x2 - x3 + x4 = 4
-- 2*x1 - x2 + 2*x3 = 6
-- x1 + x2 - x4 = -1
--
--seqIn is the square matrix 1 2 0 -1
-- 1 1 -1 1
-- 2 -1 2 0
-- 1 1 0 -1
--
--seqIn = {{1,2,0,-1,},{1,1,-1,1},{2,-1,2,0},{1,1,0,-1}}
--vecIn is the column vector {{1},{4},{6},{-1}}
--
--The solutions will be determined via inversion of A and multiplication
--by B with the result returned. It is a new column vector containing the
--solutions for x. Reference, Matricees for Engineers, by Alan D. Kraus
return matrixProduct({matrixPower({seqIn, -1}), vecIn})
end function
--
--************************************************************************
--end MATRIX_PACKAGE;
|
Not Categorized, Please Help
|
|
