Added wrappers for LAPACK least-squares solvers.

There are more efficient least-squares solvers available in LAPACK
which are not used in Scipy. This commit includes only the wrappers
for LAPACK functions. In the next commit the modification of the
existing LAPACK least-squares solvers is done.
  modified:   scipy/linalg/flapack.pyf.src
  modified:   scipy/linalg/tests/test_lapack.py

scipy/linalg/flapack.pyf.src

@@ -763,6 +763,133 @@ interface
 
    end subroutine <prefix2>lasd4
 
+   subroutine <prefix2>gelsy(m,n,minmn,maxmn,nrhs,a,b,jptv,cond,r,work,lwork,info)
+
+   ! v,x,j,rank,work,info = dgelsy(a,b,jptv,cond,lwork,overwrite_a=True,overwrite_b=True)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,jptv,&cond,&r,work,&lwork,&info)
+     callprotoargument int*,int*,int*,<ctype2>*,int*,<ctype2>*,int*,int*,<ctype2>*,int*,<ctype2>*,int*,int*
+
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
+     <ftype2> dimension(m,n),intent(in,out,copy,out=v) :: a
+
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     <ftype2> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
+     intent(in,out,copy,out=x) b
+
+     <ftype2> intent(in) :: cond
+     integer intent(out,out=rank) :: r
+     integer intent(in,out,out=j),dimension(n),depend(n) :: jptv
+
+     ! LWORK is obtained by the query call
+     integer intent(in),depend(nrhs,m,n),check(lwork>=1||lwork==-1) :: lwork!=MAX(m*n+3*n+1,2*m*n+nrhs
+     <ftype2> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
+     integer intent(out)::info
+
+   end subroutine <prefix2>gelsy
+
+   subroutine <prefix2c>gelsy(m,n,minmn,maxmn,nrhs,a,b,jptv,cond,r,work,lwork,rwork,info)
+
+   ! v,x,j,rank,work,rwork,info = zgelsy(a,b,jptv,cond,lwork,overwrite_a=True,overwrite_b=True)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,jptv,&cond,&r,work,&lwork,rwork,&info)
+     callprotoargument int*,int*,int*,<ctype2c>*,int*,<ctype2c>*,int*,int*,<ctype2>*,int*,<ctype2c>*,int*,<ctype2>*,int*
+
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
+     <ftype2c> dimension(m,n),intent(in,out,copy,out=v) :: a
+
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     <ftype2c> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
+     intent(in,out,copy,out=x) b
+
+     <ftype2> intent(in) :: cond
+     integer intent(out,out=rank) :: r
+     integer intent(in,out,out=j),dimension(n),depend(n) :: jptv
+
+     ! LWORK is obtained by the query call
+     integer intent(in),depend(nrhs,m,n),check(lwork>=1||lwork==-1) :: lwork !=MAX(m*n+3*n+1,2*m*n+nrhs
+     <ftype2c> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
+     <ftype2> dimension(2*n),intent(out),depend(n) :: rwork
+     integer intent(out)::info
+
+   end subroutine <prefix2c>gelsy
+
+   subroutine <prefix2>gelsd(m,n,minmn,maxmn,nrhs,a,b,s,cond,r,work,lwork,size_iwork,iwork,info)
+
+   ! x,s,rank,work,iwork,info = dgelsd(a,b,lwork,size_iwork,cond=-1.0,overwrite_a=True,overwrite_b=True)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,s,&cond,&r,work,&lwork,iwork,&info)
+     callprotoargument int*,int*,int*,<ctype2>*,int*,<ctype2>*,int*,<ctype2>*,<ctype2>*,int*,<ctype2>*,int*,int*,int*
+
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
+     <ftype2> dimension(m,n),intent(in,copy) :: a
+
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     <ftype2> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
+     intent(in,out,copy,out=x) b
+
+     <ftype2> intent(in),optional :: cond=-1.0
+     integer intent(out,out=rank) :: r
+     <ftype2> intent(out),dimension(minmn),depend(minmn) :: s
+
+     integer intent(in),check(lwork>=1||lwork==-1) :: lwork
+     ! Impossible to calculate lwork explicitely, need to obtain it from query call first
+     ! Same for size_iwork
+     <ftype2> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
+
+     integer intent(in) :: size_iwork
+     integer intent(out),dimension(size_iwork),depend(size_iwork) :: iwork	
+     integer intent(out)::info
+
+   end subroutine <prefix2>gelsd
+
+   subroutine <prefix2c>gelsd(m,n,minmn,maxmn,nrhs,a,b,s,cond,r,work,lwork,size_rwork,rwork, size_iwork,iwork,info)
+
+   ! x,s,rank,work,rwork,iwork,info = zgelsd(a,b,lwork,size_rwork,size_iwork,cond=-1.0,overwrite_a=True,overwrite_b=True)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,s,&cond,&r,work,&lwork, rwork, iwork,&info)
+     callprotoargument int*,int*,int*,<ctype2c>*,int*,<ctype2c>*,int*,<ctype2>*,<ctype2>*,int*, <ctype2c>*,int*,<ctype2>*,int*,int*
+
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
+     <ftype2c> dimension(m,n),intent(in,copy) :: a
+
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     <ftype2c> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
+     intent(in,out,copy,out=x) b
+
+     <ftype2> intent(in),optional :: cond=-1.0
+     integer intent(out,out=rank) :: r
+     <ftype2> intent(out),dimension(minmn),depend(minmn) :: s
+
+     integer intent(in),check(lwork>=1||lwork==-1) :: lwork
+     ! Impossible to calculate lwork explicitely, need to obtain it from query call first
+     ! Same for size_rwork, size_iwork
+     <ftype2c> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
+
+     integer intent(in) :: size_rwork
+     <ftype2> intent(out),dimension(size_rwork),depend(size_rwork) :: rwork
+
+     integer intent(in) :: size_iwork
+     integer intent(out),dimension(size_iwork),depend(size_iwork) :: iwork	
+     integer intent(out)::info
+
+   end subroutine <prefix2c>gelsd
 
    subroutine <prefix2>geqp3(m,n,a,jpvt,tau,work,lwork,info)
 

