[Gsll-devel] lu solution followed by triangular matrix multiplication example

[Liam Healy]

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 at 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 at common-lisp.net
> http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel
>
>