Hello,
Here is a snippet of LU decomposition code that illustrates the decomposition, solution, and then solution check via blas based routines for triangular matrix multiplication:
;; decompose
(lu-decomposition mat per)
;; solve for x given b
(lu-solve mat per b x)
;; check by multiplying mat by x. But mat is now a product of tridiagonal matrices L and U.
;; The LUx multiplication is done in two stages as L(Ux). For the upper multiplication we
;; specify the matrix-product-triangular using the Upper and :NonUnit diagonal.
;; For the lower, we specify :Lower and :Unit diagonal
(matrix-product-triangular
mat (matrix-product-triangular mat x 1
:Upper :NoTrans :NonUnit)
1 :Lower :NoTrans :Unit))
One can off course save the matrix mat before the decomposition and then use a direct matrix multiplication.
Mirko