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

[Mirko Vukovic]

Hm,

My code snippets (tested on SBCL and CLISP+CYGWIN) work as expected.  Here
is my matrix multiplication:
(matrix-product-triangular
           mat (matrix-product-triangular mat (copy x) 1
                          :Upper :NoTrans :NonUnit)
           1 :Lower :NoTrans :Unit))

(I use (copy x) to create a working copy, so as not to obliterate x which I
need later)



On Mon, Feb 23, 2009 at 5:18 PM, Liam Healy <lhealy at common-lisp.net> wrote:

> 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
> >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.common-lisp.net/pipermail/gsll-devel/attachments/20090223/950a25fe/attachment.html>