# [Numpy-discussion] Problem in LinearAlgebra?

Rob Hooft rob at hooft.net
Sun Nov 2 08:26:06 CST 2003

```Nadav Horesh wrote:
> Many unstable problems have a stable solution if you choose the right
> algorithm. The question is if somewhere the developers decided to switch
> the equations solver, or there is a real bug. I hope that one of the
> developers will reply in this forum. I will try to look at it also,
> since it is a core component in the linear algebra package. Also if you
> suspect that your work-around is useful --- please post it here!
>

The workaround is to use "generalized_inverse" instead of
"solve_linear_equations".

The changes in the latter routine since 17.1.2 are:

@@ -269,18 +408,37 @@
bstar = Numeric.zeros((ldb,n_rhs),t)
bstar[:b.shape[0],:n_rhs] = copy.copy(b)
a,bstar = _castCopyAndTranspose(t, a, bstar)
-    lwork = 8*min(n,m) + max([2*min(m,n),max(m,n),n_rhs])
s = Numeric.zeros((min(m,n),),real_t)
-    work = Numeric.zeros((lwork,), t)
+    nlvl = max( 0, int( math.log( float(min( m,n ))/2. ) ) + 1 )
+    iwork = Numeric.zeros((3*min(m,n)*nlvl+11*min(m,n),), 'l')
if _array_kind[t] == 1: # Complex routines take different arguments
-        lapack_routine = lapack_lite.zgelss
-        rwork = Numeric.zeros((5*min(m,n)-1,), real_t)
+        lapack_routine = lapack_lite.zgelsd
+        lwork = 1
+        rwork = Numeric.zeros((lwork,), real_t)
+        work = Numeric.zeros((lwork,),t)
results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
-                        0,work,lwork,rwork,0 )
+                        0,work,-1,rwork,iwork,0 )
+        lwork = int(abs(work[0]))
+        rwork = Numeric.zeros((lwork,),real_t)
+        a_real = Numeric.zeros((m,n),real_t)
+        bstar_real = Numeric.zeros((ldb,n_rhs,),real_t)
+        results = lapack_lite.dgelsd( m, n, n_rhs, a_real, m,
bstar_real,ldb , s, rcond,
+                        0,rwork,-1,iwork,0 )
+        lrwork = int(rwork[0])
+        work = Numeric.zeros((lwork,), t)
+        rwork = Numeric.zeros((lrwork,), real_t)
+        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
+                        0,work,lwork,rwork,iwork,0 )
else:
-        lapack_routine = lapack_lite.dgelss
+        lapack_routine = lapack_lite.dgelsd
+        lwork = 1
+        work = Numeric.zeros((lwork,), t)
+        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
+                        0,work,-1,iwork,0 )
+        lwork = int(work[0])
+        work = Numeric.zeros((lwork,), t)
results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
-                        0,work,lwork,0 )
+                        0,work,lwork,iwork,0 )
if results['info'] > 0:
raise LinAlgError, 'SVD did not converge in Linear Least Squares'
resids = Numeric.array([],t)

I'm not deep enough into this to know where the new version goes wrong.

Regards,

Rob Hooft

--
Rob W.W. Hooft  ||  rob at hooft.net  ||  http://www.hooft.net/people/rob/

```