Add julia lu factorization. Fix some promotion for linear systems and fix bidiagonal and triangular solvers for multiple right hand side.

[andreasnoack]

I have added a julia implementation of the LU decomposition and adjusted the higher level functions to call that when eltype is not BlasType. The decomposition is stored in the same type as the LAPACK version. I don't think pivoting is so important when eltype is not BlasType since the result should be exact for rational matrices and for Matrix{BigFloat} you can always adjust the precision. If necessary, pivoting can be added later.

The promotion fixes should also fix #5301

This is progress towards fixing #3014

I have changed the tests to use multiple right hand side to catch the error in the Bidiagonal/Triangular solvers.

[andreasnoack]

Update:

• I have changed \ for bitarrays to return Array{Float32} which is consistent with true/true=1.0f0.
• I have also defined convert methods for SparseMatrixCSC without that don't specify the integer type. Should the integer type be dropped as parameter and just assumed to be BlasInt?
• Changed how eps is chosen in @test_approx_eps such that an eps is calculated for each argument and the largest is chosen.
• Changed Triangular and LU to be immutable. Actually, this was just to test if it worked, but I think it makes sense. You don't change the fields of these.

[andreasnoack]

The tests pass on both my Mac and Ubuntu.

[andreasnoack]

Maybe I should add a commercial for this pull request. With this you can do:

inv([big(1)//(i+j-1) for i = 1:5,j=1:5])
5x5 Array{Rational{BigInt},2}:
    25//1    -300//1     1050//1    -1400//1     630//1
  -300//1    4800//1   -18900//1    26880//1  -12600//1
  1050//1  -18900//1    79380//1  -117600//1   56700//1
 -1400//1   26880//1  -117600//1   179200//1  -88200//1
   630//1  -12600//1    56700//1   -88200//1   44100//1

[StefanKarpinski]

Ah, brilliant. Likely to be quite slow, unfortunately, but still, very cool.

[jiahao]

woot!!!!

[dmbates]

symmetrize! was written this way to follow the BLAS/LAPACK conventions. In BLAS/LAPACK only one triangle of a symmetric/Hermitian matrix is used or constructed. Thus if you use the BLAS routine *syrk to evaluate A'*A or A*A' only the upper (or lower, depending on the uplo argument) triangle is filled in. This is fine if you plan to use the result in further calls to BLAS/LAPACK routines only. However, if you want to use the result otherwise you must copy the strict upper triangle to the strict lower triangle and I called that operation symmetrize. If the name is confusing then another name should be created. But that particular operation is still needed as long as special matrix forms (Triangular, Symmetric/Hermitian) as used as general arrays.

[jiahao]

I think copytri! would be much less ambiguous than symmetrize!, since there are multiple conflicting ways to do the latter, e.g. here and here and here and here

[dmbates]

The copytri! name is fine with me.

[StefanKarpinski]

How about u2l in the name to indicate which way to goes. Otherwise I would never remember.

[andreasnoack]

Finally green.

What about the symmetrize! thing? I guess that could be done in a separate issue, which might also be more transparent than adding it to this request.

[jiahao]

Yes, let's discuss that separately.

[andreasnoack]

Okay, I think this one is ready to merge. I'll open a new issue for symmetrize!.

[ViralBShah]

Thus will probably not be much slower than other libraries in cases like BigFloat, ehere much of the time will be spent in arithmetic.

[andreasnoack]

I have generalized the LU factorization to work for rectangular matrices similarly to what LAPACK does. @jiahao, do you have any further comments before merging?

[jiahao]

I'm guess that there is a missing int-to-float conversion in the linear solver test which is causing the InexactError in the test suite.

[andreasnoack]

This is weird. I cannot find the problem and the tests pass on all machines I have access to. As I can see, the error is in the infimum norm calculation in line 587, but why?

[andreasnoack]

Could this have something to do with the new reducedim code?
cc: @timholy , @lindahua

[jiahao]

The computation of norm(M, Inf) now throws an InexactError for a matrix of Float16s.

julia> norm(convert(Matrix{Float16}, randn(5,5)), Inf)
ERROR: InexactError()
 in setindex! at array.jl:342

The test fails because of this.

[andreasnoack]

Noup. It is the Inf that is used. I have just found the problem. It is that abs returns a Matrix{Any} when the input matrix is Matrix{Float16}. I don't know why and I don't know why this doesn't happen on master.

[andreasnoack]

@jiahao Now I see what your comment was referring to. It was an intermediate error in my search for the real error.

I have devectorized norm(Matrix) which avoids some temporaries and at the same makes the tests pass. I have opened an issue for the problem that made the tests fail.

The tests are passing now, but the gcc machine timed out.

[jiahao]

I think we are good to go after the openlibm submodule hash is resolved.

[andreasnoack]

I would like to ask for some help on git here. I don't know why the change to the openlibm is included in the commit and I don't know how to get it out.

[jiahao]

🍰

[andreasnoack]

Thanks. That was fast.

[stevengj]

It is crazy not to do pivoting, even for BigFloat, because not doing so can absurdly increase the required precision.

[stevengj]

This should not have been merged without pivoting, in my opinion.

[stevengj]

(This is only really useful for BigFloat, quaternions, and other user-defined fp types, in my opinion. Rational systems overflow too quickly unless you use BigInt, but with BigInt the complexity is generally exponential due to growth of the denominators, and exponential-complexity algorithms are practically useless. So I don't see the point of merging this until it can be taken seriously for fp types.)

[jiahao]

I agree that we would want a pivoted version so that down the line we could have something that is at least somewhat useful for Float128. Anyway, we can work on this next.

[stevengj]

I would go further and say that if the pivoted version isn't implemented soon, this patch should be reverted; we shouldn't include broken functionality, and LU without partial pivoting is broken.