Thanks Holle and Mr. Gundla this is what have ended up with. Sharing it here because I think that people will find it useful:
The following is a hacked together LaTeX interface to a lua matrix computation module. The result is that matrix computations can be performed "in document". Matrixop.sty provides two macros: \matrixop<args> and \formatmatrix{matrix}. In one of its basic forms:
\matrixop{operation}{matrix1}{matrix2}
performs the computation and typesets the operation and result. For example the input on the left produces the output on the right:

Most operations provided by the matrix module are supported in both numeric and symbolic forms: addition, subtraction, multiplication, inversion, "division", determinants, vector scalar and cross products. For details put the files matrixop.sty, matrix.lua, matrixop.lua and matrixop-test.tex in the same directory (or install the .sty and put the lua files in the same directory) and compile matrixop-test.tex with lualatex.
File: matrixop.sty
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{matrixop}[2012/09/26 ver 1.0]
\RequirePackage{xparse}
\RequirePackage{fontspec}
\RequirePackage{mathtools}
\RequirePackage{geometry}
\DeclareMathOperator{\rref}{rref}
\directlua{dofile(kpse.find_file("matrixop.lua"))}
\directlua{dofile(kpse.find_file("matrix.lua"))}
\ExplSyntaxOn
\seq_new:N \g__matrop_symbops_seq
\seq_gset_split:Nnn \g__matrop_symbops_seq {,}{add,sub,mul,det,transpose,scalar,cross,add*,sub*,mul*,transpose*,cross*}
% Main Document Command: \matrixop
% - starred gives symbolic output
% - operation is first mandatory argument
% - operations are: add, sub, mul, div, invert, pow, root, det, dogauss, scalar and cross
% - starring the operation name, e.g \matrixop{mul*}... doesn't typeset operation but returns
% - the resulting matrix (not in useable form) in the macro \mout
% - root and pow require a numeric argument after the operation name giving the power/root, e.g. \matrixop{pow}[3]{matrix}
% - the final, optional, argument is the number of digits to round the output to
\NewDocumentCommand{\matrixop}{ s m o m g o }
{
\IfBooleanTF{#1}
{ % if star
\seq_if_in:NnTF \g__matrop_symbops_seq {#2}
{
\IfNoValueTF{#5}
{ % if unary
\directlua{matropsymb("#2","#4")}
}
{ % if binary
\directlua{matropsymb("#2","#4","#5")}
}
}
{
\msg_error:nnx {matrop} {Operation~not~compatible~with~symbolic} {#2}
}
}
{ % if no star
% if not root or pow
\IfNoValueTF{#3}
{ % if not binary
\IfNoValueTF{#5}
{ % if unary
% if not round
\IfNoValueTF{#6}
% then unary op, no round
{\directlua{matrop("#2","#4")}}
% then unary op, with round
{\directlua{matrop("#2","#4","#6")}}
}
{ % if binary
% if not round
\IfNoValueTF{#6}
% then binary, no round
{\directlua{matrop("#2","#4","#5")}}
% then binary, with round
{\directlua{matrop("#2","#4","#5","#6")}}
}
}
{ % root or pow
% if not round
\IfNoValueTF{#6}
% then root or pow, no round
{\directlua{matrop("#2","#3","#4")}}
% then root or pow, with round
{\directlua{matrop("#2","#3","#4","#6")}}
}
}
}
% This takes a macro that contains a matrix, and typesets the result.
% - first arg is optional to specify column justification
\NewDocumentCommand{\formatmatrix}{ O{r} m}
{
\matrixop_format:nn {#1}{#2}
}
% helper macro for the above
\cs_new:Npn \matrixop_format:nn #1#2
{
\tl_set:NV \l_tmpa_tl {#2}
\tl_replace_all:Nnn \l_tmpa_tl {)(}{\\}
\tl_replace_all:Nnn \l_tmpa_tl {,}{&}
\tl_set:Nx \l_tmpa_tl {\tl_tail:N \l_tmpa_tl}
\tl_replace_all:Nnn \l_tmpa_tl {)}{}
\begin{bmatrix*}[#1]
\tl_use:N \l_tmpa_tl
\end{bmatrix*}
}
\ExplSyntaxOff
File: matrixop.lua
function matropsymb (...)
-- desired operation is passed as first arg
local op = arg[1]
local matrix = require 'matrix'
local symbol = matrix.symbol
local m1, m2, mout
-- determine whether starred
local test = op
if string.match(test,"*") then
-- if so, then set flag and drop star
op = string.sub(op,1,-2)
star = true
else
star = false
end -- if string
-- if 2 args, must be unary
if #arg == 2 then
m1 = matrix{unpack(CreateMatrix(arg[2]))}:replace(symbol)
mout = matrix[op](m1)
-- if 3 args, must be binary
elseif #arg == 3 then
m1 = matrix{unpack(CreateMatrix(arg[2]))}:replace(symbol)
m2 = matrix{unpack(CreateMatrix(arg[3]))}:replace(symbol)
mout = matrix[op](m1,m2)--:replace(symbol)
end
-- if no star, then typeset operation and result
if not star then
if op == "mul" then
tex.sprint( matrix.latex(m1,"r").."\\cdot")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "add" then
tex.sprint( matrix.latex(m1,"r").."+")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "sub" then
tex.sprint( matrix.latex(m1,"r").."-")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "det" then
tex.sprint( "\\det"..matrix.latex(m1,"r").."=")
tex.sprint(mout)
elseif op == "transpose" then
tex.sprint(matrix.latex(m1,"r").."^{T}".."\\mkern -8mu=")
tex.sprint(matrix.latex(mout,"r"))
elseif op == "scalar" then
m1 = matrix.transpose(m1)
m2 = matrix.transpose(m2)
tex.sprint(matrix.latex(m1,"r").."\\cdot"..matrix.latex(m2,"r").."=")
tex.sprint(mout)
elseif op == "cross" then
tex.sprint(matrix.latex(m1,"r").."\\times"..matrix.latex(m2,"r").."=")
tex.sprint(matrix.latex(mout,"r"))
end -- if
end -- if not star
-- store output matrix in macro.
-- tex.print("\\gdef\\matout{"..ReturnMatrix(mout).."}")
if op ~= "det" then
if op ~="scalar" then
tex.print("\\gdef\\matout{"..ReturnMatrix(mout).."}")
end
end
end -- function
function matrop (...)
-- desired operation is passed as first arg
local op = arg[1]
local matrix = require 'matrix'
local m1, m2, m3, mout
-- determine whether starred
local test = op
if string.match(test,"*") then
op = string.sub(op,1,-2)
star = true
else
star = false
end
-- to test whether args are numbers or matrices
local test2 = tonumber(arg[2])
local test3 = tonumber(arg[3])
if #arg == 2 then
-- if 2 args, then unary
-- dogauss, alters original matrix so handled differently
m1 = matrix{unpack(CreateMatrix(arg[2]))}
if op == "dogauss" then
mout = matrix.copy(m1)
matrix[op](mout)
else
mout = matrix[op](m1)
end -- if gauss
-- several cases for 3 args
elseif #arg == 3 then
-- if second arg is a number-> root or pow
if test2 then
m1 = matrix{unpack(CreateMatrix(arg[3]))}
mout = matrix[op](m1,tonumber(arg[2]))
else
-- if not then if arg 3 is a number-> unary with round
-- as above, dogauss handled differently
if test3 then
m1 = matrix{unpack(CreateMatrix(arg[2]))}
if op == "dogauss" then
mout = matrix.copy(m1)
matrix[op](mout)
mout = matrix.round(mout,arg[3])
else
m2 = matrix[op](m1)
mout = matrix.round(m2,arg[3])
end -- if dogauss
-- if arg 3 is nan-> binary no round
else
m1 = matrix{unpack(CreateMatrix(arg[2]))}
m2 = matrix{unpack(CreateMatrix(arg[3]))}
mout = matrix[op](m1,m2)
end -- if test3
end -- if test2
else -- pow or root with round or binary with round
if test2 then
m1 = matrix{unpack(CreateMatrix(arg[3]))}
m2 = matrix[op](m1,tonumber(arg[2]))
mout = matrix.round(m2,arg[4])
else
m1 = matrix{unpack(CreateMatrix(arg[2]))}
m2 = matrix{unpack(CreateMatrix(arg[3]))}
m3 = matrix[op](m1,m2)
mout = matrix.round(m3,arg[4])
end -- if test2
end -- if #arg
-- no star-> typeset operation and result
if not star then
if op == "mul" then
tex.sprint( matrix.latex(m1,"r").."\\cdot")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "add" then
tex.sprint( matrix.latex(m1,"r").."+")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "div" then
tex.sprint( matrix.latex(m1,"r").."\\cdot")
tex.sprint( matrix.latex(m2,"r").."^{-1}".."\\mkern -8mu=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "sub" then
tex.sprint( matrix.latex(m1,"r").."-")
tex.sprint( matrix.latex(m2,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "invert" then
tex.sprint( matrix.latex(m1,"r").."^{-1}".."\\mkern -8mu=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "det" then
tex.sprint( "\\det"..matrix.latex(m1,"r").."=")
tex.sprint(mout)
elseif op == "dogauss" then
tex.sprint("\\rref"..matrix.latex(m1,"r").."=")
tex.sprint( matrix.latex(mout,"r"))
elseif op == "root" then
tex.sprint(matrix.latex(m1,"r").."^{\\frac{1}{"..arg[2].."}}".."\\mkern -8mu=")
tex.sprint(matrix.latex(mout,"r"))
elseif op == "transpose" then
tex.sprint(matrix.latex(m1,"r").."^{T}".."\\mkern -8mu=")
tex.sprint(matrix.latex(mout,"r"))
elseif op == "pow" then
tex.sprint(matrix.latex(m1,"r").."^{"..arg[2].."}".."\\mkern -8mu=")
tex.sprint(matrix.latex(mout,"r"))
elseif op == "scalar" then
m1 = matrix.transpose(m1)
m2 = matrix.transpose(m2)
tex.sprint(matrix.latex(m1,"r").."\\cdot"..matrix.latex(m2,"r").."=")
tex.sprint(mout)
elseif op == "cross" then
tex.sprint(matrix.latex(m1,"r").."\\times"..matrix.latex(m2,"r").."=")
tex.sprint(matrix.latex(mout,"r"))
end -- if op
end -- if not star
tex.print("\\gdef\\matout{"..ReturnMatrix(mout).."}")
end -- function
-- builds matrix table from input passed by tex
function CreateMatrix(str)
rows={}
-- get the elements between the braces
-- and execute the following function
string.gsub(str, "%((.-)%)",
function(sub)
-- for debugging
-- tex.print("Row: "..sub.."\\\\")
-- split the string at ','
elements = string.explode(sub,",")
row={}
for i,e in ipairs(elements) do
-- remove spaces (not really necessary)
e = string.gsub(e," ","")
-- insert the element to the row-table
table.insert(row,e)
-- for debugging
-- tex.print("Element :"..e.."\\\\")
end
-- insert the row-table to the rows-table
table.insert(rows,row)
end
)
return rows
end
-- returns output matrix in format that can be passed back
function ReturnMatrix(t)
if not (type(t)=="table") then
return t
else
local str = ""
for i=1,#t do
str=str.."("..t[i][1]
for j=2,#t[1] do
str=str..","..t[i][j]
end
str=str..")"
end
return str
end
end
File: matrix.lua
This file can be found here: https://github.com/davidm/lua-matrix
I modified the file matrix.lua slightly. In particular, I changed the matrix.latex function. For the macro to work you will need to edit that function to the following (lines 898-918) in the source:
--// matrix.latex ( mtx [, align] )
-- LaTeX output
function matrix.latex( mtx, align )
-- align : option to align the elements
-- c = center; l = left; r = right
-- \usepackage{dcolumn}; D{.}{,}{-1}; aligns number by . replaces it with ,
local align = align or "c"
local str = "\\begin{bmatrix*}[r]"
local getstr = matrix.type( mtx ) == "tensor" and tensor_tostring or number_tostring
for i = 1,#mtx do
str = str.."\t"..getstr(mtx[i][1])
for j = 2,#mtx[1] do
str = str.." & "..getstr(mtx[i][j])
end
-- close line
if i == #mtx then
str = str.."\n"
else
str = str.." \\\\\n"
end
end
return str.."\\end{bmatrix*}"
end
File: matrixop-test.tex
\documentclass{article}
\usepackage{matrixop}
\usepackage{showexpl}
\setlength{\parindent}{0pt}
\begin{document}
\begin{center}{\Huge Matrix operations in \TeX}\end{center}
\textbf{matrixop.sty} provides two commands: \verb=\matrixop= and \verb=\formatmatrix=. The package requires that the document be compiled with Lua\LaTeX. Matrices are passed to the macro either directly or in macros. A matrix is passed as bracket delimited rows with comma separated entries, e.g. $(1,2,3)(4,5,6)$ would represent a $2\times 3$ matrix. The \verb=\matrixop= macro itself takes a variety of forms.
\begin{enumerate}
\item \verb=\matrixop<args>= Performs the operation and typesets the operation and result.
\item \verb=\matrixop*<args>= Produces symbolic rather than calculated results.
\item \verb=\matrixop<args>[n]= Rounds the entries of the resulting matrix to $n$ decimal places.
\item \verb=\matrixop{op*}<args>= Performs the operation, but does not typeset the operation or result.
\end{enumerate}
In all cases, excluding symbolic determinant or symbolic dot product, the resulting matrix is stored in the macro \verb=\matout= in the same format as input. Thus, \verb=\matout= can be reused as input. Since \verb=\matout= is overwritten after each operation, it may be necessary to save it via, e.g. \verb=\let\mymat\matout=, so that it can be used later. The supported operations are
\begin{itemize}
\item Binary:
\begin{itemize}
\item add
\item sub
\item mul
\item div
\item scalar (dot product)
\item cross (cross product)
\end{itemize}
\item Unary:
\begin{itemize}
\item det
\item invert
\item dogauss (returns $\rref$)
\item root
\item pow
\item transpose
\end{itemize}
\end{itemize}
The operations that support symbolic output are: add, mul, sub, det, transpose, scalar, and cross.
\section{Usage}
\def\matA{(1,2,3)(4,5,6)(7,0,9)}
\def\matB{(1,2,-1)(2,2,4)(1,3,-3)}
\def\matC{(1,1,2)(2,3,4)(5,6,7)}
\def\vA{(1)(2)(3)}
\def\vB{(4)(5)(6)}
\def\vC{(a)(b)(c)}
\def\matD{(a,b,c)(d,e,f)(g,h,i)}
\def\matE{(r,s,t)(u,v,w)(x,y,z)}
\begin{LTXexample}[pos=r]
\def\matA{(1,2,3)(4,5,6)(7,0,9)}
\def\matB{(1,2,-1)(2,2,4)(1,3,-3)}
\[\matrixop{mul}{\matA}{\matB}\]
\end{LTXexample}
For symbolic output, we may use (symbolic can be used with numeric input as well)
\begin{LTXexample}[pos=b, wide=true]
\def\matD{(a,b,c)(d,e,f)(g,h,i)}
\def\matE{(r,s,t)(u,v,w)(x,y,z)}
\[\matrixop*{mul}{\matD}{\matE}\]
\end{LTXexample}
If we prefer that the output is not typeset, then a star can be appended to the operation name within the first argument. In this case the matrix (or number) resulting from the operation is stored in the \verb=\matout= macro and can be invoked as such before the next operation is executed, or stored for later use. With the following code:
\begin{LTXexample}[pos=r,rframe=]
\matrixop{invert*}{\matA} % output stored in \matout
\let\ainv\matout % save for later
\matrixop{invert*}{\matB} % etc.
\let\binv\matout
\matrixop{mul*}{\binv}{\ainv}[5]
\let\bainv\matout
\matrixop{mul*}{\matA}{\matB}
\let\abprod\matout
\matrixop{mul*}{\abprod}{\bainv}[3]
\end{LTXexample}
% this is just here because of grouping introduced by showexpl, ie. the definitions
% above don't leave the scope of the environment
\matrixop{invert*}{\matA}
\let\ainv\matout
\matrixop{invert*}{\matB}
\let\binv\matout
\matrixop{mul*}{\binv}{\ainv}[5]
\let\bainv\matout
\matrixop{mul*}{\matA}{\matB}
\let\abprod\matout
\matrixop{mul*}{\abprod}{\bainv}[3]
We could then type,
\begin{LTXexample}[pos=b]
If $A$ and $B$ are as given below,
\[
A=\formatmatrix{\matA}\qquad B=\formatmatrix{\matB}
\]
Then the product $AB$
\[
AB=\formatmatrix{\abprod}
\]
when multiplied (on the right) by the matrix $B^{-1}A^{-1}$
\[
B^{-1}A^{-1}=\formatmatrix{\bainv}
\]
should produce the identity matrix,
\[
AB(B^{-1}A^{-1})=\formatmatrix{\abprod}\cdot\formatmatrix{\bainv}=\formatmatrix{\matout}
\]
\end{LTXexample}
Some more examples below
\begin{LTXexample}[pos=b]
$\matrixop*{add}{\matD}{\matE}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{mul}{\matA}{\matB}$
\end{LTXexample}
\begin{LTXexample}[pos=b,wide=true]
$\matrixop*{mul}{\matD}{\matE}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{div}{\matA}{\matB}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{sub}{\matA}{\matB}$
\end{LTXexample}
\begin{LTXexample}[pos=b]
$\matrixop*{sub}{\matD}{\matE}$
\end{LTXexample}
\begin{LTXexample}[pos=b]
$\matrixop{invert}{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{invert}{\matA}[4]$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{det}{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=b]
$\matrixop*{det}{\matE}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{dogauss}{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=b]
$\matrixop{root}[5]{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{root}[5]{\matA}[4]$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{transpose}{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{pow}[3]{\matA}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{scalar}{\vA}{\vB}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop*{scalar}{\vA}{\vC}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop{cross}{\vA}{\vB}$
\end{LTXexample}
\begin{LTXexample}[pos=r]
$\matrixop*{cross}{\vA}{\vC}$
\end{LTXexample}
\end{document}