scipy/linalg/tests/test_lapack.py

@@ -131,6 +131,106 @@ def test_clapack(self):
             # clapack module is empty
             pass
 
+class TestLeastSquaresSolvers(TestCase):
+
+    def test_gelsd(self):
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gelsd, = get_lapack_funcs(('gelsd',), (a1, b1))
+            #x,s,rank,work,iwork,info = *gelsd(a,b,lwork,size_iwork,[cond,overwrite_a,overwrite_b])
+            x, s, rank, work, iwork, info = gelsd(a1, b1, -1, 1)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+            iwork_size = iwork[0].real.astype(np.int)
+
+            x, s, rank, work, iwork, info = gelsd(a1, b1, lwork, iwork_size, -1, False, False)
+            assert_allclose(x[:-1], np.array([-14.333333333333323, 14.999999999999991], dtype=dtype), rtol=10*np.finfo(dtype).eps)
+            assert_allclose(s, np.array([12.596017180511966, 0.583396253199685], dtype=dtype), rtol=10*np.finfo(dtype).eps)
+
+        for dtype in COMPLEX_DTYPES:
+            a1 = np.array([[1.0+4.0j,2.0],
+                          [4.0+0.5j,5.0-3.0j],
+                          [7.0-2.0j,8.0+0.7j]], dtype=dtype)
+            b1 = np.array([16.0, 17.0+2.0j, 20.0-4.0j], dtype=dtype)
+            gelsd, = get_lapack_funcs(('gelsd',), (a1, b1))
+            #x,s,rank,work,rwork,iwork,info = cgelsd(a,b,lwork,size_rwork,size_iwork,[cond,overwrite_a,overwrite_b])
+            x, s, rank, work, rwork, iwork, info = gelsd(a1, b1, -1, 1,1)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+            rwork_size = rwork[0].real.astype(np.int)
+            iwork_size = iwork[0].real.astype(np.int)
+
+            x, s, rank, work, rwork, iwork, info = gelsd(a1, b1, lwork, rwork_size, iwork_size, -1, False, False)
+            assert_allclose(x[:-1], np.array([1.161753632288328-1.901075709391912j, 1.735882340522193+1.521240901196909j],
+                            dtype=dtype), rtol=10*np.finfo(dtype).eps)
+            assert_allclose(s, np.array([13.035514762572043, 4.337666985231382], dtype=dtype), rtol=10*np.finfo(dtype).eps)
+
+    def test_gelss(self):
+
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gelss, = get_lapack_funcs(('gelss',), (a1, b1))
+            #v,x,s,rank,work,info = dgelss(a,b,[cond,lwork,overwrite_a,overwrite_b])
+            v,x,s,rank,work,info = gelss(a1, b1,-1, -1, False, False)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+
+            v,x,s,rank,work,info = gelss(a1, b1,-1,lwork, False, False)
+            assert_allclose(x[:-1], np.array([-14.333333333333323, 14.999999999999991], dtype=dtype),rtol=10*np.finfo(dtype).eps)
+            assert_allclose(s, np.array([12.596017180511966, 0.583396253199685], dtype=dtype), rtol=10*np.finfo(dtype).eps)
+
+        for dtype in COMPLEX_DTYPES:
+            a1 = np.array([[1.0+4.0j,2.0],
+                          [4.0+0.5j,5.0-3.0j],
+                          [7.0-2.0j,8.0+0.7j]], dtype=dtype)
+            b1 = np.array([16.0, 17.0+2.0j, 20.0-4.0j], dtype=dtype)
+            gelss, = get_lapack_funcs(('gelss',), (a1, b1))
+            #v,x,s,rank,work,info = cgelss(a,b,[cond,lwork,overwrite_a,overwrite_b])
+            v, x, s, rank, work, info = gelss(a1, b1,-1, -1, False, False)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+
+            v,x,s,rank,work,info = gelss(a1, b1,-1,lwork, False, False)
+            assert_allclose(x[:-1], np.array([1.161753632288328-1.901075709391912j, 1.735882340522193+1.521240901196909j],
+                            dtype=dtype), rtol=10*np.finfo(dtype).eps)
+            assert_allclose(s, np.array([13.035514762572043, 4.337666985231382], dtype=dtype), rtol=10*np.finfo(dtype).eps)
+
+    def test_gelsy(self):
+
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gelsy, = get_lapack_funcs(('gelsy',), (a1, b1))
+            #x,s,rank,work,iwork,info = *gelsd(a,b,lwork,size_iwork,[cond,overwrite_a,overwrite_b])
+            jptv = np.zeros((a1.shape[1],1), dtype=np.int32)
+            v,x,j,rank,work,info = gelsy(a1, b1, jptv, np.finfo(dtype).eps,-1)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+
+            jptv = np.zeros((a1.shape[1],1), dtype=np.int32)
+            v,x,j,rank,work,info = gelsy(a1, b1, jptv, np.finfo(dtype).eps, lwork, False, False)
+            assert_allclose(x[:-1], np.array([-14.333333333333323, 14.999999999999991], dtype=dtype),rtol=10*np.finfo(dtype).eps)
+            #assert_allclose(s, np.array([ 0.66666669,], dtype=dtype ) )
+
+        for dtype in COMPLEX_DTYPES:
+            a1 = np.array([[1.0+4.0j,2.0],
+                          [4.0+0.5j,5.0-3.0j],
+                          [7.0-2.0j,8.0+0.7j]], dtype=dtype)
+            b1 = np.array([16.0, 17.0+2.0j, 20.0-4.0j], dtype=dtype)
+            gelsy, = get_lapack_funcs(('gelsy',), (a1, b1))
+            #x,s,rank,work,rwork,iwork,info = cgelsd(a,b,lwork,size_rwork,size_iwork,[cond,overwrite_a,overwrite_b])
+            jptv = np.zeros((a1.shape[1],1), dtype=np.int32)
+            v,x,j,rank,work,rwork,info = gelsy(a1, b1, jptv, np.finfo(dtype).eps,-1)  # Request of sizes
+            lwork = work[0].real.astype(np.int)
+
+            jptv = np.zeros((a1.shape[1],1), dtype=np.int32)
+            v,x,j,rank,work,rwork,info = gelsy(a1, b1, jptv, np.finfo(dtype).eps, lwork, False, False)
+            assert_allclose(x[:-1], np.array([1.161753632288328-1.901075709391912j, 1.735882340522193+1.521240901196909j],
+                            dtype=dtype), rtol=10*np.finfo(dtype).eps)
+            #assert_allclose( s, np.array([ 106.12267169], dtype=dtype ) )
 
 class TestRegression(TestCase):