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
>
>