Thanks for the example. This actually reveals a bug.
Taking a look at this example, I realize that LU functions needed to be simplified. So, I have modified LU-solve and LU-decomposition in a way similar to what I did with cholesky recently; that is, I unified LU-solve and LU-solvex (solve in-place) into a single function, and now the final argument of each is optional, with the default being to make the correct vector for LU-decomposition and to solve in-place for LU-solve (if it's T, it will make the correct vector and put the solution there). I think this is more convenient for users.
All this works well. I then decided to add a couple tests (one for double-float and one for complex-double-float) based on your example multiplying back the matrices with the solution. Lo and behold, the answers come out reversed!
Expected (#(42.730000000000004d0 -17.24d0 43.309999999999995d0)) but saw (#(43.309999999999995d0 -17.24d0 42.730000000000004d0))
I confirmed that this happens with both SBCL and CCL, so it's not the CL compiler that's a problem. It also happens with the previous version; since it's not a regression, I've checked this in with a note in status.text that there's some bug here. The only thing I know about it is that it's somewhere in the matrix-product-triangular calls, because I verified that it computes the right answer with matrix-product.
Liam
On Mon, Feb 23, 2009 at 9:12 AM, Mirko Vukovic mirko.vukovic@gmail.com wrote:
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
Gsll-devel mailing list Gsll-devel@common-lisp.net http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel