Hm,<br><br>My code snippets (tested on SBCL and CLISP+CYGWIN) work as expected.  Here is my matrix multiplication:<br>(matrix-product-triangular<br>           mat (matrix-product-triangular mat (copy x) 1<br>                          :Upper :NoTrans :NonUnit)<br>
           1 :Lower :NoTrans :Unit))<br><br>(I use (copy x) to create a working copy, so as not to obliterate x which I need later)<br><br><br><br><div class="gmail_quote">On Mon, Feb 23, 2009 at 5:18 PM, Liam Healy <span dir="ltr"><<a href="mailto:lhealy@common-lisp.net">lhealy@common-lisp.net</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Thanks for the example.  This actually reveals a bug.<br>
<br>
Taking a look at this example, I realize that LU functions needed to<br>
be simplified.  So, I have modified LU-solve and LU-decomposition in a<br>
way similar to what I did with cholesky recently; that is, I unified<br>
LU-solve and LU-solvex (solve in-place) into a single function, and<br>
now the final argument of each is optional, with the default being to<br>
make the correct vector for LU-decomposition and to solve in-place for<br>
LU-solve (if it's T, it will make the correct vector and put the<br>
solution there).  I think this is more convenient for users.<br>
<br>
All this works well.  I then decided to add a couple tests (one for<br>
double-float and one for complex-double-float) based on your example<br>
multiplying back the matrices with the solution.  Lo and behold, the<br>
answers come out reversed!<br>
<br>
Expected (#(42.730000000000004d0 -17.24d0 43.309999999999995d0))<br>
but saw (#(43.309999999999995d0 -17.24d0 42.730000000000004d0))<br>
<br>
I confirmed that this happens with both SBCL and CCL, so it's not the<br>
CL compiler that's a problem.  It also happens with the previous<br>
version; since it's not a regression, I've checked this in with a note<br>
in status.text that there's some bug here.  The only thing I know<br>
about it is that it's somewhere in the matrix-product-triangular<br>
calls, because I verified that it computes the right answer with<br>
matrix-product.<br>
<br>
Liam<br>
<div><div></div><div class="Wj3C7c"><br>
On Mon, Feb 23, 2009 at 9:12 AM, Mirko Vukovic <<a href="mailto:mirko.vukovic@gmail.com">mirko.vukovic@gmail.com</a>> wrote:<br>
> Hello,<br>
><br>
> Here is a snippet of LU decomposition code that illustrates the<br>
> decomposition, solution, and then solution check via blas based routines for<br>
> triangular matrix multiplication:<br>
><br>
> ;; decompose<br>
> (lu-decomposition mat per)<br>
> ;; solve for x given b<br>
> (lu-solve mat per b x)<br>
> ;; check by multiplying mat by x.  But mat is now a product of tridiagonal<br>
> matrices L and U.<br>
> ;; The LUx multiplication is done in two stages as L(Ux).  For the upper<br>
> multiplication we<br>
> ;; specify the matrix-product-triangular using the Upper and :NonUnit<br>
> diagonal.<br>
> ;; For the lower, we specify :Lower and :Unit diagonal<br>
> (matrix-product-triangular<br>
>           mat (matrix-product-triangular mat x 1<br>
>                          :Upper :NoTrans :NonUnit)<br>
>           1 :Lower :NoTrans :Unit))<br>
><br>
> One can off course save the matrix mat before the decomposition and then use<br>
> a direct matrix multiplication.<br>
><br>
> Mirko<br>
><br>
><br>
</div></div>> _______________________________________________<br>
> Gsll-devel mailing list<br>
> <a href="mailto:Gsll-devel@common-lisp.net">Gsll-devel@common-lisp.net</a><br>
> <a href="http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel" target="_blank">http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel</a><br>
><br>
><br>
</blockquote></div><br>