NumPy/SciPy LAPACK version
Hi all, This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module. The routines of most interest to me are: DGELSD DGGGLM DGGLSE I've also found that SciPy's sqrtm is broken:
a=matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) sqrtm(a) array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j, 3.8064764 +1.03420519e-08j], [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j, 5.3595916 -2.06841037e-08j], [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j, 6.9127068 +1.03420519e-08j]]) sqrtm(a)*sqrtm(a) array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j, 14.48926259 +7.87335527e-08j], [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j, 28.72522212 -2.21716697e-07j], [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j, 47.78551529 +1.42983144e-07j]])
So not even close... (and yes, it does deserve a bug report if one doesn't already exist) -Kevin
Hi, Did you try numpy.dot(sqrtm(a), sqrtm(a)) ? Matthieu 2007/7/16, Kevin Jacobs <jacobs@bioinformed.com> <bioinformed@gmail.com>:
Hi all,
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
I've also found that SciPy's sqrtm is broken:
a=matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) sqrtm(a) array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j, 3.8064764 +1.03420519e-08j], [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j, 5.3595916 -2.06841037e-08j], [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j, 6.9127068 +1.03420519e-08j]]) sqrtm(a)*sqrtm(a) array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j, 14.48926259 +7.87335527e-08j], [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j, 28.72522212 -2.21716697e-07j], [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j, 47.78551529 +1.42983144e-07j]])
So not even close... (and yes, it does deserve a bug report if one doesn't already exist)
-Kevin
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Oups, sorry, I missed the 'a=matrix'... 2007/7/16, Matthieu Brucher <matthieu.brucher@gmail.com>:
Hi,
Did you try numpy.dot(sqrtm(a), sqrtm(a)) ?
Matthieu
2007/7/16, Kevin Jacobs <jacobs@bioinformed.com> < bioinformed@gmail.com>:
Hi all,
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
I've also found that SciPy's sqrtm is broken:
a=matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) sqrtm(a) array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j, 3.8064764 +1.03420519e-08j], [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j, 5.3595916 -2.06841037e-08j], [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j, 6.9127068 +1.03420519e-08j]]) sqrtm(a)*sqrtm(a) array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j, 14.48926259 +7.87335527e-08j], [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j, 28.72522212 -2.21716697e-07j], [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j, 47.78551529 +1.42983144e-07j]])
So not even close... (and yes, it does deserve a bug report if one doesn't already exist)
-Kevin
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Kevin, the problem appears to be that sqrtm() gives back an array, rather than a matrix:
import scipy.linalg as sl a = s.matrix([[59, 64, 69],[64, 72, 80],[69, 80, 91]]) type(a) <class 'numpy.core.defmatrix.matrix'> a matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) a*a - N.dot(a,a) matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) b = sl.sqrtm(a) type(b) <type 'numpy.ndarray'> N.dot(b,b) array([[ 59. +1.85288457e-22j, 64. -6.61744490e-23j, 69. +1.85288457e-22j], [ 64. -2.64697796e-23j, 72. -3.70576914e-22j, 80. -5.29395592e-23j], [ 69. +1.85288457e-22j, 80. -1.32348898e-22j, 91. +2.38228016e-22j]]) b*b - N.dot(b,b) array([[-33.91512711 +1.03595922e-07j, -44.57413548 -1.82329475e-07j, -54.51073741 +7.87335527e-08j], [-44.57413548 -1.82329475e-07j, -48.15108664 +4.04046172e-07j, -51.27477788 -2.21716697e-07j], [-54.51073741 +7.87335527e-08j, -51.27477788 -2.21716697e-07j, -43.21448471 +1.42983144e-07j]])
This certainly is a slightly unexpected behaviour. Hanno "Kevin Jacobs <jacobs@bioinformed.com>" <bioinformed@gmail.com> said:
------=_Part_59405_32758974.1184593945795 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-Disposition: inline
Hi all,
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
I've also found that SciPy's sqrtm is broken:
a=matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) sqrtm(a) array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j, 3.8064764 +1.03420519e-08j], [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j, 5.3595916 -2.06841037e-08j], [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j, 6.9127068 +1.03420519e-08j]]) sqrtm(a)*sqrtm(a) array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j, 14.48926259 +7.87335527e-08j], [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j, 28.72522212 -2.21716697e-07j], [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j, 47.78551529 +1.42983144e-07j]])
So not even close... (and yes, it does deserve a bug report if one doesn't already exist)
-Kevin
------=_Part_59405_32758974.1184593945795 Content-Type: text/html; charset=ISO-8859-1 Content-Transfer-Encoding: 7bit Content-Disposition: inline
Hi all,<br><br>This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.<br><br>The routines of most interest to me are: <br>DGELSD<br>DGGGLM<br>DGGLSE<br><br>I've also found that SciPy's sqrtm is broken:<br><br>>>> a=matrix([[59, 64, 69],<br> [64, 72, 80],<br> [69, 80, 91]])<br>>>> sqrtm(a)<br>array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j,<br> 3.8064764 +1.03420519e-08j],<br> [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j,<br> 5.3595916 -2.06841037e-08j],<br> [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j,<br> 6.9127068 +1.03420519e-08j]])<br>>>> sqrtm(a)*sqrtm(a)<br>array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j,<br> 14.48926259 +7.87335527e-08j],<br> [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j,<br> 28.72522212 -2.21716697e-07j],<br> [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j,<br> 47.78551529 +1.42983144e-07j]])<br><br>So not even close... (and yes, it does deserve a bug report if one doesn't already exist)<br><br>-Kevin<br><br>
------=_Part_59405_32758974.1184593945795--
-- Hanno Klemm klemm@phys.ethz.ch
Mea culpa on the msqrt example, however I still think it is wrong to get a complex square-root back when a real valued result is expected and exists. -Kevin On 7/16/07, Hanno Klemm <klemm@phys.ethz.ch> wrote:
Kevin,
the problem appears to be that sqrtm() gives back an array, rather than a matrix:
import scipy.linalg as sl a = s.matrix([[59, 64, 69],[64, 72, 80],[69, 80, 91]]) type(a) <class 'numpy.core.defmatrix.matrix'> a matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) a*a - N.dot(a,a) matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) b = sl.sqrtm(a) type(b) <type 'numpy.ndarray'> N.dot(b,b) array([[ 59. +1.85288457e-22j, 64. -6.61744490e-23j, 69. +1.85288457e-22j], [ 64. -2.64697796e-23j, 72. -3.70576914e-22j, 80. -5.29395592e-23j], [ 69. +1.85288457e-22j, 80. -1.32348898e-22j, 91. +2.38228016e-22j]]) b*b - N.dot(b,b) array([[-33.91512711 +1.03595922e-07j, -44.57413548 -1.82329475e-07j, -54.51073741 +7.87335527e-08j], [-44.57413548 -1.82329475e-07j, -48.15108664 +4.04046172e-07j, -51.27477788 -2.21716697e-07j], [-54.51073741 +7.87335527e-08j, -51.27477788 -2.21716697e-07j, -43.21448471 +1.42983144e-07j]])
This certainly is a slightly unexpected behaviour.
Hanno
"Kevin Jacobs <jacobs@bioinformed.com>" <bioinformed@gmail.com> said:
------=_Part_59405_32758974.1184593945795 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-Disposition: inline
Hi all,
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
I've also found that SciPy's sqrtm is broken:
a=matrix([[59, 64, 69], [64, 72, 80], [69, 80, 91]]) sqrtm(a) array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j, 3.8064764 +1.03420519e-08j], [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j, 5.3595916 -2.06841037e-08j], [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j, 6.9127068 +1.03420519e-08j]]) sqrtm(a)*sqrtm(a) array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j, 14.48926259 +7.87335527e-08j], [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j, 28.72522212 -2.21716697e-07j], [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j, 47.78551529 +1.42983144e-07j]])
So not even close... (and yes, it does deserve a bug report if one doesn't already exist)
-Kevin
------=_Part_59405_32758974.1184593945795 Content-Type: text/html; charset=ISO-8859-1 Content-Transfer-Encoding: 7bit Content-Disposition: inline
Hi all,<br><br>This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.<br><br>The routines of most interest to me are: <br>DGELSD<br>DGGGLM<br>DGGLSE<br><br>I've also found that SciPy's sqrtm is broken:<br><br>>>> a=matrix([[59, 64, 69],<br> [64, 72, 80],<br> [69, 80, 91]])<br>>>> sqrtm(a)<br>array([[ 5.0084801 +1.03420519e-08j, 4.40747825 -2.06841037e-08j,<br> 3.8064764 +1.03420519e-08j],<br> [ 4.40747825 -2.06841037e-08j, 4.88353492 +4.13682075e-08j,<br> 5.3595916 -2.06841037e-08j],<br> [ 3.8064764 +1.03420519e-08j, 5.3595916 -2.06841037e-08j,<br> 6.9127068 +1.03420519e-08j]])<br>>>> sqrtm(a)*sqrtm(a)<br>array([[ 25.08487289 +1.03595922e-07j, 19.42586452 -1.82329475e-07j,<br> 14.48926259 +7.87335527e-08j],<br> [ 19.42586452 -1.82329475e-07j, 23.84891336 +4.04046172e-07j,<br> 28.72522212 -2.21716697e-07j],<br> [ 14.48926259 +7.87335527e-08j, 28.72522212 -2.21716697e-07j,<br> 47.78551529 +1.42983144e-07j]])<br><br>So not even close... (and yes, it does deserve a bug report if one doesn't already exist)<br><br>-Kevin<br><br>
------=_Part_59405_32758974.1184593945795--
-- Hanno Klemm klemm@phys.ethz.ch
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Kevin Jacobs <jacobs@bioinformed.com> wrote:
Mea culpa on the msqrt example, however I still think it is wrong to get a complex square-root back when a real valued result is expected and exists.
No, in floating point you accumulate error. Those 1e-22j's are part of the actual result. Some systems like MATLAB implicitly silent such small imaginary components; we don't. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On 7/16/07, Robert Kern <robert.kern@gmail.com> wrote:
Kevin Jacobs <jacobs@bioinformed.com> wrote:
Mea culpa on the msqrt example, however I still think it is wrong to get a complex square-root back when a real valued result is expected and exists.
No, in floating point you accumulate error. Those 1e-22j's are part of the actual result. Some systems like MATLAB implicitly silent such small imaginary components; we don't.
The problem is that the given matrix has a conditon number of about 10**17 and is almost singular. A singular value decomposition works fine, but apparently the sqrtm call suffers from roundoff and takes the sqrt of a negative number. Sqrtm returns real results in better conditioned cases. In [2]: sqrtm(eye(2)) Out[2]: array([[ 1., 0.], [ 0., 1.]]) Perhaps we aren't using the best method here. Chuck
On 7/16/07, Charles R Harris <charlesr.harris@gmail.com> wrote:
On 7/16/07, Robert Kern <robert.kern@gmail.com> wrote:
Kevin Jacobs <jacobs@bioinformed.com> wrote:
Mea culpa on the msqrt example, however I still think it is wrong to get a complex square-root back when a real valued result is expected and exists.
No, in floating point you accumulate error. Those 1e-22j's are part of the actual result. Some systems like MATLAB implicitly silent such small imaginary components; we don't.
The problem is that the given matrix has a conditon number of about 10**17 and is almost singular. A singular value decomposition works fine, but apparently the sqrtm call suffers from roundoff and takes the sqrt of a negative number. Sqrtm returns real results in better conditioned cases.
In [2]: sqrtm(eye(2)) Out[2]: array([[ 1., 0.], [ 0., 1.]])
Perhaps we aren't using the best method here.
Here is a better conditioned example:
a array([[ 1. , 0.5 , 0.3333, 0.25 ], [ 0.5 , 0.3333, 0.25 , 0.2 ], [ 0.3333, 0.25 , 0.2 , 0.1667], [ 0.25 , 0.2 , 0.1667, 0.1429]]) b=sqrtm(a) dot(b,b) array([[ 1. +0.j, 0.5 +0.j, 0.3333+0.j, 0.25 +0.j], [ 0.5 +0.j, 0.3333+0.j, 0.25 +0.j, 0.2 +0.j], [ 0.3333+0.j, 0.25 +0.j, 0.2 +0.j, 0.1667+0.j], [ 0.25 +0.j, 0.2 +0.j, 0.1667+0.j, 0.1429+0.j]]) dot(b,b)-a array([[ -1.99840144e-15+0.j, -9.43689571e-16+0.j, -5.55111512e-16+0.j, -5.55111512e-16+0.j], [ -1.05471187e-15+0.j, -5.55111512e-17+0.j, 5.55111512e-17+0.j, 0.00000000e+00+0.j], [ -6.66133815e-16+0.j, 1.11022302e-16+0.j, 1.66533454e-16+0.j, 1.11022302e-16+0.j], [ -5.55111512e-16+0.j, 1.11022302e-16+0.j, 1.38777878e-16+0.j, 8.32667268e-17+0.j]])
Also verified the results against svd and eigenvalue methods for computing msqrt. I suppose I'm just used to seeing msqrt() implemented completely using real valued arithmetic. -Kevin
On 7/16/07, Kevin Jacobs <jacobs@bioinformed.com> <bioinformed@gmail.com> wrote:
On 7/16/07, Charles R Harris <charlesr.harris@gmail.com> wrote:
On 7/16/07, Robert Kern < robert.kern@gmail.com> wrote:
Kevin Jacobs <jacobs@bioinformed.com> wrote:
Mea culpa on the msqrt example, however I still think it is wrong to get a complex square-root back when a real valued result is expected and exists.
No, in floating point you accumulate error. Those 1e-22j's are part of the actual result. Some systems like MATLAB implicitly silent such small imaginary components; we don't.
The problem is that the given matrix has a conditon number of about 10**17 and is almost singular. A singular value decomposition works fine, but apparently the sqrtm call suffers from roundoff and takes the sqrt of a negative number. Sqrtm returns real results in better conditioned cases.
In [2]: sqrtm(eye(2)) Out[2]: array([[ 1., 0.], [ 0., 1.]])
Perhaps we aren't using the best method here.
Here is a better conditioned example:
a array([[ 1. , 0.5 , 0.3333, 0.25 ], [ 0.5 , 0.3333, 0.25 , 0.2 ], [ 0.3333, 0.25 , 0.2 , 0.1667], [ 0.25 , 0.2 , 0.1667, 0.1429]]) b=sqrtm(a) dot(b,b) array([[ 1. +0.j, 0.5 +0.j, 0.3333+0.j, 0.25 +0.j], [ 0.5 +0.j, 0.3333+0.j, 0.25 +0.j, 0.2 +0.j], [ 0.3333+0.j, 0.25 +0.j, 0.2 +0.j, 0.1667+0.j], [ 0.25 +0.j, 0.2 +0.j, 0.1667+0.j, 0.1429+0.j]]) dot(b,b)-a array([[ -1.99840144e-15+0.j, -9.43689571e-16+0.j, -5.55111512e-16+0.j, -5.55111512e-16+0.j], [ -1.05471187e-15+0.j, -5.55111512e-17+0.j , 5.55111512e-17+0.j, 0.00000000e+00+0.j], [ -6.66133815e-16+0.j, 1.11022302e-16+0.j, 1.66533454e-16+0.j, 1.11022302e-16+0.j], [ -5.55111512e-16+0.j, 1.11022302e-16+0.j , 1.38777878e-16+0.j, 8.32667268e-17+0.j]])
Also verified the results against svd and eigenvalue methods for computing msqrt. I suppose I'm just used to seeing msqrt() implemented completely using real valued arithmetic.
Hmm, I get a real result for this, although the result is wildly incorrect. Sqrtm isn't part of numpy, where are you getting it from? Mine is coming from pylab and looks remarkably buggy. Chuck
On 7/16/07, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hmm,
I get a real result for this, although the result is wildly incorrect. Sqrtm isn't part of numpy, where are you getting it from? Mine is coming from pylab and looks remarkably buggy.
from scipy.linalg import sqrtm I'm posting on the NumPy list since I already subscribe here, although my questions are more topical for the SciPy folks. I'm hoping the overlap is considerable and I'll be forgiven (since adding yet another mailing list is sure to capsize my poor mailbox). Thanks, -Kevin
On 7/16/07, Kevin Jacobs <jacobs@bioinformed.com> <bioinformed@gmail.com> wrote:
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
STEGR Thanks for all of the feedback on sqrtm. Can anyone comment on the suitability of adding LAPACK 3.0 functions to SciPy.LinAlg? I need to do the work regardless, but being able to contribute it back would be very nice. -Kevin
Kevin Jacobs <jacobs@bioinformed.com> wrote:
On 7/16/07, *Kevin Jacobs <jacobs@bioinformed.com <mailto:jacobs@bioinformed.com>>* <bioinformed@gmail.com <mailto:bioinformed@gmail.com>> wrote:
This is a bit of a SciPy question, but I thought I'd ask here since I'm already subscribed. I'd like to add some new LAPACK bindings to SciPy and was wondering if there was a minimum version requirement for LAPACK, since it would be ideal if I could use some of the newer 3.0 features. In addition to using some block methods only added in 3.0, it is very convenient to use the WORK=-1 for space queries instead of reimplementing the underlying logic in the calc_work module.
The routines of most interest to me are: DGELSD DGGGLM DGGLSE
STEGR
Thanks for all of the feedback on sqrtm. Can anyone comment on the suitability of adding LAPACK 3.0 functions to SciPy.LinAlg? I need to do the work regardless, but being able to contribute it back would be very nice.
And we'd certainly appreciate the contribution. I'm tentatively going to say yes, we should start requiring LAPACK 3.0 unless if there is some very important platform that only comes with an older LAPACK. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On 7/16/07, Robert Kern <robert.kern@gmail.com> wrote:
And we'd certainly appreciate the contribution. I'm tentatively going to say yes, we should start requiring LAPACK 3.0 unless if there is some very important platform that only comes with an older LAPACK.
Great! The added benefit is that the calc_work module effectively goes away, at least in its current form. I've checked that the current versions of the major platform-optimized math libraries all use 3.0 or greater, including Sun's math performance library, Intel MKL, and AMD's ACML. If one is willing to forgo vendor-supplied platform optimizations, then ATLAS+LAPACK or just Netlib LAPACK should suffice for many users. -Kevin
participants (5)
-
Charles R Harris -
Hanno Klemm -
Kevin Jacobs <jacobs＠bioinformed.com> -
Matthieu Brucher -
Robert Kern