# LAPACK exception with gels!

**URL:** https://discourse.julialang.org/t/lapack-exception-with-gels/43620
**Category:** Numerics
**Tags:** linearalgebra, lapack
**Created:** [July 24, 2020, 2:46pm UTC](https://discourse.julialang.org/t/lapack-exception-with-gels/43620 "2020-07-24T14:46:43Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)
#### Post date: [July 24, 2020, 2:46pm UTC](https://discourse.julialang.org/t/lapack-exception-with-gels/43620/1 "2020-07-24T14:46:43Z")

</div>

If I run the following program then Julia quits with an exception. I suspect this has something to do with rounding in the factorization. The matrix B below can be replaced with a matrix of zeros of the same size and the problem persists, but if A is replaced with a matrix of zeros then we get a solution. What is the best fix and what is the easiest fix for problems like this? Thanks!

To be clear: I understand that there is a discontinuity at zeros versus close to zero. But even if I set the remaining elements to exact zeros, I’m not getting something like a generalized inverse of A times B.

```julia
using LinearAlgebra
import LinearAlgebra.LAPACK.gels!

A = [0.06044523175910898 0.05417624043156454 -2.6133655153109947e-9 -2.6133655153109947e-9; 0.05417624043156454 0.1355075110293397 -6.536641412174804e-9 -6.536641412174804e-9; -2.6133655153109947e-9 -6.5366414139095275e-9 -4.0441353880288355e-8 -4.0441353880288355e-8; -2.6133655153109947e-9 -6.5366414139095275e-9 -4.0441353880288355e-8 -4.0441353880288355e-8]

B = [-0.08117075255109502 -0.04961506097021058 0.018683685925334124 0.018683688132780544 0.018683675632683132 0.01868368133744242 0.01868368337844292 0.018683710556390808 0.01868369486710419; -6.93889390390725e-18 -0.06474891160657195 0.009249843663413143 0.00924984475626679 0.00924983856776912 0.009249841392058273 0.009249842402508547 0.00924985585765707 0.009249848090268504; 1.1153959005309217e-25 5.84905020402288e-25 8.170386227774643e-9 8.170387188473825e-9 8.170381712421108e-9 8.170384213190477e-9 8.17038512532588e-9 8.170396991872857e-9 8.170390125098772e-9; 1.1479437019748901e-41 0.0 8.271806125530277e-25 8.271806125530277e-25 8.271806125530277e-25 8.271806125530277e-25 0.0 8.271806125530277e-25 8.271806125530277e-25]

gels!('N',A,B)

```

---

<div class="post-metadata">

### Author: ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)
#### Post date: [July 28, 2020, 7:43pm UTC](https://discourse.julialang.org/t/lapack-exception-with-gels/43620/2 "2020-07-28T19:43:28Z")

</div>

The matrix `A` is numerically singular

```julia
julia> cond(A)
2.2306564779458504e38

```

so you probably want to use a method that determines the rank as part of the solve and gives a solution to the rank deficient problem. In LAPACK that is `xgelsy` whereas `xgels` assumes that `A` has full rank. You can also get this with Julia’s API, i.e.

```julia
julia> ldiv!(qr!(copy(A), Val(true)), copy(B), eps())[1] ≈ gelsy!(copy(A), copy(B))[1]
true

```

where the `eps()` is the tolerance used in the determination of the rank.
