uncomplicate / neanderthal

Fast Clojure Matrix Library
http://neanderthal.uncomplicate.org
Eclipse Public License 1.0
1.06k stars 56 forks source link

mm! doesn't seem to work with symmetric matrices #67

Closed kxygk closed 5 years ago

kxygk commented 5 years ago

Maybe I'm having some issue with my reading comprehension, or there is a bug

In the docs for mm! it says

If called with only 2 matrices, a and b, multiplies matrix a by a matrix b, and puts the result in the one that is a GE matrix. In this case, exactly one of a or b has to be, and one not be a GE matrix, but TR, or SY, or a matrix type that supports in-place multiplication.

I apologize in advance b/c I'm having a bit of trouble cooking up a simple example (b/c I don't know how to write a symmetric matrix by hand), but hopefully this illustrates the issue:

> (def A (dge [[43.0 36.0 38.0 90.0]
        [21.0 98.0 55.0 48.0]
        [72.0 13.0 98.0 12.0]
        [28.0 38.0 73.0 20.0]]))
#'linearsystems-part2.core/A
> (first-column-reflector A)
#RealUploMatrix[double, type:sy, mxn:4x4, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ↓       ┓    
   →       0.38    *       *       *            
   →       0.32    0.83    *       *            
   →       0.34   -0.17    0.82    *            
   →       0.80   -0.41   -0.44   -0.03         
   ┗                                       ┛    
> (mm! (first-column-reflector A)
        A)
Execution error (ExceptionInfo) at uncomplicate.neanderthal.internal.host.mkl.DoubleSYEngine/mm (mkl.clj:1692).
In-place mm! is not supported for SY matrices.

Basically I've got some code in (first-column-reflector ..) that spits out a symmetric matrix and yet I can't multiply in-place a GE and SY matrix like the docs say I should be able to

Or am I doing something wrong?

blueberry commented 5 years ago

I'm sorry, this is an error in the documentation. SY matrices support out of place mm (AB -> C) and NOT in-place mm (AB -> A or AB -> B). Various triangular matrices support in-place mm...

kxygk commented 5 years ago

Oh that's a shame. I will need to rework my algorithm. In that case, how would I properly do mm into a temporary matrix and then overwrite the data in my GE/input matrix? (to get the same end result)

blueberry commented 5 years ago

I don't have idea what your algorithm does, but I guess:

  1. AB -> C1
  2. C1D -> C2
  3. C2E -> C1 etc...
kxygk commented 5 years ago

I mean how in my example would I overwrite the data of the input matrix? with mm! I was hoping to do AB -> A, but now I need and intermediary matric C and to do AB -> C and then C -> A. I'm just not clear on what's the write way to do that

In what I'm writing the A is a sub matrix of a larger matrix which Im trying to update . It's an implementation of the Householder QR decomposition I'm trying to write out and the submarines are recursively being reduced to the upper triangular form.. it's all purely educational. I'm trying to learn to use your library and BLAS :)