From 72a0c5631042ae16e3509dcd626a3a3544b28c4d Mon Sep 17 00:00:00 2001
From: Fabian Pedregosa <fabian.pedregosa@inria.fr>
Date: Mon, 22 Nov 2010 14:56:37 +0100
Subject: [PATCH] Add support for np.float32 matrices in lars_path.
Some single-precision routines have had to be ported from cblas. Later
on we could probably remove most of them since new scipy has
triangular_solve function and cholesky_delete could be implemented in
Python.
---
scikits/learn/src/cblas/ATL_srefcopy.c | 148 +++
scikits/learn/src/cblas/ATL_srefrot.c | 170 ++++
scikits/learn/src/cblas/ATL_srefrotg.c | 146 +++
scikits/learn/src/cblas/ATL_sreftrsv.c | 213 ++++
scikits/learn/src/cblas/ATL_sreftrsvLNN.c | 94 ++
scikits/learn/src/cblas/ATL_sreftrsvLNU.c | 94 ++
scikits/learn/src/cblas/ATL_sreftrsvLTN.c | 96 ++
scikits/learn/src/cblas/ATL_sreftrsvLTU.c | 96 ++
scikits/learn/src/cblas/ATL_sreftrsvUNN.c | 95 ++
scikits/learn/src/cblas/ATL_sreftrsvUNU.c | 95 ++
scikits/learn/src/cblas/ATL_sreftrsvUTN.c | 95 ++
scikits/learn/src/cblas/ATL_sreftrsvUTU.c | 95 ++
scikits/learn/src/cblas/atlas_ssysinfo.h | 0
scikits/learn/src/cblas/cblas_scopy.c | 52 +
scikits/learn/src/cblas/cblas_srot.c | 60 ++
scikits/learn/src/cblas/cblas_srotg.c | 42 +
scikits/learn/src/cblas/cblas_strsv.c | 87 ++
scikits/learn/utils/arrayfuncs.c | 1100 +++++++--------------
scikits/learn/utils/arrayfuncs.pyx | 50 +-
scikits/learn/utils/src/cholesky_delete.c | 36 +-
20 files changed, 2125 insertions(+), 739 deletions(-)
create mode 100644 scikits/learn/src/cblas/ATL_srefcopy.c
create mode 100644 scikits/learn/src/cblas/ATL_srefrot.c
create mode 100644 scikits/learn/src/cblas/ATL_srefrotg.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsv.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvLNN.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvLNU.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvLTN.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvLTU.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvUNN.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvUNU.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvUTN.c
create mode 100644 scikits/learn/src/cblas/ATL_sreftrsvUTU.c
create mode 100644 scikits/learn/src/cblas/atlas_ssysinfo.h
create mode 100644 scikits/learn/src/cblas/cblas_scopy.c
create mode 100644 scikits/learn/src/cblas/cblas_srot.c
create mode 100644 scikits/learn/src/cblas/cblas_srotg.c
create mode 100644 scikits/learn/src/cblas/cblas_strsv.c
diff --git a/scikits/learn/src/cblas/ATL_srefcopy.c b/scikits/learn/src/cblas/ATL_srefcopy.c
new file mode 100644
index 0000000000..fec830e4fc
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_srefcopy.c
@@ -0,0 +1,148 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflevel1.h"
+
+void ATL_srefcopy
+(
+ const int N,
+ const float * X,
+ const int INCX,
+ float * Y,
+ const int INCY
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_srefcopy copies the entries of an n-vector x into an n-vector y.
+ *
+ * Arguments
+ * =========
+ *
+ * N (input) const int
+ * On entry, N specifies the length of the vector x. N must be
+ * at least zero. Unchanged on exit.
+ *
+ * X (input) const float *
+ * On entry, X points to the first entry to be accessed of an
+ * incremented array of size equal to or greater than
+ * ( 1 + ( n - 1 ) * abs( INCX ) ) * sizeof( float ),
+ * that contains the vector x. Unchanged on exit.
+ *
+ * INCX (input) const int
+ * On entry, INCX specifies the increment for the elements of X.
+ * INCX must not be zero. Unchanged on exit.
+ *
+ * Y (input/output) float *
+ * On entry, Y points to the first entry to be accessed of an
+ * incremented array of size equal to or greater than
+ * ( 1 + ( n - 1 ) * abs( INCY ) ) * sizeof( float ),
+ * that contains the vector y. On exit, the entries of the in-
+ * cremented array X are copied into the entries of the incre-
+ * mented array Y.
+ *
+ * INCY (input) const int
+ * On entry, INCY specifies the increment for the elements of Y.
+ * INCY must not be zero. Unchanged on exit.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float x0, x1, x2, x3, x4, x5, x6, x7;
+ float * StX;
+ register int i;
+ int nu;
+ const int incX2 = 2 * INCX, incY2 = 2 * INCY,
+ incX3 = 3 * INCX, incY3 = 3 * INCY,
+ incX4 = 4 * INCX, incY4 = 4 * INCY,
+ incX5 = 5 * INCX, incY5 = 5 * INCY,
+ incX6 = 6 * INCX, incY6 = 6 * INCY,
+ incX7 = 7 * INCX, incY7 = 7 * INCY,
+ incX8 = 8 * INCX, incY8 = 8 * INCY;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ if( N > 0 )
+ {
+ if( ( nu = ( N >> 3 ) << 3 ) != 0 )
+ {
+ StX = (float *)X + nu * INCX;
+
+ do
+ {
+ x0 = (*X); x4 = X[incX4]; x1 = X[INCX ]; x5 = X[incX5];
+ x2 = X[incX2]; x6 = X[incX6]; x3 = X[incX3]; x7 = X[incX7];
+
+ *Y = x0; Y[incY4] = x4; Y[INCY ] = x1; Y[incY5] = x5;
+ Y[incY2] = x2; Y[incY6] = x6; Y[incY3] = x3; Y[incY7] = x7;
+
+ X += incX8;
+ Y += incY8;
+
+ } while( X != StX );
+ }
+
+ for( i = N - nu; i != 0; i-- )
+ {
+ x0 = (*X);
+ *Y = x0;
+
+ X += INCX;
+ Y += INCY;
+ }
+ }
+/*
+ * End of ATL_srefcopy
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_srefrot.c b/scikits/learn/src/cblas/ATL_srefrot.c
new file mode 100644
index 0000000000..e1a1d948ff
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_srefrot.c
@@ -0,0 +1,170 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflevel1.h"
+
+void ATL_srefrot
+(
+ const int N,
+ float * X,
+ const int INCX,
+ float * Y,
+ const int INCY,
+ const float C,
+ const float S
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_srefrot applies a plane rotation to the two n-vectors x and y.
+ * This routine computes:
+ *
+ * [ x_i ] [ c s ] [ x_i ]
+ * [ y_i ] = [ -s c ] [ y_i ] for all i = 1 .. n.
+ *
+ * If n <= 0 or if c = 1 and s = 0, this subroutine returns immediately.
+ *
+ * Arguments
+ * =========
+ *
+ * N (input) const int
+ * On entry, N specifies the length of the vector x. N must be
+ * at least zero. Unchanged on exit.
+ *
+ * X (input/output) float *
+ * On entry, X points to the first entry to be accessed of an
+ * incremented array of size equal to or greater than
+ * ( 1 + ( n - 1 ) * abs( INCX ) ) * sizeof( float ),
+ * that contains the vector x. On exit, the entries of the in-
+ * cremented array X are rotated with the entries of the incre-
+ * mented array Y.
+ *
+ * INCX (input) const int
+ * On entry, INCX specifies the increment for the elements of X.
+ * INCX must not be zero. Unchanged on exit.
+ *
+ * Y (input/output) float *
+ * On entry, Y points to the first entry to be accessed of an
+ * incremented array of size equal to or greater than
+ * ( 1 + ( n - 1 ) * abs( INCY ) ) * sizeof( float ),
+ * that contains the vector y. On exit, the entries of the in-
+ * cremented array Y are rotated with the entries of the incre-
+ * mented array X.
+ *
+ * INCY (input) const int
+ * On entry, INCY specifies the increment for the elements of Y.
+ * INCY must not be zero. Unchanged on exit.
+ *
+ * C (input) const float
+ * On entry, C specifies the scalar c definiting the plane rota-
+ * tion. Unchanged on exit.
+ *
+ * S (input) const float
+ * On entry, S specifies the scalar s definiting the plane rota-
+ * tion. Unchanged on exit.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float x0, x1, x2, x3, y0, y1, y2, y3;
+ register const float co = C, si = S;
+ float * StX;
+ register int i;
+ int nu;
+ const int incX2 = 2 * INCX, incY2 = 2 * INCY,
+ incX3 = 3 * INCX, incY3 = 3 * INCY,
+ incX4 = 4 * INCX, incY4 = 4 * INCY;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ if( ( N > 0 ) && !Msone( co, si ) )
+ {
+ if( ( nu = ( N >> 2 ) << 2 ) != 0 )
+ {
+ StX = (float *)X + nu * INCX;
+
+ do
+ {
+ x0 = (*X); y0 = (*Y);
+ x1 = X[INCX ]; y1 = Y[INCY ];
+ x2 = X[incX2]; y2 = Y[incY2];
+ x3 = X[incX3]; y3 = Y[incY3];
+
+ *X = co * x0 + si * y0; *Y = co * y0 - si * x0;
+ X[INCX ] = co * x1 + si * y1; Y[INCY ] = co * y1 - si * x1;
+ X[incX2] = co * x2 + si * y2; Y[incY2] = co * y2 - si * x2;
+ X[incX3] = co * x3 + si * y3; Y[incY3] = co * y3 - si * x3;
+
+ X += incX4;
+ Y += incY4;
+
+ } while( X != StX );
+ }
+
+ for( i = N - nu; i != 0; i-- )
+ {
+ x0 = (*X);
+ y0 = (*Y);
+
+ *X = co * x0 + si * y0;
+ *Y = co * y0 - si * x0;
+
+ X += INCX;
+ Y += INCY;
+ }
+ }
+/*
+ * End of ATL_srefrot
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_srefrotg.c b/scikits/learn/src/cblas/ATL_srefrotg.c
new file mode 100644
index 0000000000..4692a03754
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_srefrotg.c
@@ -0,0 +1,146 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflevel1.h"
+
+void ATL_srefrotg
+(
+ float * A,
+ float * B,
+ float * C,
+ float * S
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_srefrotg constructs a Givens plane rotation. Given the scalars a
+ * and b, this routine computes the following quantities:
+ *
+ * sigma = sgn(a) if |a| > |b|, sgn(b) otherwise;
+ * r = sigma * sqrt( a^2 + b^2 );
+ * c = a / r if r <> 0, 1 otherwise;
+ * s = b / r if r <> 0, 0 otherwise.
+ *
+ * The numbers c, s and r then satisfy the matrix equation:
+ *
+ * [ c s ] [ a ] [ r ]
+ * [ -s c ] [ b ] = [ 0 ].
+ *
+ * The introduction of sigma is not essential to the computation of a
+ * Givens rotation matrix, but it permits later stable reconstruction of
+ * c and s from just one stored number. For this purpose, this routine
+ * also computes
+ *
+ * s if |a| > |b|,
+ * z = 1 / c if |b| >= |a| and c <> 0,
+ * 1 if c = 0.
+ *
+ * This subroutine returns r overwriting a, and z overwriting b, as well
+ * as returning c and s. If one later wishes to reconstruct c and s from
+ * z, it can be done as follows:
+ *
+ * if z = 1, set c = 0 and s = 1,
+ * if |z| < 1, set c = sqrt(1 - z^2) and s = z,
+ * if |z| > 1, set c = 1 / z and s = sqrt(1 - c^2).
+ *
+ * See ``Basic Linear Algebra Subprograms for Fortran Usage'' by C. Law-
+ * son, R. Hanson, D. Kincaid and F. Krogh, ACM Transactions on Mathema-
+ * tical Software, 1979, 5(3) pp 308-323, for further information.
+ *
+ * Arguments
+ * =========
+ *
+ * A (input/output) float *
+ * On entry, A specifies the scalar a. On exit, A is overwritten
+ * by the scalar r defined above.
+ *
+ * B (input/output) float *
+ * On entry, B specifies the scalar b. On exit, B is overwritten
+ * by the scalar z defined above.
+ *
+ * C (output) float *
+ * On exit, C specifies the scalar c defined above.
+ *
+ * S (output) float *
+ * On exit, S specifies the scalar s defined above.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float absa, absb, roe, scale, r, tmpa, tmpb, z;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ absa = Msabs( *A ); absb = Msabs( *B );
+ roe = ( absa > absb ? (*A) : (*B) ); scale = absa + absb;
+
+ if( scale != ATL_sZERO )
+ {
+ tmpa = (*A) / scale; tmpb = (*B) / scale;
+ if( roe < ATL_sZERO )
+ { r = - scale * sqrt( ( tmpa * tmpa ) + ( tmpb * tmpb ) ); }
+ else
+ { r = scale * sqrt( ( tmpa * tmpa ) + ( tmpb * tmpb ) ); }
+ *C = (*A) / r; *S = (*B) / r; z = ATL_sONE;
+ if( absa > absb ) { z = *S; }
+ if( ( absb >= absa ) && ( (*C) != ATL_sZERO ) )
+ { z = ATL_sONE / (*C); }
+ }
+ else { *C = ATL_sONE; *S = r = z = ATL_sZERO; }
+
+ *A = r; *B = z;
+/*
+ * End of ATL_srefrotg
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsv.c b/scikits/learn/src/cblas/ATL_sreftrsv.c
new file mode 100644
index 0000000000..f77d179215
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsv.c
@@ -0,0 +1,213 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsv
+(
+ const enum ATLAS_UPLO UPLO,
+ const enum ATLAS_TRANS TRANS,
+ const enum ATLAS_DIAG DIAG,
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsv solves one of the systems of equations
+ *
+ * A * x = b, or A' * x = b,
+ *
+ * where b and x are n-element vectors and A is an n by n unit, or non-
+ * unit, upper or lower triangular matrix.
+ *
+ * No test for singularity or near-singularity is included in this
+ * routine. Such tests must be performed before calling this routine.
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) const enum ATLAS_UPLO
+ * On entry, UPLO specifies whether the matrix is an upper or
+ * lower triangular matrix as follows:
+ *
+ * UPLO = AtlasUpper A is an upper triangular matrix.
+ *
+ * UPLO = AtlasLower A is a lower triangular matrix.
+ *
+ * Unchanged on exit.
+ *
+ * TRANS (input) const enum ATLAS_TRANS
+ * On entry, TRANS specifies the equations to be solved as fol-
+ * lows:
+ *
+ * TRANS = AtlasNoTrans A * x = b,
+ *
+ * TRANS = AtlasConj A * x = b,
+ *
+ * TRANS = AtlasTrans A' * x = b,
+ *
+ * TRANS = AtlasTrans A' * x = b.
+ *
+ * Unchanged on exit.
+ *
+ * DIAG (input) const enum ATLAS_DIAG
+ * On entry, DIAG specifies whether or not A is unit triangu-
+ * lar as follows:
+ *
+ * DIAG = AtlasUnit A is assumed to be unit triangular,
+ *
+ * DIAG = AtlasNonUnit A is not assumed to be unit trian-
+ * gular.
+ *
+ * Unchanged on exit.
+ *
+ * N (input) const int
+ * On entry, N specifies the order of the matrix A. N must be at
+ * least zero. Unchanged on exit.
+ *
+ * A (input) const float *
+ * On entry, A points to an array of size equal to or greater
+ * than LDA * n * sizeof( float ). Before entry with
+ * UPLO = AtlasUpper, the leading n by n upper triangular part
+ * of the array A must contain the upper triangular matrix and
+ * the strictly lower triangular part of A is not referenced.
+ * Before entry with UPLO = AtlasLower, the leading n by n lower
+ * triangular part of the array A must contain the lower trian-
+ * gular matrix and the strictly upper triangular part of A is
+ * not referenced. Unchanged on exit.
+ *
+ * Note that when DIAG = AtlasUnit, the diagonal elements of A
+ * are not referenced either, but are assumed to be unity.
+ *
+ * LDA (input) const int
+ * On entry, LDA specifies the leading dimension of A as decla-
+ * red in the calling (sub) program. LDA must be at least
+ * MAX( 1, n ). Unchanged on exit.
+ *
+ * X (input/output) float *
+ * On entry, X points to the first entry to be accessed of an
+ * incremented array of size equal to or greater than
+ * ( 1 + ( n - 1 ) * abs( INCX ) ) * sizeof( float ),
+ * that contains the vector x. Before entry, the incremented ar-
+ * ray X must contain the n element right-hand side vector b. On
+ * exit, X is overwritten with the solution vector x.
+ *
+ * INCX (input) const int
+ * On entry, INCX specifies the increment for the elements of X.
+ * INCX must not be zero. Unchanged on exit.
+ *
+ * ---------------------------------------------------------------------
+ */
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ if( N == 0 ) return;
+
+ if( UPLO == AtlasUpper )
+ {
+ if( ( TRANS == AtlasNoTrans ) || ( TRANS == AtlasConj ) )
+ {
+ if( DIAG == AtlasNonUnit )
+ {
+ ATL_sreftrsvUNN( N, A, LDA, X, INCX );
+ }
+ else
+ {
+ ATL_sreftrsvUNU( N, A, LDA, X, INCX );
+ }
+ }
+ else
+ {
+ if( DIAG == AtlasNonUnit )
+ {
+ ATL_sreftrsvUTN( N, A, LDA, X, INCX );
+ }
+ else
+ {
+ ATL_sreftrsvUTU( N, A, LDA, X, INCX );
+ }
+ }
+ }
+ else
+ {
+ if( ( TRANS == AtlasNoTrans ) || ( TRANS == AtlasConj ) )
+ {
+ if( DIAG == AtlasNonUnit )
+ {
+ ATL_sreftrsvLNN( N, A, LDA, X, INCX );
+ }
+ else
+ {
+ ATL_sreftrsvLNU( N, A, LDA, X, INCX );
+ }
+ }
+ else
+ {
+ if( DIAG == AtlasNonUnit )
+ {
+ ATL_sreftrsvLTN( N, A, LDA, X, INCX );
+ }
+ else
+ {
+ ATL_sreftrsvLTU( N, A, LDA, X, INCX );
+ }
+ }
+ }
+/*
+ * End of ATL_sreftrsv
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvLNN.c b/scikits/learn/src/cblas/ATL_sreftrsvLNN.c
new file mode 100644
index 0000000000..109a688848
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvLNN.c
@@ -0,0 +1,94 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvLNN
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvLNN( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasLower, AtlasNoTrans, AtlasNonUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx, ldap1 = LDA + 1;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = 0, jaj = 0, jx = 0; j < N; j++, jaj += ldap1, jx += INCX )
+ {
+ X[jx] /= A[jaj]; t0 = X[jx];
+ for( i = j+1, iaij = jaj+1, ix = jx + INCX;
+ i < N; i++, iaij += 1, ix += INCX ) { X[ix] -= t0 * A[iaij]; }
+ }
+/*
+ * End of ATL_sreftrsvLNN
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvLNU.c b/scikits/learn/src/cblas/ATL_sreftrsvLNU.c
new file mode 100644
index 0000000000..f97e51be13
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvLNU.c
@@ -0,0 +1,94 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvLNU
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvLNU( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasLower, AtlasNoTrans, AtlasUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx, ldap1 = LDA + 1;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = 0, jaj = 0, jx = 0; j < N; j++, jaj += ldap1, jx += INCX )
+ {
+ t0 = X[jx];
+ for( i = j+1, iaij = jaj+1, ix = jx + INCX;
+ i < N; i++, iaij += 1, ix += INCX ) { X[ix] -= t0 * A[iaij]; }
+ }
+/*
+ * End of ATL_sreftrsvLNU
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvLTN.c b/scikits/learn/src/cblas/ATL_sreftrsvLTN.c
new file mode 100644
index 0000000000..923f84e709
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvLTN.c
@@ -0,0 +1,96 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvLTN
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvLTN( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasLower, AtlasTrans, AtlasNonUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx, ldap1 = LDA + 1;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = N-1, jaj = (N-1)*(ldap1), jx = (N-1)*INCX;
+ j >= 0; j--, jaj -= ldap1, jx -= INCX )
+ {
+ t0 = X[jx];
+ for( i = j+1, iaij = 1+jaj, ix = jx + INCX;
+ i < N; i++, iaij += 1, ix += INCX ) { t0 -= A[iaij] * X[ix]; }
+ t0 /= A[jaj]; X[jx] = t0;
+ }
+/*
+ * End of ATL_sreftrsvLTN
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvLTU.c b/scikits/learn/src/cblas/ATL_sreftrsvLTU.c
new file mode 100644
index 0000000000..40288f92aa
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvLTU.c
@@ -0,0 +1,96 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvLTU
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvLTU( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasLower, AtlasTrans, AtlasUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx, ldap1 = LDA + 1;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = N-1, jaj = (N-1)*(ldap1), jx = (N-1)*INCX;
+ j >= 0; j--, jaj -= ldap1, jx -= INCX )
+ {
+ t0 = X[jx];
+ for( i = j+1, iaij = 1+jaj, ix = jx + INCX;
+ i < N; i++, iaij += 1, ix += INCX ) { t0 -= A[iaij] * X[ix]; }
+ X[jx] = t0;
+ }
+/*
+ * End of ATL_sreftrsvLTU
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvUNN.c b/scikits/learn/src/cblas/ATL_sreftrsvUNN.c
new file mode 100644
index 0000000000..6e4b465d09
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvUNN.c
@@ -0,0 +1,95 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvUNN
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvUNN( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasUpper, AtlasNoTrans, AtlasNonUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = N-1, jaj = (N-1)*LDA, jx = (N-1)*INCX;
+ j >= 0; j--, jaj -= LDA, jx -= INCX )
+ {
+ X[jx] /= A[j+jaj]; t0 = X[jx];
+ for( i = 0, iaij = jaj, ix = 0; i < j; i++, iaij += 1, ix += INCX )
+ { X[ix] -= t0 * A[iaij]; }
+ }
+/*
+ * End of ATL_sreftrsvUNN
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvUNU.c b/scikits/learn/src/cblas/ATL_sreftrsvUNU.c
new file mode 100644
index 0000000000..a3c2acaf76
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvUNU.c
@@ -0,0 +1,95 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvUNU
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvUNU( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasUpper, AtlasNoTrans, AtlasUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = N-1, jaj = (N-1)*LDA, jx = (N-1)*INCX;
+ j >= 0; j--, jaj -= LDA, jx -= INCX )
+ {
+ t0 = X[jx];
+ for( i = 0, iaij = jaj, ix = 0; i < j; i++, iaij += 1, ix += INCX )
+ { X[ix] -= t0 * A[iaij]; }
+ }
+/*
+ * End of ATL_sreftrsvUNU
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvUTN.c b/scikits/learn/src/cblas/ATL_sreftrsvUTN.c
new file mode 100644
index 0000000000..6a689fa0eb
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvUTN.c
@@ -0,0 +1,95 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvUTN
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvUTN( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasUpper, AtlasTrans, AtlasNonUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = 0, jaj = 0,jx = 0; j < N; j++, jaj += LDA, jx += INCX )
+ {
+ t0 = X[jx];
+ for( i = 0, iaij = jaj, ix = 0; i < j; i++, iaij += 1, ix += INCX )
+ { t0 -= A[iaij] * X[ix]; }
+ t0 /= A[iaij]; X[jx] = t0;
+ }
+/*
+ * End of ATL_sreftrsvUTN
+ */
+}
diff --git a/scikits/learn/src/cblas/ATL_sreftrsvUTU.c b/scikits/learn/src/cblas/ATL_sreftrsvUTU.c
new file mode 100644
index 0000000000..244700aa6e
--- /dev/null
+++ b/scikits/learn/src/cblas/ATL_sreftrsvUTU.c
@@ -0,0 +1,95 @@
+/* ---------------------------------------------------------------------
+ *
+ * -- Automatically Tuned Linear Algebra Software (ATLAS)
+ * (C) Copyright 2000 All Rights Reserved
+ *
+ * -- ATLAS routine -- Version 3.9.24 -- December 25, 2000
+ *
+ * Author : Antoine P. Petitet
+ * Originally developed at the University of Tennessee,
+ * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * -- Copyright notice and Licensing terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in
+ * the documentation and/or other materials provided with the distri-
+ * bution.
+ * 3. The name of the University, the ATLAS group, or the names of its
+ * contributors may not be used to endorse or promote products deri-
+ * ved from this software without specific written permission.
+ *
+ * -- Disclaimer:
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
+ * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
+ * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
+ * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * Include files
+ */
+#include "atlas_refmisc.h"
+#include "atlas_reflvl2.h"
+#include "atlas_reflevel2.h"
+
+void ATL_sreftrsvUTU
+(
+ const int N,
+ const float * A,
+ const int LDA,
+ float * X,
+ const int INCX
+)
+{
+/*
+ * Purpose
+ * =======
+ *
+ * ATL_sreftrsvUTU( ... )
+ *
+ * <=>
+ *
+ * ATL_sreftrsv( AtlasUpper, AtlasTrans, AtlasUnit, ... )
+ *
+ * See ATL_sreftrsv for details.
+ *
+ * ---------------------------------------------------------------------
+ */
+/*
+ * .. Local Variables ..
+ */
+ register float t0;
+ int i, iaij, ix, j, jaj, jx;
+/* ..
+ * .. Executable Statements ..
+ *
+ */
+ for( j = 0, jaj = 0, jx = 0; j < N; j++, jaj += LDA, jx += INCX )
+ {
+ t0 = X[jx];
+ for( i = 0, iaij = jaj, ix = 0; i < j; i++, iaij += 1, ix += INCX )
+ { t0 -= A[iaij] * X[ix]; }
+ X[jx] = t0;
+ }
+/*
+ * End of ATL_sreftrsvUTU
+ */
+}
diff --git a/scikits/learn/src/cblas/atlas_ssysinfo.h b/scikits/learn/src/cblas/atlas_ssysinfo.h
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/scikits/learn/src/cblas/cblas_scopy.c b/scikits/learn/src/cblas/cblas_scopy.c
new file mode 100644
index 0000000000..38feca50a6
--- /dev/null
+++ b/scikits/learn/src/cblas/cblas_scopy.c
@@ -0,0 +1,52 @@
+/*
+ * Automatically Tuned Linear Algebra Software v3.9.25
+ * (C) Copyright 1999 R. Clint Whaley
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * 3. The name of the ATLAS group or the names of its contributers may
+ * not be used to endorse or promote products derived from this
+ * software without specific written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+
+#define SREAL
+#include "atlas_misc.h"
+#ifdef ATL_USEPTHREADS
+ #include "atlas_ptalias1.h"
+#endif
+#include "atlas_level1.h"
+#include "cblas.h"
+
+void cblas_scopy(const int N, const float *X, const int incX,
+ float *Y, const int incY)
+{
+ if (N > 0)
+ {
+ if (incX < 0)
+ {
+ if (incY < 0) ATL_scopy(N, X, -incX, Y, -incY);
+ else ATL_scopy(N, X+(1-N)*incX, incX, Y, incY);
+ }
+ else if (incY < 0) ATL_scopy(N, X+(N-1)*incX, -incX, Y, -incY);
+ else ATL_scopy(N, X, incX, Y, incY);
+ }
+}
diff --git a/scikits/learn/src/cblas/cblas_srot.c b/scikits/learn/src/cblas/cblas_srot.c
new file mode 100644
index 0000000000..b2365d3866
--- /dev/null
+++ b/scikits/learn/src/cblas/cblas_srot.c
@@ -0,0 +1,60 @@
+/*
+ * Automatically Tuned Linear Algebra Software v3.9.25
+ * (C) Copyright 1999 R. Clint Whaley
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * 3. The name of the ATLAS group or the names of its contributers may
+ * not be used to endorse or promote products derived from this
+ * software without specific written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+
+#define SREAL
+#include "atlas_misc.h"
+#ifdef ATL_USEPTHREADS
+ #include "atlas_ptalias1.h"
+#endif
+#include "atlas_level1.h"
+#include "cblas.h"
+
+void cblas_srot(const int N, float *X, const int incX,
+ float *Y, const int incY, const float c, const float s)
+{
+ float *x = X, *y = Y;
+ int incx = incX, incy = incY;
+
+ if (N > 0)
+ {
+ if (incX < 0)
+ {
+ if (incY < 0) { incx = -incx; incy = -incy; }
+ else x += -incX * ((N-1));
+ }
+ else if (incY < 0)
+ {
+ incy = -incy;
+ incx = -incx;
+ x += (N-1)*(incX);
+ }
+ ATL_srot(N, x, incx, y, incy, c, s);
+ }
+}
diff --git a/scikits/learn/src/cblas/cblas_srotg.c b/scikits/learn/src/cblas/cblas_srotg.c
new file mode 100644
index 0000000000..56bb01757c
--- /dev/null
+++ b/scikits/learn/src/cblas/cblas_srotg.c
@@ -0,0 +1,42 @@
+/*
+ * Automatically Tuned Linear Algebra Software v3.9.25
+ * (C) Copyright 1999 R. Clint Whaley
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * 3. The name of the ATLAS group or the names of its contributers may
+ * not be used to endorse or promote products derived from this
+ * software without specific written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+
+#define SREAL
+#include "atlas_misc.h"
+#ifdef ATL_USEPTHREADS
+ #include "atlas_ptalias1.h"
+#endif
+#include "atlas_level1.h"
+#include "cblas.h"
+
+void cblas_srotg(float *a, float *b, float *c, float *s)
+{
+ Mjoin(PATL,rotg)(a, b, c, s);
+}
diff --git a/scikits/learn/src/cblas/cblas_strsv.c b/scikits/learn/src/cblas/cblas_strsv.c
new file mode 100644
index 0000000000..367ffded02
--- /dev/null
+++ b/scikits/learn/src/cblas/cblas_strsv.c
@@ -0,0 +1,87 @@
+/*
+ * Automatically Tuned Linear Algebra Software v3.9.25
+ * (C) Copyright 1999 R. Clint Whaley
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions, and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * 3. The name of the ATLAS group or the names of its contributers may
+ * not be used to endorse or promote products derived from this
+ * software without specific written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+
+#define SREAL
+#include "atlas_misc.h"
+#include "cblas.h"
+#ifdef ATL_USEPTHREADS
+ #include "atlas_ptalias2.h"
+#endif
+#include "atlas_level2.h"
+
+void cblas_strsv(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
+ const enum CBLAS_TRANSPOSE TA, const enum CBLAS_DIAG Diag,
+ const int N, const float *A, const int lda,
+ float *X, const int incX)
+{
+ int info = 2000;
+ enum CBLAS_UPLO uplo;
+ enum CBLAS_TRANSPOSE ta;
+ #define x X
+
+#ifndef NoCblasErrorChecks
+ if (Order != CblasColMajor && Order != CblasRowMajor)
+ info = cblas_errprn(1, info, "Order must be %d or %d, but is set to %d",
+ CblasRowMajor, CblasColMajor, Order);
+ if (Uplo != CblasUpper && Uplo != CblasLower)
+ info = cblas_errprn(2, info, "UPLO must be %d or %d, but is set to %d",
+ CblasUpper, CblasLower, Uplo);
+ if (TA != CblasNoTrans && TA != CblasTrans && TA != CblasConjTrans)
+ info = cblas_errprn(3, info,
+ "TransA must be %d, %d or %d, but is set to %d",
+ CblasNoTrans, CblasTrans, CblasConjTrans, TA);
+ if (Diag != CblasUnit && Diag != CblasNonUnit)
+ info = cblas_errprn(4, info, "DIAG must be %d or %d, but is set to %d",
+ CblasUnit, CblasNonUnit, Diag);
+
+ if (N < 0) info = cblas_errprn(5, info,
+ "N cannot be less than zero; is set to %d.", N);
+ if (lda < N || lda < 1)
+ info = cblas_errprn(7, info, "lda must be >= MAX(N,1): lda=%d N=%d",
+ lda, N);
+ if (!incX) info = cblas_errprn(9, info,
+ "incX cannot be zero; is set to %d.", incX);
+ if (info != 2000)
+ {
+ cblas_xerbla(info, "cblas_strsv", "");
+ return;
+ }
+#endif
+ if (incX < 0) x += (1-N)*incX;
+ if (Order == CblasColMajor)
+ ATL_strsv(Uplo, TA, Diag, N, A, lda, x, incX);
+ else
+ {
+ uplo = ( (Uplo == CblasUpper) ? CblasLower : CblasUpper );
+ if (TA == CblasNoTrans) ta = CblasTrans;
+ else ta = CblasNoTrans;
+ ATL_strsv(uplo, ta, Diag, N, A, lda, x, incX);
+ }
+}
diff --git a/scikits/learn/utils/arrayfuncs.c b/scikits/learn/utils/arrayfuncs.c
index 6e65702eda..068f615421 100644
--- a/scikits/learn/utils/arrayfuncs.c
+++ b/scikits/learn/utils/arrayfuncs.c
@@ -1,4 +1,4 @@
-/* Generated by Cython 0.12.1 on Fri Nov 5 15:03:01 2010 */
+/* Generated by Cython 0.12.1 on Mon Nov 22 14:49:19 2010 */
#define PY_SSIZE_T_CLEAN
#include "Python.h"
@@ -434,34 +434,6 @@ static void __Pyx_RaiseArgtupleInvalid(const char* func_name, int exact,
static int __Pyx_ParseOptionalKeywords(PyObject *kwds, PyObject **argnames[], PyObject *kwds2, PyObject *values[], Py_ssize_t num_pos_args, const char* function_name); /*proto*/
-/* Run-time type information about structs used with buffers */
-struct __Pyx_StructField_;
-
-typedef struct {
- const char* name; /* for error messages only */
- struct __Pyx_StructField_* fields;
- size_t size; /* sizeof(type) */
- char typegroup; /* _R_eal, _C_omplex, Signed _I_nt, _U_nsigned int, _S_truct, _P_ointer, _O_bject */
-} __Pyx_TypeInfo;
-
-typedef struct __Pyx_StructField_ {
- __Pyx_TypeInfo* type;
- const char* name;
- size_t offset;
-} __Pyx_StructField;
-
-typedef struct {
- __Pyx_StructField* field;
- size_t parent_offset;
-} __Pyx_BufFmt_StackElem;
-
-
-static CYTHON_INLINE void __Pyx_SafeReleaseBuffer(Py_buffer* info);
-static int __Pyx_GetBufferAndValidate(Py_buffer* buf, PyObject* obj, __Pyx_TypeInfo* dtype, int flags, int nd, int cast, __Pyx_BufFmt_StackElem* stack);
-
-static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb); /*proto*/
-static CYTHON_INLINE void __Pyx_ErrFetch(PyObject **type, PyObject **value, PyObject **tb); /*proto*/
-
static CYTHON_INLINE void __Pyx_RaiseNeedMoreValuesError(Py_ssize_t index);
static CYTHON_INLINE void __Pyx_RaiseTooManyValuesError(void);
@@ -477,21 +449,14 @@ static void __Pyx_UnpackTupleError(PyObject *, Py_ssize_t index); /*proto*/
static int __Pyx_ArgTypeTest(PyObject *obj, PyTypeObject *type, int none_allowed,
const char *name, int exact); /*proto*/
-#if PY_MAJOR_VERSION < 3
-static int __Pyx_GetBuffer(PyObject *obj, Py_buffer *view, int flags);
-static void __Pyx_ReleaseBuffer(Py_buffer *view);
-#else
-#define __Pyx_GetBuffer PyObject_GetBuffer
-#define __Pyx_ReleaseBuffer PyBuffer_Release
-#endif
-
-Py_ssize_t __Pyx_zeros[] = {0, 0};
-Py_ssize_t __Pyx_minusones[] = {-1, -1};
static PyObject *__Pyx_Import(PyObject *name, PyObject *from_list); /*proto*/
static PyObject *__Pyx_GetName(PyObject *dict, PyObject *name); /*proto*/
+static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb); /*proto*/
+static CYTHON_INLINE void __Pyx_ErrFetch(PyObject **type, PyObject **value, PyObject **tb); /*proto*/
+
static void __Pyx_Raise(PyObject *type, PyObject *value, PyObject *tb); /*proto*/
#if CYTHON_CCOMPLEX
@@ -643,7 +608,6 @@ static CYTHON_INLINE PyObject *__pyx_f_5numpy_get_array_base(PyArrayObject *); /
static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *, Py_ssize_t); /*proto*/
static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(double *, Py_ssize_t); /*proto*/
-static __Pyx_TypeInfo __Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE = { "scikits.learn.utils.arrayfuncs.DOUBLE", NULL, sizeof(__pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE), 'R' };
#define __Pyx_MODULE_NAME "scikits.learn.utils.arrayfuncs"
int __pyx_module_is_main_scikits__learn__utils__arrayfuncs = 0;
@@ -652,15 +616,16 @@ static PyObject *__pyx_builtin_ValueError;
static PyObject *__pyx_builtin_range;
static PyObject *__pyx_builtin_RuntimeError;
static char __pyx_k_1[] = "Unsupported dtype for array X";
-static char __pyx_k_2[] = "ndarray is not C contiguous";
-static char __pyx_k_3[] = "ndarray is not Fortran contiguous";
-static char __pyx_k_4[] = "Non-native byte order not supported";
-static char __pyx_k_5[] = "unknown dtype code in numpy.pxd (%d)";
-static char __pyx_k_6[] = "Format string allocated too short, see comment in numpy.pxd";
-static char __pyx_k_7[] = "Format string allocated too short.";
-static char __pyx_k_8[] = "\nSmall collection of auxiliary functions that operate on arrays\n\n";
-static char __pyx_k_9[] = "min_pos (line 42)";
-static char __pyx_k_10[] = "solve_triangular (line 73)";
+static char __pyx_k_2[] = "Unsupported or inconsistent dtype in arrays X, y";
+static char __pyx_k_3[] = "ndarray is not C contiguous";
+static char __pyx_k_4[] = "ndarray is not Fortran contiguous";
+static char __pyx_k_5[] = "Non-native byte order not supported";
+static char __pyx_k_6[] = "unknown dtype code in numpy.pxd (%d)";
+static char __pyx_k_7[] = "Format string allocated too short, see comment in numpy.pxd";
+static char __pyx_k_8[] = "Format string allocated too short.";
+static char __pyx_k_9[] = "\nSmall collection of auxiliary functions that operate on arrays\n\n";
+static char __pyx_k_10[] = "min_pos (line 47)";
+static char __pyx_k_11[] = "solve_triangular (line 78)";
static char __pyx_k__B[] = "B";
static char __pyx_k__H[] = "H";
static char __pyx_k__I[] = "I";
@@ -713,13 +678,14 @@ static char __pyx_k__RuntimeError[] = "RuntimeError";
static char __pyx_k__solve_triangular[] = "solve_triangular";
static PyObject *__pyx_kp_s_1;
static PyObject *__pyx_kp_u_10;
-static PyObject *__pyx_kp_u_2;
+static PyObject *__pyx_kp_u_11;
+static PyObject *__pyx_kp_s_2;
static PyObject *__pyx_kp_u_3;
static PyObject *__pyx_kp_u_4;
static PyObject *__pyx_kp_u_5;
static PyObject *__pyx_kp_u_6;
static PyObject *__pyx_kp_u_7;
-static PyObject *__pyx_kp_u_9;
+static PyObject *__pyx_kp_u_8;
static PyObject *__pyx_n_s__L;
static PyObject *__pyx_n_s__RuntimeError;
static PyObject *__pyx_n_s__ValueError;
@@ -756,7 +722,7 @@ static PyObject *__pyx_n_s__type_num;
static PyObject *__pyx_n_s__y;
static PyObject *__pyx_int_15;
-/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":42
+/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":47
*
*
* def min_pos(np.ndarray X): # <<<<<<<<<<<<<<
@@ -765,7 +731,7 @@ static PyObject *__pyx_int_15;
*/
static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *__pyx_self, PyObject *__pyx_v_X); /*proto*/
-static char __pyx_doc_7scikits_5learn_5utils_10arrayfuncs_min_pos[] = "\n Find the minimum value of an array over positivie values\n ";
+static char __pyx_doc_7scikits_5learn_5utils_10arrayfuncs_min_pos[] = "\n Find the minimum value of an array over positivie values\n\n Returns a huge value if none of the values are positive\n ";
static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *__pyx_self, PyObject *__pyx_v_X) {
PyObject *__pyx_r = NULL;
PyObject *__pyx_t_1 = NULL;
@@ -775,28 +741,28 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
__Pyx_RefNannySetupContext("min_pos");
__pyx_self = __pyx_self;
__Pyx_INCREF((PyObject *)__pyx_v_X);
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 42; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":46
- * Find the minimum value of an array over positivie values
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":53
+ * Returns a huge value if none of the values are positive
* """
* if X.dtype.name == 'float32': # <<<<<<<<<<<<<<
* return _float_min_pos(<float *> X.data, X.size)
* elif X.dtype.name == 'float64':
*/
- __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 46; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 53; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
- __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 46; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 53; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float32), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 46; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float32), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 53; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
- __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 46; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 53; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
if (__pyx_t_3) {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":47
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":54
* """
* if X.dtype.name == 'float32':
* return _float_min_pos(<float *> X.data, X.size) # <<<<<<<<<<<<<<
@@ -804,11 +770,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
* return _double_min_pos(<double *> X.data, X.size)
*/
__Pyx_XDECREF(__pyx_r);
- __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__size); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__size); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 54; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
- __pyx_t_4 = __Pyx_PyIndex_AsSsize_t(__pyx_t_1); if (unlikely((__pyx_t_4 == (Py_ssize_t)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_4 = __Pyx_PyIndex_AsSsize_t(__pyx_t_1); if (unlikely((__pyx_t_4 == (Py_ssize_t)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 54; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __pyx_t_1 = PyFloat_FromDouble(__pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(((float *)((PyArrayObject *)__pyx_v_X)->data), __pyx_t_4)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyFloat_FromDouble(__pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(((float *)((PyArrayObject *)__pyx_v_X)->data), __pyx_t_4)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 54; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_r = __pyx_t_1;
__pyx_t_1 = 0;
@@ -816,26 +782,26 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
goto __pyx_L5;
}
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":48
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":55
* if X.dtype.name == 'float32':
* return _float_min_pos(<float *> X.data, X.size)
* elif X.dtype.name == 'float64': # <<<<<<<<<<<<<<
* return _double_min_pos(<double *> X.data, X.size)
* else:
*/
- __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 55; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
- __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 55; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float64), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float64), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 55; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
- __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 55; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
if (__pyx_t_3) {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":49
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":56
* return _float_min_pos(<float *> X.data, X.size)
* elif X.dtype.name == 'float64':
* return _double_min_pos(<double *> X.data, X.size) # <<<<<<<<<<<<<<
@@ -843,11 +809,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
* raise ValueError('Unsupported dtype for array X')
*/
__Pyx_XDECREF(__pyx_r);
- __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__size); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyObject_GetAttr(__pyx_v_X, __pyx_n_s__size); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 56; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
- __pyx_t_4 = __Pyx_PyIndex_AsSsize_t(__pyx_t_1); if (unlikely((__pyx_t_4 == (Py_ssize_t)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_4 = __Pyx_PyIndex_AsSsize_t(__pyx_t_1); if (unlikely((__pyx_t_4 == (Py_ssize_t)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 56; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __pyx_t_1 = PyFloat_FromDouble(__pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(((double *)((PyArrayObject *)__pyx_v_X)->data), __pyx_t_4)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyFloat_FromDouble(__pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(((double *)((PyArrayObject *)__pyx_v_X)->data), __pyx_t_4)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 56; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_r = __pyx_t_1;
__pyx_t_1 = 0;
@@ -856,24 +822,24 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
}
/*else*/ {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":51
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":58
* return _double_min_pos(<double *> X.data, X.size)
* else:
* raise ValueError('Unsupported dtype for array X') # <<<<<<<<<<<<<<
*
*
*/
- __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 51; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_INCREF(((PyObject *)__pyx_kp_s_1));
PyTuple_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_kp_s_1));
__Pyx_GIVEREF(((PyObject *)__pyx_kp_s_1));
- __pyx_t_2 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_1, NULL); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 51; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_2 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_1, NULL); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_Raise(__pyx_t_2, 0, 0);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 51; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_L5:;
@@ -891,7 +857,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_min_pos(PyObject *
return __pyx_r;
}
-/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":54
+/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":61
*
*
* cdef float _float_min_pos (float *X, Py_ssize_t size): # <<<<<<<<<<<<<<
@@ -910,7 +876,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
int __pyx_t_5;
__Pyx_RefNannySetupContext("_float_min_pos");
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":56
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":63
* cdef float _float_min_pos (float *X, Py_ssize_t size):
* cdef Py_ssize_t i
* cdef float min_val = DBL_MAX # <<<<<<<<<<<<<<
@@ -919,7 +885,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
*/
__pyx_v_min_val = DBL_MAX;
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":57
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":64
* cdef Py_ssize_t i
* cdef float min_val = DBL_MAX
* for i in range(size): # <<<<<<<<<<<<<<
@@ -930,7 +896,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
__pyx_v_i = __pyx_t_2;
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":58
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":65
* cdef float min_val = DBL_MAX
* for i in range(size):
* if X[i] > 0. and X[i] < min_val: # <<<<<<<<<<<<<<
@@ -946,7 +912,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
}
if (__pyx_t_5) {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":59
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":66
* for i in range(size):
* if X[i] > 0. and X[i] < min_val:
* min_val = X[i] # <<<<<<<<<<<<<<
@@ -959,7 +925,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
__pyx_L5:;
}
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":60
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":67
* if X[i] > 0. and X[i] < min_val:
* min_val = X[i]
* return min_val # <<<<<<<<<<<<<<
@@ -975,7 +941,7 @@ static float __pyx_f_7scikits_5learn_5utils_10arrayfuncs__float_min_pos(float *
return __pyx_r;
}
-/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":63
+/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":70
*
*
* cdef double _double_min_pos (double *X, Py_ssize_t size): # <<<<<<<<<<<<<<
@@ -994,7 +960,7 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
int __pyx_t_5;
__Pyx_RefNannySetupContext("_double_min_pos");
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":65
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":72
* cdef double _double_min_pos (double *X, Py_ssize_t size):
* cdef Py_ssize_t i
* cdef np.float64_t min_val = FLT_MAX # <<<<<<<<<<<<<<
@@ -1003,7 +969,7 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
*/
__pyx_v_min_val = FLT_MAX;
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":66
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":73
* cdef Py_ssize_t i
* cdef np.float64_t min_val = FLT_MAX
* for i in range(size): # <<<<<<<<<<<<<<
@@ -1014,7 +980,7 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
__pyx_v_i = __pyx_t_2;
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":67
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":74
* cdef np.float64_t min_val = FLT_MAX
* for i in range(size):
* if X[i] > 0. and X[i] < min_val: # <<<<<<<<<<<<<<
@@ -1030,7 +996,7 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
}
if (__pyx_t_5) {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":68
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":75
* for i in range(size):
* if X[i] > 0. and X[i] < min_val:
* min_val = X[i] # <<<<<<<<<<<<<<
@@ -1043,12 +1009,12 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
__pyx_L5:;
}
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":69
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":76
* if X[i] > 0. and X[i] < min_val:
* min_val = X[i]
* return min_val # <<<<<<<<<<<<<<
*
- * @cython.boundscheck(False)
+ * def solve_triangular (np.ndarray X, np.ndarray y):
*/
__pyx_r = __pyx_v_min_val;
goto __pyx_L0;
@@ -1059,31 +1025,28 @@ static double __pyx_f_7scikits_5learn_5utils_10arrayfuncs__double_min_pos(doubl
return __pyx_r;
}
-/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":73
- * @cython.boundscheck(False)
- * @cython.wraparound(False)
- * def solve_triangular ( # <<<<<<<<<<<<<<
- * np.ndarray[DOUBLE, ndim=2] X,
- * np.ndarray[DOUBLE, ndim=1] y
+/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":78
+ * return min_val
+ *
+ * def solve_triangular (np.ndarray X, np.ndarray y): # <<<<<<<<<<<<<<
+ * """
+ * Solves a triangular system (overwrites y)
*/
static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static char __pyx_doc_7scikits_5learn_5utils_10arrayfuncs_solve_triangular[] = "\n Solves a triangular system (overwrites y)\n\n TODO: option for upper, lower, etc.\n ";
+static char __pyx_doc_7scikits_5learn_5utils_10arrayfuncs_solve_triangular[] = "\n Solves a triangular system (overwrites y)\n\n Note: The lapack function to solve triangular systems was added to\n scipy v0.9. Remove this when we stop supporting earlier versions.\n ";
static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
PyArrayObject *__pyx_v_X = 0;
PyArrayObject *__pyx_v_y = 0;
int __pyx_v_lda;
- Py_buffer __pyx_bstruct_y;
- Py_ssize_t __pyx_bstride_0_y = 0;
- Py_ssize_t __pyx_bshape_0_y = 0;
- Py_buffer __pyx_bstruct_X;
- Py_ssize_t __pyx_bstride_0_X = 0;
- Py_ssize_t __pyx_bstride_1_X = 0;
- Py_ssize_t __pyx_bshape_0_X = 0;
- Py_ssize_t __pyx_bshape_1_X = 0;
PyObject *__pyx_r = NULL;
- int __pyx_t_1;
- size_t __pyx_t_2;
+ PyObject *__pyx_t_1 = NULL;
+ PyObject *__pyx_t_2 = NULL;
+ int __pyx_t_3;
+ int __pyx_t_4;
+ int __pyx_t_5;
+ int __pyx_t_6;
+ size_t __pyx_t_7;
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s__X,&__pyx_n_s__y,0};
__Pyx_RefNannySetupContext("solve_triangular");
__pyx_self = __pyx_self;
@@ -1105,11 +1068,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__y);
if (likely(values[1])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 73; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 78; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
}
if (unlikely(kw_args > 0)) {
- if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "solve_triangular") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 73; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "solve_triangular") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 78; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
__pyx_v_X = ((PyArrayObject *)values[0]);
__pyx_v_y = ((PyArrayObject *)values[1]);
@@ -1121,78 +1084,179 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
- __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 73; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 78; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.solve_triangular");
return NULL;
__pyx_L4_argument_unpacking_done:;
- __pyx_bstruct_X.buf = NULL;
- __pyx_bstruct_y.buf = NULL;
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 74; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 75; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- {
- __Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_X, (PyObject*)__pyx_v_X, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 73; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_INCREF((PyObject *)__pyx_v_X);
+ __Pyx_INCREF((PyObject *)__pyx_v_y);
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 78; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 78; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":87
+ * cdef int lda
+ *
+ * if X.dtype.name == 'float64' and y.dtype.name == 'float64': # <<<<<<<<<<<<<<
+ * lda = <int> X.strides[0] / sizeof(double)
+ *
+ */
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_X), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float64), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ if (__pyx_t_3) {
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_y), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float64), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_4 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_4 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_5 = __pyx_t_4;
+ } else {
+ __pyx_t_5 = __pyx_t_3;
}
- __pyx_bstride_0_X = __pyx_bstruct_X.strides[0]; __pyx_bstride_1_X = __pyx_bstruct_X.strides[1];
- __pyx_bshape_0_X = __pyx_bstruct_X.shape[0]; __pyx_bshape_1_X = __pyx_bstruct_X.shape[1];
- {
- __Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_y, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 73; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (__pyx_t_5) {
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":88
+ *
+ * if X.dtype.name == 'float64' and y.dtype.name == 'float64':
+ * lda = <int> X.strides[0] / sizeof(double) # <<<<<<<<<<<<<<
+ *
+ * cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ */
+ __pyx_t_6 = ((int)(__pyx_v_X->strides[0]));
+ __pyx_t_7 = (sizeof(double));
+ if (unlikely(__pyx_t_7 == 0)) {
+ PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 88; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_v_lda = (__pyx_t_6 / __pyx_t_7);
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":92
+ * cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ * CblasNonUnit, <int> X.shape[0], <double *> X.data,
+ * lda, <double *> y.data, 1); # <<<<<<<<<<<<<<
+ *
+ * elif X.dtype.name == 'float32' and y.dtype.name == 'float32':
+ */
+ cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, ((int)(__pyx_v_X->dimensions[0])), ((double *)__pyx_v_X->data), __pyx_v_lda, ((double *)__pyx_v_y->data), 1);
+ goto __pyx_L6;
}
- __pyx_bstride_0_y = __pyx_bstruct_y.strides[0];
- __pyx_bshape_0_y = __pyx_bstruct_y.shape[0];
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":82
- * TODO: option for upper, lower, etc.
- * """
- * cdef int lda = <int> X.strides[0] / sizeof(double) # <<<<<<<<<<<<<<
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":94
+ * lda, <double *> y.data, 1);
+ *
+ * elif X.dtype.name == 'float32' and y.dtype.name == 'float32': # <<<<<<<<<<<<<<
+ * lda = <int> X.strides[0] / sizeof(float)
+ *
+ */
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_X), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float32), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_5 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_5 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ if (__pyx_t_5) {
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_y), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float32), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_4 = __pyx_t_3;
+ } else {
+ __pyx_t_4 = __pyx_t_5;
+ }
+ if (__pyx_t_4) {
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":95
+ *
+ * elif X.dtype.name == 'float32' and y.dtype.name == 'float32':
+ * lda = <int> X.strides[0] / sizeof(float) # <<<<<<<<<<<<<<
*
- * cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ * cblas_strsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ */
+ __pyx_t_6 = ((int)(__pyx_v_X->strides[0]));
+ __pyx_t_7 = (sizeof(float));
+ if (unlikely(__pyx_t_7 == 0)) {
+ PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 95; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_v_lda = (__pyx_t_6 / __pyx_t_7);
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":99
+ * cblas_strsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ * CblasNonUnit, <int> X.shape[0], <float *> X.data,
+ * lda, <float *> y.data, 1); # <<<<<<<<<<<<<<
+ * else:
+ * raise ValueError ('Unsupported or inconsistent dtype in arrays X, y')
*/
- __pyx_t_1 = ((int)(__pyx_v_X->strides[0]));
- __pyx_t_2 = (sizeof(double));
- if (unlikely(__pyx_t_2 == 0)) {
- PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ cblas_strsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, ((int)(__pyx_v_X->dimensions[0])), ((float *)__pyx_v_X->data), __pyx_v_lda, ((float *)__pyx_v_y->data), 1);
+ goto __pyx_L6;
}
- __pyx_v_lda = (__pyx_t_1 / __pyx_t_2);
+ /*else*/ {
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":86
- * cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
- * CblasNonUnit, <int> X.shape[0], <double *> X.data,
- * lda, <double *> y.data, 1); # <<<<<<<<<<<<<<
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":101
+ * lda, <float *> y.data, 1);
+ * else:
+ * raise ValueError ('Unsupported or inconsistent dtype in arrays X, y') # <<<<<<<<<<<<<<
*
*
*/
- cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, ((int)(__pyx_v_X->dimensions[0])), ((double *)__pyx_v_X->data), __pyx_v_lda, ((double *)__pyx_v_y->data), 1);
+ __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 101; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_INCREF(((PyObject *)__pyx_kp_s_2));
+ PyTuple_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_kp_s_2));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_s_2));
+ __pyx_t_2 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_1, NULL); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 101; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __Pyx_Raise(__pyx_t_2, 0, 0);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 101; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_L6:;
__pyx_r = Py_None; __Pyx_INCREF(Py_None);
goto __pyx_L0;
__pyx_L1_error:;
- { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
- __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_y);
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_X);
- __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
+ __Pyx_XDECREF(__pyx_t_1);
+ __Pyx_XDECREF(__pyx_t_2);
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.solve_triangular");
__pyx_r = NULL;
- goto __pyx_L2;
__pyx_L0:;
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_y);
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_X);
- __pyx_L2:;
+ __Pyx_DECREF((PyObject *)__pyx_v_X);
+ __Pyx_DECREF((PyObject *)__pyx_v_y);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
-/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":89
+/* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":104
+ *
*
+ * def cholesky_delete (np.ndarray L, int go_out): # <<<<<<<<<<<<<<
*
- * def cholesky_delete ( # <<<<<<<<<<<<<<
- * np.ndarray[DOUBLE, ndim=2] L,
- * int go_out):
+ * cdef int n = <int> L.shape[0]
*/
static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
@@ -1201,14 +1265,12 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete(Py
int __pyx_v_go_out;
int __pyx_v_n;
int __pyx_v_m;
- Py_buffer __pyx_bstruct_L;
- Py_ssize_t __pyx_bstride_0_L = 0;
- Py_ssize_t __pyx_bstride_1_L = 0;
- Py_ssize_t __pyx_bshape_0_L = 0;
- Py_ssize_t __pyx_bshape_1_L = 0;
PyObject *__pyx_r = NULL;
- int __pyx_t_1;
- size_t __pyx_t_2;
+ PyObject *__pyx_t_1 = NULL;
+ PyObject *__pyx_t_2 = NULL;
+ int __pyx_t_3;
+ int __pyx_t_4;
+ size_t __pyx_t_5;
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s__L,&__pyx_n_s__go_out,0};
__Pyx_RefNannySetupContext("cholesky_delete");
__pyx_self = __pyx_self;
@@ -1230,79 +1292,138 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete(Py
values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__go_out);
if (likely(values[1])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 89; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
}
if (unlikely(kw_args > 0)) {
- if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "cholesky_delete") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 89; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "cholesky_delete") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
__pyx_v_L = ((PyArrayObject *)values[0]);
- __pyx_v_go_out = __Pyx_PyInt_AsInt(values[1]); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __pyx_v_go_out = __Pyx_PyInt_AsInt(values[1]); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
} else if (PyTuple_GET_SIZE(__pyx_args) != 2) {
goto __pyx_L5_argtuple_error;
} else {
__pyx_v_L = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 0));
- __pyx_v_go_out = __Pyx_PyInt_AsInt(PyTuple_GET_ITEM(__pyx_args, 1)); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __pyx_v_go_out = __Pyx_PyInt_AsInt(PyTuple_GET_ITEM(__pyx_args, 1)); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
- __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 89; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.cholesky_delete");
return NULL;
__pyx_L4_argument_unpacking_done:;
- __pyx_bstruct_L.buf = NULL;
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_L), __pyx_ptype_5numpy_ndarray, 1, "L", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 90; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- {
- __Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_L, (PyObject*)__pyx_v_L, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 89; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- }
- __pyx_bstride_0_L = __pyx_bstruct_L.strides[0]; __pyx_bstride_1_L = __pyx_bstruct_L.strides[1];
- __pyx_bshape_0_L = __pyx_bstruct_L.shape[0]; __pyx_bshape_1_L = __pyx_bstruct_L.shape[1];
+ __Pyx_INCREF((PyObject *)__pyx_v_L);
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_L), __pyx_ptype_5numpy_ndarray, 1, "L", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 104; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":93
- * int go_out):
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":106
+ * def cholesky_delete (np.ndarray L, int go_out):
*
* cdef int n = <int> L.shape[0] # <<<<<<<<<<<<<<
- * cdef int m = <int> L.strides[0] / sizeof (double)
- * c_cholesky_delete (m, n, <double *> L.data, go_out)
+ * cdef int m
+ *
*/
__pyx_v_n = ((int)(__pyx_v_L->dimensions[0]));
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":94
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":109
+ * cdef int m
*
- * cdef int n = <int> L.shape[0]
- * cdef int m = <int> L.strides[0] / sizeof (double) # <<<<<<<<<<<<<<
- * c_cholesky_delete (m, n, <double *> L.data, go_out)
- */
- __pyx_t_1 = ((int)(__pyx_v_L->strides[0]));
- __pyx_t_2 = (sizeof(double));
- if (unlikely(__pyx_t_2 == 0)) {
- PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ * if L.dtype.name == 'float64': # <<<<<<<<<<<<<<
+ * m = <int> L.strides[0] / sizeof (double)
+ * double_cholesky_delete (m, n, <double *> L.data, go_out)
+ */
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_L), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 109; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 109; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float64), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 109; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 109; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ if (__pyx_t_3) {
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":110
+ *
+ * if L.dtype.name == 'float64':
+ * m = <int> L.strides[0] / sizeof (double) # <<<<<<<<<<<<<<
+ * double_cholesky_delete (m, n, <double *> L.data, go_out)
+ * elif L.dtype.name == 'float32':
+ */
+ __pyx_t_4 = ((int)(__pyx_v_L->strides[0]));
+ __pyx_t_5 = (sizeof(double));
+ if (unlikely(__pyx_t_5 == 0)) {
+ PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_v_m = (__pyx_t_4 / __pyx_t_5);
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":111
+ * if L.dtype.name == 'float64':
+ * m = <int> L.strides[0] / sizeof (double)
+ * double_cholesky_delete (m, n, <double *> L.data, go_out) # <<<<<<<<<<<<<<
+ * elif L.dtype.name == 'float32':
+ * m = <int> L.strides[0] / sizeof (float)
+ */
+ double_cholesky_delete(__pyx_v_m, __pyx_v_n, ((double *)__pyx_v_L->data), __pyx_v_go_out);
+ goto __pyx_L6;
}
- __pyx_v_m = (__pyx_t_1 / __pyx_t_2);
- /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":95
- * cdef int n = <int> L.shape[0]
- * cdef int m = <int> L.strides[0] / sizeof (double)
- * c_cholesky_delete (m, n, <double *> L.data, go_out) # <<<<<<<<<<<<<<
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":112
+ * m = <int> L.strides[0] / sizeof (double)
+ * double_cholesky_delete (m, n, <double *> L.data, go_out)
+ * elif L.dtype.name == 'float32': # <<<<<<<<<<<<<<
+ * m = <int> L.strides[0] / sizeof (float)
+ * float_cholesky_delete (m, n, <float *> L.data, go_out)
+ */
+ __pyx_t_1 = PyObject_GetAttr(((PyObject *)__pyx_v_L), __pyx_n_s__dtype); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__name); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_RichCompare(__pyx_t_2, ((PyObject *)__pyx_n_s__float32), Py_EQ); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_3 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_3 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ if (__pyx_t_3) {
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":113
+ * double_cholesky_delete (m, n, <double *> L.data, go_out)
+ * elif L.dtype.name == 'float32':
+ * m = <int> L.strides[0] / sizeof (float) # <<<<<<<<<<<<<<
+ * float_cholesky_delete (m, n, <float *> L.data, go_out)
+ *
*/
- cholesky_delete(__pyx_v_m, __pyx_v_n, ((double *)__pyx_v_L->data), __pyx_v_go_out);
+ __pyx_t_4 = ((int)(__pyx_v_L->strides[0]));
+ __pyx_t_5 = (sizeof(float));
+ if (unlikely(__pyx_t_5 == 0)) {
+ PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 113; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_v_m = (__pyx_t_4 / __pyx_t_5);
+
+ /* "/home/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":114
+ * elif L.dtype.name == 'float32':
+ * m = <int> L.strides[0] / sizeof (float)
+ * float_cholesky_delete (m, n, <float *> L.data, go_out) # <<<<<<<<<<<<<<
+ *
+ */
+ float_cholesky_delete(__pyx_v_m, __pyx_v_n, ((float *)__pyx_v_L->data), __pyx_v_go_out);
+ goto __pyx_L6;
+ }
+ __pyx_L6:;
__pyx_r = Py_None; __Pyx_INCREF(Py_None);
goto __pyx_L0;
__pyx_L1_error:;
- { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
- __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_L);
- __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
+ __Pyx_XDECREF(__pyx_t_1);
+ __Pyx_XDECREF(__pyx_t_2);
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.cholesky_delete");
__pyx_r = NULL;
- goto __pyx_L2;
__pyx_L0:;
- __Pyx_SafeReleaseBuffer(&__pyx_bstruct_L);
- __pyx_L2:;
+ __Pyx_DECREF((PyObject *)__pyx_v_L);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
@@ -1437,9 +1558,9 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyObject *__pyx_v_self, Py_buf
*/
__pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 205; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_4);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_2));
- PyTuple_SET_ITEM(__pyx_t_4, 0, ((PyObject *)__pyx_kp_u_2));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_2));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_3));
+ PyTuple_SET_ITEM(__pyx_t_4, 0, ((PyObject *)__pyx_kp_u_3));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_3));
__pyx_t_5 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_4, NULL); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 205; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
@@ -1483,9 +1604,9 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyObject *__pyx_v_self, Py_buf
*/
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 209; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_3));
- PyTuple_SET_ITEM(__pyx_t_5, 0, ((PyObject *)__pyx_kp_u_3));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_3));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_4));
+ PyTuple_SET_ITEM(__pyx_t_5, 0, ((PyObject *)__pyx_kp_u_4));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_4));
__pyx_t_4 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_5, NULL); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 209; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
@@ -1760,9 +1881,9 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyObject *__pyx_v_self, Py_buf
*/
__pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 247; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_4);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_4));
- PyTuple_SET_ITEM(__pyx_t_4, 0, ((PyObject *)__pyx_kp_u_4));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_4));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_5));
+ PyTuple_SET_ITEM(__pyx_t_4, 0, ((PyObject *)__pyx_kp_u_5));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_5));
__pyx_t_5 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_4, NULL); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 247; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
@@ -2004,7 +2125,7 @@ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyObject *__pyx_v_self, Py_buf
*/
__pyx_t_5 = PyInt_FromLong(__pyx_v_t); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 266; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
- __pyx_t_4 = PyNumber_Remainder(((PyObject *)__pyx_kp_u_5), __pyx_t_5); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 266; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_4 = PyNumber_Remainder(((PyObject *)__pyx_kp_u_6), __pyx_t_5); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 266; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 266; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -2513,9 +2634,9 @@ static CYTHON_INLINE char *__pyx_f_5numpy__util_dtypestring(PyArray_Descr *__pyx
*/
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 786; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_6));
- PyTuple_SET_ITEM(__pyx_t_5, 0, ((PyObject *)__pyx_kp_u_6));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_6));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_7));
+ PyTuple_SET_ITEM(__pyx_t_5, 0, ((PyObject *)__pyx_kp_u_7));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_7));
__pyx_t_3 = PyObject_Call(__pyx_builtin_RuntimeError, __pyx_t_5, NULL); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 786; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
@@ -2570,9 +2691,9 @@ static CYTHON_INLINE char *__pyx_f_5numpy__util_dtypestring(PyArray_Descr *__pyx
*/
__pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 790; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_3);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_4));
- PyTuple_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_kp_u_4));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_4));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_5));
+ PyTuple_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_kp_u_5));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_5));
__pyx_t_5 = PyObject_Call(__pyx_builtin_ValueError, __pyx_t_3, NULL); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 790; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
@@ -2679,9 +2800,9 @@ static CYTHON_INLINE char *__pyx_f_5numpy__util_dtypestring(PyArray_Descr *__pyx
*/
__pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 810; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_3);
- __Pyx_INCREF(((PyObject *)__pyx_kp_u_7));
- PyTuple_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_kp_u_7));
- __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_7));
+ __Pyx_INCREF(((PyObject *)__pyx_kp_u_8));
+ PyTuple_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_kp_u_8));
+ __Pyx_GIVEREF(((PyObject *)__pyx_kp_u_8));
__pyx_t_5 = PyObject_Call(__pyx_builtin_RuntimeError, __pyx_t_3, NULL); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 810; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
@@ -3029,7 +3150,7 @@ static CYTHON_INLINE char *__pyx_f_5numpy__util_dtypestring(PyArray_Descr *__pyx
* f += 1
* else:
*/
- __pyx_t_3 = PyNumber_Remainder(((PyObject *)__pyx_kp_u_5), __pyx_v_t); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 831; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_3 = PyNumber_Remainder(((PyObject *)__pyx_kp_u_6), __pyx_v_t); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 831; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 831; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_5);
@@ -3253,7 +3374,7 @@ static void __pyx_init_filenames(void); /*proto*/
static struct PyModuleDef __pyx_moduledef = {
PyModuleDef_HEAD_INIT,
__Pyx_NAMESTR("arrayfuncs"),
- __Pyx_DOCSTR(__pyx_k_8), /* m_doc */
+ __Pyx_DOCSTR(__pyx_k_9), /* m_doc */
-1, /* m_size */
__pyx_methods /* m_methods */,
NULL, /* m_reload */
@@ -3266,13 +3387,14 @@ static struct PyModuleDef __pyx_moduledef = {
static __Pyx_StringTabEntry __pyx_string_tab[] = {
{&__pyx_kp_s_1, __pyx_k_1, sizeof(__pyx_k_1), 0, 0, 1, 0},
{&__pyx_kp_u_10, __pyx_k_10, sizeof(__pyx_k_10), 0, 1, 0, 0},
- {&__pyx_kp_u_2, __pyx_k_2, sizeof(__pyx_k_2), 0, 1, 0, 0},
+ {&__pyx_kp_u_11, __pyx_k_11, sizeof(__pyx_k_11), 0, 1, 0, 0},
+ {&__pyx_kp_s_2, __pyx_k_2, sizeof(__pyx_k_2), 0, 0, 1, 0},
{&__pyx_kp_u_3, __pyx_k_3, sizeof(__pyx_k_3), 0, 1, 0, 0},
{&__pyx_kp_u_4, __pyx_k_4, sizeof(__pyx_k_4), 0, 1, 0, 0},
{&__pyx_kp_u_5, __pyx_k_5, sizeof(__pyx_k_5), 0, 1, 0, 0},
{&__pyx_kp_u_6, __pyx_k_6, sizeof(__pyx_k_6), 0, 1, 0, 0},
{&__pyx_kp_u_7, __pyx_k_7, sizeof(__pyx_k_7), 0, 1, 0, 0},
- {&__pyx_kp_u_9, __pyx_k_9, sizeof(__pyx_k_9), 0, 1, 0, 0},
+ {&__pyx_kp_u_8, __pyx_k_8, sizeof(__pyx_k_8), 0, 1, 0, 0},
{&__pyx_n_s__L, __pyx_k__L, sizeof(__pyx_k__L), 0, 0, 1, 1},
{&__pyx_n_s__RuntimeError, __pyx_k__RuntimeError, sizeof(__pyx_k__RuntimeError), 0, 0, 1, 1},
{&__pyx_n_s__ValueError, __pyx_k__ValueError, sizeof(__pyx_k__ValueError), 0, 0, 1, 1},
@@ -3310,8 +3432,8 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
{0, 0, 0, 0, 0, 0, 0}
};
static int __Pyx_InitCachedBuiltins(void) {
- __pyx_builtin_ValueError = __Pyx_GetName(__pyx_b, __pyx_n_s__ValueError); if (!__pyx_builtin_ValueError) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 51; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __pyx_builtin_range = __Pyx_GetName(__pyx_b, __pyx_n_s__range); if (!__pyx_builtin_range) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_builtin_ValueError = __Pyx_GetName(__pyx_b, __pyx_n_s__ValueError); if (!__pyx_builtin_ValueError) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_builtin_range = __Pyx_GetName(__pyx_b, __pyx_n_s__range); if (!__pyx_builtin_range) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 64; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_builtin_RuntimeError = __Pyx_GetName(__pyx_b, __pyx_n_s__RuntimeError); if (!__pyx_builtin_RuntimeError) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 786; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
return 0;
__pyx_L1_error:;
@@ -3364,7 +3486,7 @@ PyMODINIT_FUNC PyInit_arrayfuncs(void)
#endif
/*--- Module creation code ---*/
#if PY_MAJOR_VERSION < 3
- __pyx_m = Py_InitModule4(__Pyx_NAMESTR("arrayfuncs"), __pyx_methods, __Pyx_DOCSTR(__pyx_k_8), 0, PYTHON_API_VERSION);
+ __pyx_m = Py_InitModule4(__Pyx_NAMESTR("arrayfuncs"), __pyx_methods, __Pyx_DOCSTR(__pyx_k_9), 0, PYTHON_API_VERSION);
#else
__pyx_m = PyModule_Create(&__pyx_moduledef);
#endif
@@ -3418,14 +3540,14 @@ PyMODINIT_FUNC PyInit_arrayfuncs(void)
__pyx_t_3 = __Pyx_GetAttrString(__pyx_t_2, "__doc__");
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
- if (PyDict_SetItem(__pyx_t_1, ((PyObject *)__pyx_kp_u_9), __pyx_t_3) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (PyDict_SetItem(__pyx_t_1, ((PyObject *)__pyx_kp_u_10), __pyx_t_3) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = PyObject_GetAttr(__pyx_m, __pyx_n_s__solve_triangular); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = __Pyx_GetAttrString(__pyx_t_3, "__doc__");
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
- if (PyDict_SetItem(__pyx_t_1, ((PyObject *)__pyx_kp_u_10), __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (PyDict_SetItem(__pyx_t_1, ((PyObject *)__pyx_kp_u_11), __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
if (PyObject_SetAttr(__pyx_m, __pyx_n_s____test__, ((PyObject *)__pyx_t_1)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
@@ -3589,477 +3711,6 @@ bad:
return -1;
}
-static CYTHON_INLINE int __Pyx_IsLittleEndian(void) {
- unsigned int n = 1;
- return *(unsigned char*)(&n) != 0;
-}
-
-typedef struct {
- __Pyx_StructField root;
- __Pyx_BufFmt_StackElem* head;
- size_t fmt_offset;
- int new_count, enc_count;
- int is_complex;
- char enc_type;
- char packmode;
-} __Pyx_BufFmt_Context;
-
-static void __Pyx_BufFmt_Init(__Pyx_BufFmt_Context* ctx,
- __Pyx_BufFmt_StackElem* stack,
- __Pyx_TypeInfo* type) {
- stack[0].field = &ctx->root;
- stack[0].parent_offset = 0;
- ctx->root.type = type;
- ctx->root.name = "buffer dtype";
- ctx->root.offset = 0;
- ctx->head = stack;
- ctx->head->field = &ctx->root;
- ctx->fmt_offset = 0;
- ctx->head->parent_offset = 0;
- ctx->packmode = '@';
- ctx->new_count = 1;
- ctx->enc_count = 0;
- ctx->enc_type = 0;
- ctx->is_complex = 0;
- while (type->typegroup == 'S') {
- ++ctx->head;
- ctx->head->field = type->fields;
- ctx->head->parent_offset = 0;
- type = type->fields->type;
- }
-}
-
-static int __Pyx_BufFmt_ParseNumber(const char** ts) {
- int count;
- const char* t = *ts;
- if (*t < '0' || *t > '9') {
- return -1;
- } else {
- count = *t++ - '0';
- while (*t >= '0' && *t < '9') {
- count *= 10;
- count += *t++ - '0';
- }
- }
- *ts = t;
- return count;
-}
-
-static void __Pyx_BufFmt_RaiseUnexpectedChar(char ch) {
- char msg[] = {ch, 0};
- PyErr_Format(PyExc_ValueError, "Unexpected format string character: '%s'", msg);
-}
-
-static const char* __Pyx_BufFmt_DescribeTypeChar(char ch, int is_complex) {
- switch (ch) {
- case 'b': return "'char'";
- case 'B': return "'unsigned char'";
- case 'h': return "'short'";
- case 'H': return "'unsigned short'";
- case 'i': return "'int'";
- case 'I': return "'unsigned int'";
- case 'l': return "'long'";
- case 'L': return "'unsigned long'";
- case 'q': return "'long long'";
- case 'Q': return "'unsigned long long'";
- case 'f': return (is_complex ? "'complex float'" : "'float'");
- case 'd': return (is_complex ? "'complex double'" : "'double'");
- case 'g': return (is_complex ? "'complex long double'" : "'long double'");
- case 'T': return "a struct";
- case 'O': return "Python object";
- case 'P': return "a pointer";
- case 0: return "end";
- default: return "unparseable format string";
- }
-}
-
-static size_t __Pyx_BufFmt_TypeCharToStandardSize(char ch, int is_complex) {
- switch (ch) {
- case '?': case 'c': case 'b': case 'B': return 1;
- case 'h': case 'H': return 2;
- case 'i': case 'I': case 'l': case 'L': return 4;
- case 'q': case 'Q': return 8;
- case 'f': return (is_complex ? 8 : 4);
- case 'd': return (is_complex ? 16 : 8);
- case 'g': {
- PyErr_SetString(PyExc_ValueError, "Python does not define a standard format string size for long double ('g')..");
- return 0;
- }
- case 'O': case 'P': return sizeof(void*);
- default:
- __Pyx_BufFmt_RaiseUnexpectedChar(ch);
- return 0;
- }
-}
-
-static size_t __Pyx_BufFmt_TypeCharToNativeSize(char ch, int is_complex) {
- switch (ch) {
- case 'c': case 'b': case 'B': return 1;
- case 'h': case 'H': return sizeof(short);
- case 'i': case 'I': return sizeof(int);
- case 'l': case 'L': return sizeof(long);
- #ifdef HAVE_LONG_LONG
- case 'q': case 'Q': return sizeof(PY_LONG_LONG);
- #endif
- case 'f': return sizeof(float) * (is_complex ? 2 : 1);
- case 'd': return sizeof(double) * (is_complex ? 2 : 1);
- case 'g': return sizeof(long double) * (is_complex ? 2 : 1);
- case 'O': case 'P': return sizeof(void*);
- default: {
- __Pyx_BufFmt_RaiseUnexpectedChar(ch);
- return 0;
- }
- }
-}
-
-typedef struct { char c; short x; } __Pyx_st_short;
-typedef struct { char c; int x; } __Pyx_st_int;
-typedef struct { char c; long x; } __Pyx_st_long;
-typedef struct { char c; float x; } __Pyx_st_float;
-typedef struct { char c; double x; } __Pyx_st_double;
-typedef struct { char c; long double x; } __Pyx_st_longdouble;
-typedef struct { char c; void *x; } __Pyx_st_void_p;
-#ifdef HAVE_LONG_LONG
-typedef struct { char c; PY_LONG_LONG x; } __Pyx_s_long_long;
-#endif
-
-static size_t __Pyx_BufFmt_TypeCharToAlignment(char ch, int is_complex) {
- switch (ch) {
- case '?': case 'c': case 'b': case 'B': return 1;
- case 'h': case 'H': return sizeof(__Pyx_st_short) - sizeof(short);
- case 'i': case 'I': return sizeof(__Pyx_st_int) - sizeof(int);
- case 'l': case 'L': return sizeof(__Pyx_st_long) - sizeof(long);
-#ifdef HAVE_LONG_LONG
- case 'q': case 'Q': return sizeof(__Pyx_s_long_long) - sizeof(PY_LONG_LONG);
-#endif
- case 'f': return sizeof(__Pyx_st_float) - sizeof(float);
- case 'd': return sizeof(__Pyx_st_double) - sizeof(double);
- case 'g': return sizeof(__Pyx_st_longdouble) - sizeof(long double);
- case 'P': case 'O': return sizeof(__Pyx_st_void_p) - sizeof(void*);
- default:
- __Pyx_BufFmt_RaiseUnexpectedChar(ch);
- return 0;
- }
-}
-
-static size_t __Pyx_BufFmt_TypeCharToGroup(char ch, int is_complex) {
- switch (ch) {
- case 'c': case 'b': case 'h': case 'i': case 'l': case 'q': return 'I';
- case 'B': case 'H': case 'I': case 'L': case 'Q': return 'U';
- case 'f': case 'd': case 'g': return (is_complex ? 'C' : 'R');
- case 'O': return 'O';
- case 'P': return 'P';
- default: {
- __Pyx_BufFmt_RaiseUnexpectedChar(ch);
- return 0;
- }
- }
-}
-
-static void __Pyx_BufFmt_RaiseExpected(__Pyx_BufFmt_Context* ctx) {
- if (ctx->head == NULL || ctx->head->field == &ctx->root) {
- const char* expected;
- const char* quote;
- if (ctx->head == NULL) {
- expected = "end";
- quote = "";
- } else {
- expected = ctx->head->field->type->name;
- quote = "'";
- }
- PyErr_Format(PyExc_ValueError,
- "Buffer dtype mismatch, expected %s%s%s but got %s",
- quote, expected, quote,
- __Pyx_BufFmt_DescribeTypeChar(ctx->enc_type, ctx->is_complex));
- } else {
- __Pyx_StructField* field = ctx->head->field;
- __Pyx_StructField* parent = (ctx->head - 1)->field;
- PyErr_Format(PyExc_ValueError,
- "Buffer dtype mismatch, expected '%s' but got %s in '%s.%s'",
- field->type->name, __Pyx_BufFmt_DescribeTypeChar(ctx->enc_type, ctx->is_complex),
- parent->type->name, field->name);
- }
-}
-
-static int __Pyx_BufFmt_ProcessTypeChunk(__Pyx_BufFmt_Context* ctx) {
- char group;
- size_t size, offset;
- if (ctx->enc_type == 0) return 0;
- group = __Pyx_BufFmt_TypeCharToGroup(ctx->enc_type, ctx->is_complex);
- do {
- __Pyx_StructField* field = ctx->head->field;
- __Pyx_TypeInfo* type = field->type;
-
- if (ctx->packmode == '@' || ctx->packmode == '^') {
- size = __Pyx_BufFmt_TypeCharToNativeSize(ctx->enc_type, ctx->is_complex);
- } else {
- size = __Pyx_BufFmt_TypeCharToStandardSize(ctx->enc_type, ctx->is_complex);
- }
- if (ctx->packmode == '@') {
- int align_at = __Pyx_BufFmt_TypeCharToAlignment(ctx->enc_type, ctx->is_complex);
- int align_mod_offset;
- if (align_at == 0) return -1;
- align_mod_offset = ctx->fmt_offset % align_at;
- if (align_mod_offset > 0) ctx->fmt_offset += align_at - align_mod_offset;
- }
-
- if (type->size != size || type->typegroup != group) {
- if (type->typegroup == 'C' && type->fields != NULL) {
- /* special case -- treat as struct rather than complex number */
- size_t parent_offset = ctx->head->parent_offset + field->offset;
- ++ctx->head;
- ctx->head->field = type->fields;
- ctx->head->parent_offset = parent_offset;
- continue;
- }
-
- __Pyx_BufFmt_RaiseExpected(ctx);
- return -1;
- }
-
- offset = ctx->head->parent_offset + field->offset;
- if (ctx->fmt_offset != offset) {
- PyErr_Format(PyExc_ValueError,
- "Buffer dtype mismatch; next field is at offset %"PY_FORMAT_SIZE_T"d "
- "but %"PY_FORMAT_SIZE_T"d expected", ctx->fmt_offset, offset);
- return -1;
- }
-
- ctx->fmt_offset += size;
-
- --ctx->enc_count; /* Consume from buffer string */
-
- /* Done checking, move to next field, pushing or popping struct stack if needed */
- while (1) {
- if (field == &ctx->root) {
- ctx->head = NULL;
- if (ctx->enc_count != 0) {
- __Pyx_BufFmt_RaiseExpected(ctx);
- return -1;
- }
- break; /* breaks both loops as ctx->enc_count == 0 */
- }
- ctx->head->field = ++field;
- if (field->type == NULL) {
- --ctx->head;
- field = ctx->head->field;
- continue;
- } else if (field->type->typegroup == 'S') {
- size_t parent_offset = ctx->head->parent_offset + field->offset;
- if (field->type->fields->type == NULL) continue; /* empty struct */
- field = field->type->fields;
- ++ctx->head;
- ctx->head->field = field;
- ctx->head->parent_offset = parent_offset;
- break;
- } else {
- break;
- }
- }
- } while (ctx->enc_count);
- ctx->enc_type = 0;
- ctx->is_complex = 0;
- return 0;
-}
-
-static int __Pyx_BufFmt_FirstPack(__Pyx_BufFmt_Context* ctx) {
- if (ctx->enc_type != 0 || ctx->packmode != '@') {
- PyErr_SetString(PyExc_ValueError, "Buffer packing mode currently only allowed at beginning of format string (this is a defect)");
- return -1;
- }
- return 0;
-}
-
-static const char* __Pyx_BufFmt_CheckString(__Pyx_BufFmt_Context* ctx, const char* ts) {
- int got_Z = 0;
- while (1) {
- switch(*ts) {
- case 0:
- if (ctx->enc_type != 0 && ctx->head == NULL) {
- __Pyx_BufFmt_RaiseExpected(ctx);
- return NULL;
- }
- if (__Pyx_BufFmt_ProcessTypeChunk(ctx) == -1) return NULL;
- if (ctx->head != NULL) {
- __Pyx_BufFmt_RaiseExpected(ctx);
- return NULL;
- }
- return ts;
- case ' ':
- case 10:
- case 13:
- ++ts;
- break;
- case '<':
- if (!__Pyx_IsLittleEndian()) {
- PyErr_SetString(PyExc_ValueError, "Little-endian buffer not supported on big-endian compiler");
- return NULL;
- }
- if (__Pyx_BufFmt_FirstPack(ctx) == -1) return NULL;
- ctx->packmode = '=';
- ++ts;
- break;
- case '>':
- case '!':
- if (__Pyx_IsLittleEndian()) {
- PyErr_SetString(PyExc_ValueError, "Big-endian buffer not supported on little-endian compiler");
- return NULL;
- }
- if (__Pyx_BufFmt_FirstPack(ctx) == -1) return NULL;
- ctx->packmode = '=';
- ++ts;
- break;
- case '=':
- case '@':
- case '^':
- if (__Pyx_BufFmt_FirstPack(ctx) == -1) return NULL;
- ctx->packmode = *ts++;
- break;
- case 'T': /* substruct */
- {
- int i;
- const char* ts_after_sub;
- int struct_count = ctx->new_count;
- ctx->new_count = 1;
- ++ts;
- if (*ts != '{') {
- PyErr_SetString(PyExc_ValueError, "Buffer acquisition: Expected '{' after 'T'");
- return NULL;
- }
- ++ts;
- ts_after_sub = ts;
- for (i = 0; i != struct_count; ++i) {
- ts_after_sub = __Pyx_BufFmt_CheckString(ctx, ts);
- if (!ts_after_sub) return NULL;
- }
- ts = ts_after_sub;
- }
- break;
- case '}': /* end of substruct; either repeat or move on */
- ++ts;
- return ts;
- case 'x':
- if (__Pyx_BufFmt_ProcessTypeChunk(ctx) == -1) return NULL;
- ctx->fmt_offset += ctx->new_count;
- ctx->new_count = 1;
- ctx->enc_count = 0;
- ctx->enc_type = 0;
- ++ts;
- break;
- case 'Z':
- got_Z = 1;
- ++ts;
- if (*ts != 'f' && *ts != 'd' && *ts != 'g') {
- __Pyx_BufFmt_RaiseUnexpectedChar('Z');
- return NULL;
- } /* fall through */
- case 'c': case 'b': case 'B': case 'h': case 'H': case 'i': case 'I':
- case 'l': case 'L': case 'q': case 'Q':
- case 'f': case 'd': case 'g':
- case 'O':
- if (ctx->enc_type == *ts && got_Z == ctx->is_complex) {
- /* Continue pooling same type */
- ctx->enc_count += ctx->new_count;
- } else {
- /* New type */
- if (__Pyx_BufFmt_ProcessTypeChunk(ctx) == -1) return NULL;
- ctx->enc_count = ctx->new_count;
- ctx->enc_type = *ts;
- ctx->is_complex = got_Z;
- }
- ++ts;
- ctx->new_count = 1;
- got_Z = 0;
- break;
- default:
- {
- ctx->new_count = __Pyx_BufFmt_ParseNumber(&ts);
- if (ctx->new_count == -1) { /* First char was not a digit */
- char msg[2] = { *ts, 0 };
- PyErr_Format(PyExc_ValueError,
- "Does not understand character buffer dtype format string ('%s')", msg);
- return NULL;
- }
- }
-
- }
- }
-}
-
-static CYTHON_INLINE void __Pyx_ZeroBuffer(Py_buffer* buf) {
- buf->buf = NULL;
- buf->obj = NULL;
- buf->strides = __Pyx_zeros;
- buf->shape = __Pyx_zeros;
- buf->suboffsets = __Pyx_minusones;
-}
-
-static int __Pyx_GetBufferAndValidate(Py_buffer* buf, PyObject* obj, __Pyx_TypeInfo* dtype, int flags, int nd, int cast, __Pyx_BufFmt_StackElem* stack) {
- if (obj == Py_None) {
- __Pyx_ZeroBuffer(buf);
- return 0;
- }
- buf->buf = NULL;
- if (__Pyx_GetBuffer(obj, buf, flags) == -1) goto fail;
- if (buf->ndim != nd) {
- PyErr_Format(PyExc_ValueError,
- "Buffer has wrong number of dimensions (expected %d, got %d)",
- nd, buf->ndim);
- goto fail;
- }
- if (!cast) {
- __Pyx_BufFmt_Context ctx;
- __Pyx_BufFmt_Init(&ctx, stack, dtype);
- if (!__Pyx_BufFmt_CheckString(&ctx, buf->format)) goto fail;
- }
- if ((unsigned)buf->itemsize != dtype->size) {
- PyErr_Format(PyExc_ValueError,
- "Item size of buffer (%"PY_FORMAT_SIZE_T"d byte%s) does not match size of '%s' (%"PY_FORMAT_SIZE_T"d byte%s)",
- buf->itemsize, (buf->itemsize > 1) ? "s" : "",
- dtype->name,
- dtype->size, (dtype->size > 1) ? "s" : "");
- goto fail;
- }
- if (buf->suboffsets == NULL) buf->suboffsets = __Pyx_minusones;
- return 0;
-fail:;
- __Pyx_ZeroBuffer(buf);
- return -1;
-}
-
-static CYTHON_INLINE void __Pyx_SafeReleaseBuffer(Py_buffer* info) {
- if (info->buf == NULL) return;
- if (info->suboffsets == __Pyx_minusones) info->suboffsets = NULL;
- __Pyx_ReleaseBuffer(info);
-}
-
-static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb) {
- PyObject *tmp_type, *tmp_value, *tmp_tb;
- PyThreadState *tstate = PyThreadState_GET();
-
- tmp_type = tstate->curexc_type;
- tmp_value = tstate->curexc_value;
- tmp_tb = tstate->curexc_traceback;
- tstate->curexc_type = type;
- tstate->curexc_value = value;
- tstate->curexc_traceback = tb;
- Py_XDECREF(tmp_type);
- Py_XDECREF(tmp_value);
- Py_XDECREF(tmp_tb);
-}
-
-static CYTHON_INLINE void __Pyx_ErrFetch(PyObject **type, PyObject **value, PyObject **tb) {
- PyThreadState *tstate = PyThreadState_GET();
- *type = tstate->curexc_type;
- *value = tstate->curexc_value;
- *tb = tstate->curexc_traceback;
-
- tstate->curexc_type = 0;
- tstate->curexc_value = 0;
- tstate->curexc_traceback = 0;
-}
-
-
static CYTHON_INLINE void __Pyx_RaiseNeedMoreValuesError(Py_ssize_t index) {
PyErr_Format(PyExc_ValueError,
#if PY_VERSION_HEX < 0x02050000
@@ -4143,30 +3794,6 @@ static int __Pyx_ArgTypeTest(PyObject *obj, PyTypeObject *type, int none_allowed
return 0;
}
-#if PY_MAJOR_VERSION < 3
-static int __Pyx_GetBuffer(PyObject *obj, Py_buffer *view, int flags) {
- #if PY_VERSION_HEX >= 0x02060000
- if (Py_TYPE(obj)->tp_flags & Py_TPFLAGS_HAVE_NEWBUFFER)
- return PyObject_GetBuffer(obj, view, flags);
- #endif
- if (PyObject_TypeCheck(obj, __pyx_ptype_5numpy_ndarray)) return __pyx_pf_5numpy_7ndarray___getbuffer__(obj, view, flags);
- else {
- PyErr_Format(PyExc_TypeError, "'%100s' does not have the buffer interface", Py_TYPE(obj)->tp_name);
- return -1;
- }
-}
-
-static void __Pyx_ReleaseBuffer(Py_buffer *view) {
- PyObject* obj = view->obj;
- if (obj) {
-if (PyObject_TypeCheck(obj, __pyx_ptype_5numpy_ndarray)) __pyx_pf_5numpy_7ndarray___releasebuffer__(obj, view);
- Py_DECREF(obj);
- view->obj = NULL;
- }
-}
-
-#endif
-
static PyObject *__Pyx_Import(PyObject *name, PyObject *from_list) {
PyObject *__import__ = 0;
PyObject *empty_list = 0;
@@ -4208,6 +3835,33 @@ static PyObject *__Pyx_GetName(PyObject *dict, PyObject *name) {
return result;
}
+static CYTHON_INLINE void __Pyx_ErrRestore(PyObject *type, PyObject *value, PyObject *tb) {
+ PyObject *tmp_type, *tmp_value, *tmp_tb;
+ PyThreadState *tstate = PyThreadState_GET();
+
+ tmp_type = tstate->curexc_type;
+ tmp_value = tstate->curexc_value;
+ tmp_tb = tstate->curexc_traceback;
+ tstate->curexc_type = type;
+ tstate->curexc_value = value;
+ tstate->curexc_traceback = tb;
+ Py_XDECREF(tmp_type);
+ Py_XDECREF(tmp_value);
+ Py_XDECREF(tmp_tb);
+}
+
+static CYTHON_INLINE void __Pyx_ErrFetch(PyObject **type, PyObject **value, PyObject **tb) {
+ PyThreadState *tstate = PyThreadState_GET();
+ *type = tstate->curexc_type;
+ *value = tstate->curexc_value;
+ *tb = tstate->curexc_traceback;
+
+ tstate->curexc_type = 0;
+ tstate->curexc_value = 0;
+ tstate->curexc_traceback = 0;
+}
+
+
#if PY_MAJOR_VERSION < 3
static void __Pyx_Raise(PyObject *type, PyObject *value, PyObject *tb) {
Py_XINCREF(type);
diff --git a/scikits/learn/utils/arrayfuncs.pyx b/scikits/learn/utils/arrayfuncs.pyx
index eaffb35208..d75f6cc447 100644
--- a/scikits/learn/utils/arrayfuncs.pyx
+++ b/scikits/learn/utils/arrayfuncs.pyx
@@ -28,13 +28,18 @@ cdef extern from "cblas.h":
int N, double *A, int lda, double *X,
int incX)
+ void cblas_strsv(CBLAS_ORDER Order, CBLAS_UPLO Uplo,
+ CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
+ int N, float *A, int lda, float *X,
+ int incX)
cdef extern from "float.h":
cdef double DBL_MAX
cdef float FLT_MAX
cdef extern from "src/cholesky_delete.c":
- int c_cholesky_delete "cholesky_delete" (int m, int n, double *L, int go_out)
+ int double_cholesky_delete (int m, int n, double *L, int go_out)
+ int float_cholesky_delete (int m, int n, float *L, int go_out)
ctypedef np.float64_t DOUBLE
@@ -70,26 +75,41 @@ cdef double _double_min_pos (double *X, Py_ssize_t size):
min_val = X[i]
return min_val
-def solve_triangular (
- np.ndarray[DOUBLE, ndim=2] X,
- np.ndarray[DOUBLE, ndim=1] y
- ):
+def solve_triangular (np.ndarray X, np.ndarray y):
"""
Solves a triangular system (overwrites y)
- TODO: option for upper, lower, etc.
+ Note: The lapack function to solve triangular systems was added to
+ scipy v0.9. Remove this when we stop supporting earlier versions.
"""
- cdef int lda = <int> X.strides[0] / sizeof(double)
+ cdef int lda
+
+ if X.dtype.name == 'float64' and y.dtype.name == 'float64':
+ lda = <int> X.strides[0] / sizeof(double)
+
+ cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ CblasNonUnit, <int> X.shape[0], <double *> X.data,
+ lda, <double *> y.data, 1);
- cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
- CblasNonUnit, <int> X.shape[0], <double *> X.data,
- lda, <double *> y.data, 1);
+ elif X.dtype.name == 'float32' and y.dtype.name == 'float32':
+ lda = <int> X.strides[0] / sizeof(float)
+ cblas_strsv (CblasRowMajor, CblasLower, CblasNoTrans,
+ CblasNonUnit, <int> X.shape[0], <float *> X.data,
+ lda, <float *> y.data, 1);
+ else:
+ raise ValueError ('Unsupported or inconsistent dtype in arrays X, y')
-def cholesky_delete (
- np.ndarray[DOUBLE, ndim=2] L,
- int go_out):
+
+def cholesky_delete (np.ndarray L, int go_out):
cdef int n = <int> L.shape[0]
- cdef int m = <int> L.strides[0] / sizeof (double)
- c_cholesky_delete (m, n, <double *> L.data, go_out)
+ cdef int m
+
+ if L.dtype.name == 'float64':
+ m = <int> L.strides[0] / sizeof (double)
+ double_cholesky_delete (m, n, <double *> L.data, go_out)
+ elif L.dtype.name == 'float32':
+ m = <int> L.strides[0] / sizeof (float)
+ float_cholesky_delete (m, n, <float *> L.data, go_out)
+
diff --git a/scikits/learn/utils/src/cholesky_delete.c b/scikits/learn/utils/src/cholesky_delete.c
index fe029de4b8..7117ba0a5c 100644
--- a/scikits/learn/utils/src/cholesky_delete.c
+++ b/scikits/learn/utils/src/cholesky_delete.c
@@ -13,7 +13,7 @@
* TODO: put transpose as an option
*
*/
-int cholesky_delete (int m, int n, double *L, int go_out) {
+int double_cholesky_delete (int m, int n, double *L, int go_out) {
double c, s;
@@ -45,3 +45,37 @@ int cholesky_delete (int m, int n, double *L, int go_out) {
return 0;
}
+
+
+int float_cholesky_delete (int m, int n, float *L, int go_out) {
+
+ float c, s;
+
+ /* delete row go_out */
+ float * _L = L + (go_out * m);
+ int i;
+ for (i = go_out; i < n - 1; ++i) {
+ cblas_scopy (i + 2, _L + m , 1, _L, 1);
+ _L += m;
+ }
+
+ _L = L + (go_out * m);
+ for (i=go_out; i < n - 1; ++i) {
+
+ cblas_srotg (_L + i, _L + i + 1, &c, &s);
+ if (_L[i] < 0) {
+ /* Diagonals cannot be negative */
+ /* _L[i] = copysign(_L[i], 1.0); */
+ c = -c;
+ s = -s;
+ }
+ _L[i+1] = 0.; /* just for cleanup */
+ _L += m;
+
+ cblas_srot (n - (i + 2), _L + i, m, _L + i + 1,
+ m, c, s);
+
+ }
+
+ return 0;
+}
--
GitLab