A basic test is included at the bottom which on touch computes a matrix inverse, product with the inverse (on left and right), and prints the results. Your help in testing for an eventual version 1.0 and publishing to the Scripting Library forum would be appreciated. Thank you and enjoy!
EDIT: This is the latest version. Old versions are archived below.
CODE
// "MHLSL 3D Matrix Library" by Hewee "MH" Zetkin
// License:
// The author hereby grants unlimited permission to copy, modify, and
// redistribute this script in any form, with the exception that any
// source code listing must retain this header giving credit to the author,
// similar licensing, and a basic summary of the changes made along with the
// contributors who made them.
// Description:
// Provides basic 3D (3x3, 3x1) matrix and matrix-vector operations in LSL,
// similar to the functionality found in BLAS levels 2 and 3 (with level 1
// provided by the LSL language itself) plus matrix transposition and
// inversion.
// Release Notes:
// v0.1.1 (Dec. 18, 2008)
// - List of row-vectors matrix encoding.
// - Basic matrix functions, matrix product, matrix-vector product
// - Matrix inverse
// - Matrix inverse test harness
// - Published to SL Scripting Tips forums in hopes of testing for
// version 1.0.
// v0.1.2 (Dec. 18, 2008)
// - Fixed matrix-vector product return type.
// "Close enough" floating point value for zero comparison and such.
float EPSILON = 0.001;
// The 3x3 identity matrix in row-major order, with each row represented as a
// vector.
list MATRIX_IDENTITY =
[ <1.0, 0.0, 0.0>,
<0.0, 1.0, 0.0>,
<0.0, 0.0, 1.0> ];
// Computes the sum of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix sum of 'm' and
// 'n'.
list matrixSum(list m, list n)
{
vector mr1 = llList2Vector(m, 0);
vector mr2 = llList2Vector(m, 1);
vector mr3 = llList2Vector(m, 2);
vector nr1 = llList2Vector(n, 0);
vector nr2 = llList2Vector(n, 1);
vector nr3 = llList2Vector(n, 2);
return [ mr1+nr1, mr2+nr2, mr3+nr3 ];
}
// Computes the difference of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix difference of 'm' and
// 'n'.
list matrixDiff(list m, list n)
{
vector mr1 = llList2Vector(m, 0);
vector mr2 = llList2Vector(m, 1);
vector mr3 = llList2Vector(m, 2);
vector nr1 = llList2Vector(n, 0);
vector nr2 = llList2Vector(n, 1);
vector nr3 = llList2Vector(n, 2);
return [ mr1-nr1, mr2-nr2, mr3-nr3 ];
}
// Computes the product of a 3x3 matrix and a scalar.
// See http://en.wikipedia.org/wiki/Matrix_(mathematics)#Basic_operations
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// param s - A scalar.
// returns - A matrix encoded like 'm' that is the product of 'm' and 's'.
list matrixScalarProduct(list m, float s)
{
vector r1 = llList2Vector(m, 0);
vector r2 = llList2Vector(m, 1);
vector r3 = llList2Vector(m, 2);
return [ s*r1, s*r2, s*r3 ];
}
// Transposes a 3x3 matrix.
// See http://en.wikipedia.org/wiki/Transpose
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the transpose of 'm'.
list matrixTranspose(list m)
{
vector r1 = llList2Vector(m, 0);
vector r2 = llList2Vector(m, 1);
vector r3 = llList2Vector(m, 2);
return [ <r1.x, r2.x, r3.x>,
<r1.y, r2.y, r3.y>,
<r1.z, r2.z, r3.z> ];
}
// Computes the product of two 3x3 matrix arguments.
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// param n - A 3x3 matrix encoded like 'm';
// returns - A matrix encoded like 'm' that is the matrix product of 'm' and
// 'n'.
list matrixProduct(list m, list n)
{
vector mr1 = llList2Vector(m, 0);
vector mr2 = llList2Vector(m, 1);
vector mr3 = llList2Vector(m, 2);
vector nr1 = llList2Vector(n, 0);
vector nr2 = llList2Vector(n, 1);
vector nr3 = llList2Vector(n, 2);
// columns
vector nc1 = <nr1.x, nr2.x, nr3.x>;
vector nc2 = <nr1.y, nr2.y, nr3.y>;
vector nc3 = <nr1.z, nr2.z, nr3.z>;
// product
float p11 = mr1*nc1; float p12 = mr1*nc2; float p13 = mr1*nc3;
float p21 = mr2*nc1; float p22 = mr2*nc2; float p23 = mr2*nc3;
float p31 = mr3*nc1; float p32 = mr3*nc2; float p33 = mr3*nc3;
return [ <p11, p12, p13>,
<p21, p22, p23>,
<p31, p32, p33> ];
}
// Computes the product of a 3x3 matrix and a 3x1 matrix (vector).
// See http://en.wikipedia.org/wiki/Matrix_product
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// param v - A 3x1 matrix encoded as a vector.
// returns - A vector that is the product of 'm' and 'v'.
vector matrixVectorProduct(list m, vector v)
{
vector mr1 = llList2Vector(m, 0);
vector mr2 = llList2Vector(m, 1);
vector mr3 = llList2Vector(m, 2);
return <mr1*v, mr2*v, mr3*v>;
}
// Computes the determinant of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Determinant
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - The determinant of 'm'.
float matrixDeterminant(list m)
{
vector r1 = llList2Vector(m, 0);
vector r2 = llList2Vector(m, 1);
vector r3 = llList2Vector(m, 2);
return r1*(r2%r3);
}
// Computes the cofactor matrix of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Cofactor_(linear_algebra)#Matrix_of_cofactors
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is composed of the cofactors of
// 'm'.
list matrixCofactors(list m)
{
vector r1 = llList2Vector(m, 0);
vector r2 = llList2Vector(m, 1);
vector r3 = llList2Vector(m, 2);
float m11 = r1.x; float m12 = r1.y; float m13 = r1.z;
float m21 = r2.x; float m22 = r2.y; float m23 = r2.z;
float m31 = r3.x; float m32 = r3.y; float m33 = r3.z;
// cofactors
float c11 = m22*m33-m32*m23;
float c12 = -m21*m33+m31*m23;
float c13 = m21*m32-m31*m22;
float c21 = -m12*m33+m32*m13;
float c22 = m11*m33-m31*m13;
float c23 = -m11*m32+m31*m12;
float c31 = m12*m23-m22*m13;
float c32 = -m11*m23+m21*m13;
float c33 = m11*m22-m21*m12;
return [ <c11, c12, c13>,
<c21, c22, c23>,
<c31, c32, c33> ];
}
// Computes the adjugate matrix of the 3x3 matrix argument.
// See http://en.wikipedia.org/wiki/Cofactor_(linear_algebra)#Adjugate
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the adjugate matrix to 'm'.
list matrixAdjugate(list m)
{
return matrixTranspose(matrixCofactors(m));
}
// Inverts a 3x3 matrix.
// See http://en.wikipedia.org/wiki/Matrix_inverse
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - A matrix encoded like 'm' that is the inverse of 'm'. If 'm' is
// not invertable, returns an empty list instead.
list matrixInverse(list m)
{
float det = matrixDeterminant(m);
if (llFabs(det) < EPSILON)
{
return [];
}
list adj = matrixAdjugate(m);
vector ar1 = llList2Vector(adj, 0);
vector ar2 = llList2Vector(adj, 1);
vector ar3 = llList2Vector(adj, 2);
return [ ar1/det, ar2/det, ar3/det ];
}
// Converts a matrix to a human-readable string for debugging output.
// param m - A 3x3 matrix in row-major order, each row represented as a vector
// in the list (so 'm11' is m[0].x and 'm23' is m[1].z). In other
// words:
// [ <m11, m12, m13>,
// [ <m21, m22, m23>,
// [ <m31, m32, m33> ]
// returns - A string formatted with rows separated by newlines and row
// elements separated by ", ".
string matrixToString(list m)
{
vector r1 = llList2Vector(m, 0);
vector r2 = llList2Vector(m, 1);
vector r3 = llList2Vector(m, 2);
float m11 = r1.x; float m12 = r1.y; float m13 = r1.z;
float m21 = r2.x; float m22 = r2.y; float m23 = r2.z;
float m31 = r3.x; float m32 = r3.y; float m33 = r3.z;
return (string)m11+", "+(string)m12+", "+(string)m13+"\n"+
(string)m21+", "+(string)m22+", "+(string)m23+"\n"+
(string)m31+", "+(string)m32+", "+(string)m33;
}
// Test Harness
list randomMatrix()
{
return [ <llFrand(1.0), llFrand(1.0), llFrand(1.0)>,
<llFrand(1.0), llFrand(1.0), llFrand(1.0)>,
<llFrand(1.0), llFrand(1.0), llFrand(1.0)> ];
}
default
{
touch_start(integer nDetected)
{
if (llDetectedKey(0) != llGetOwner())
{
return;
}
list m1 = randomMatrix();
list mi1 = matrixInverse(m1);
list p1 = matrixProduct(m1, mi1);
list rp1 = matrixProduct(mi1, m1);
llOwnerSay("Random matrix:\n"+matrixToString(m1));
llOwnerSay("Inverse:\n"+matrixToString(mi1));
llOwnerSay("Product:\n"+matrixToString(p1));
llOwnerSay("Reverse product:\n"+matrixToString(rp1));
}
}
