mirror/opencv
1
0
Fork 0
mirror of https://github.com/opencv/opencv.git synced 2025-01-18 06:03:15 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

optimized lapack' SVD for noticeably better performance on small matrices

Browse Source
This commit is contained in:
Vadim Pisarevsky 2010-08-30 16:37:22 +00:00
parent fea66d9384
commit e48a456d48
19 changed files with 957 additions and 3620 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
3rdparty/include/cblas.h vendored
View File

@ -37,11 +37,28 @@ static __inline double r_sign(real *a, real *b)
return *b >= 0 ? x : -x;
}
extern const unsigned char lapack_toupper_tab[];
#define lapack_toupper(c) ((char)lapack_toupper_tab[(unsigned char)(c)])
extern const unsigned char lapack_lamch_tab[];
extern const doublereal lapack_dlamch_tab[];
extern const doublereal lapack_slamch_tab[];
static __inline logical lsame_(char *ca, char *cb)
{
return toupper(ca[0]) == toupper(cb[0]);
return lapack_toupper(ca[0]) == lapack_toupper(cb[0]);
}
static __inline doublereal dlamch_(char* cmach)
{
return lapack_dlamch_tab[lapack_lamch_tab[(unsigned char)cmach[0]]];
}
static __inline doublereal slamch_(char* cmach)
{
return lapack_slamch_tab[lapack_lamch_tab[(unsigned char)cmach[0]]];
}
static __inline integer i_nint(real *x)
{
return (integer)(*x >= 0 ? floor(*x + .5) : -floor(.5 - *x));

7
3rdparty/include/clapack.h vendored
View File

@ -3680,8 +3680,6 @@ doublereal dsecnd_();
doublereal second_();
doublereal slamch_(char *cmach);
/* Subroutine */ int slamc1_(integer *beta, integer *t, logical *rnd, logical
*ieee1);
@ -3696,8 +3694,6 @@ doublereal slamc3_(real *a, real *b);
logical *ieee, integer *emax, real *rmax);
doublereal dlamch_(char *cmach);
/* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical
*ieee1);
@ -3712,9 +3708,6 @@ doublereal dlamc3_(doublereal *a, doublereal *b);
/* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin,
logical *ieee, integer *emax, doublereal *rmax);
integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4);
#ifdef __cplusplus
}
#endif

5
3rdparty/include/f2c.h vendored
View File

@ -7,6 +7,7 @@
#ifndef F2C_INCLUDE
#define F2C_INCLUDE
#include <assert.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
@ -17,6 +18,10 @@
#include <string.h>
#include <stdio.h>
#if __SSE2__ || defined _M_X64
#include "emmintrin.h"
#endif
#ifdef __cplusplus
extern "C" {
#endif

312
3rdparty/lapack/dgemv.c vendored
View File

@ -1,312 +0,0 @@
/* dgemv.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
/* Subroutine */ int dgemv_(char *trans, integer *m, integer *n, doublereal *
alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
doublereal *beta, doublereal *y, integer *incy)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
integer i__, j, ix, iy, jx, jy, kx, ky, info;
doublereal temp;
integer lenx, leny;
extern logical lsame_(char *, char *);
extern /* Subroutine */ int xerbla_(char *, integer *);
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Arguments */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared */
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - DOUBLE PRECISION. */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the */
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- Written on 22-October-1986. */
/* Jack Dongarra, Argonne National Lab. */
/* Jeremy Du Croz, Nag Central Office. */
/* Sven Hammarling, Nag Central Office. */
/* Richard Hanson, Sandia National Labs. */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--x;
--y;
/* Function Body */
info = 0;
if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
) {
info = 1;
} else if (*m < 0) {
info = 2;
} else if (*n < 0) {
info = 3;
} else if (*lda < max(1,*m)) {
info = 6;
} else if (*incx == 0) {
info = 8;
} else if (*incy == 0) {
info = 11;
}
if (info != 0) {
xerbla_("DGEMV ", &info);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) {
return 0;
}
/* Set LENX and LENY, the lengths of the vectors x and y, and set */
/* up the start points in X and Y. */
if (lsame_(trans, "N")) {
lenx = *n;
leny = *m;
} else {
lenx = *m;
leny = *n;
}
if (*incx > 0) {
kx = 1;
} else {
kx = 1 - (lenx - 1) * *incx;
}
if (*incy > 0) {
ky = 1;
} else {
ky = 1 - (leny - 1) * *incy;
}
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
/* First form y := beta*y. */
if (*beta != 1.) {
if (*incy == 1) {
if (*beta == 0.) {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[i__] = 0.;
/* L10: */
}
} else {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[i__] = *beta * y[i__];
/* L20: */
}
}
} else {
iy = ky;
if (*beta == 0.) {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[iy] = 0.;
iy += *incy;
/* L30: */
}
} else {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[iy] = *beta * y[iy];
iy += *incy;
/* L40: */
}
}
}
}
if (*alpha == 0.) {
return 0;
}
if (lsame_(trans, "N")) {
/* Form y := alpha*A*x + y. */
jx = kx;
if (*incy == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (x[jx] != 0.) {
temp = *alpha * x[jx];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
y[i__] += temp * a[i__ + j * a_dim1];
/* L50: */
}
}
jx += *incx;
/* L60: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (x[jx] != 0.) {
temp = *alpha * x[jx];
iy = ky;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
y[iy] += temp * a[i__ + j * a_dim1];
iy += *incy;
/* L70: */
}
}
jx += *incx;
/* L80: */
}
}
} else {
/* Form y := alpha*A'*x + y. */
jy = ky;
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp = 0.;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp += a[i__ + j * a_dim1] * x[i__];
/* L90: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L100: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp = 0.;
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp += a[i__ + j * a_dim1] * x[ix];
ix += *incx;
/* L110: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L120: */
}
}
}
return 0;
/* End of DGEMV . */
} /* dgemv_ */

238
3rdparty/lapack/dgemv_custom.c vendored Normal file
View File

@ -0,0 +1,238 @@
#include "clapack.h"
/* Subroutine */ int dgemv_(char *_trans, integer *_m, integer *_n, doublereal *
_alpha, doublereal *a, integer *_lda, doublereal *x, integer *_incx,
doublereal *_beta, doublereal *y, integer *_incy)
{
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Arguments */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - DOUBLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared */
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - DOUBLE PRECISION. */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - DOUBLE PRECISION array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the */
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- Written on 22-October-1986. */
/* Jack Dongarra, Argonne National Lab. */
/* Jeremy Du Croz, Nag Central Office. */
/* Sven Hammarling, Nag Central Office. */
/* Richard Hanson, Sandia National Labs. */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* Test the input parameters. */
char trans = lapack_toupper(_trans[0]);
integer i, j, m = *_m, n = *_n, lda = *_lda, incx = *_incx, incy = *_incy;
integer leny = trans == 'N' ? m : n, lenx = trans == 'N' ? n : m;
real alpha = *_alpha, beta = *_beta;
integer info = 0;
if (trans != 'N' && trans != 'T' && trans != 'C')
info = 1;
else if (m < 0)
info = 2;
else if (n < 0)
info = 3;
else if (lda < max(1,m))
info = 6;
else if (incx == 0)
info = 8;
else if (incy == 0)
info = 11;
if (info != 0)
{
xerbla_("SGEMV ", &info);
return 0;
}
if( incy < 0 )
y -= incy*(leny - 1);
if( incx < 0 )
x -= incx*(lenx - 1);
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if( beta != 1. )
{
if( incy == 1 )
{
if( beta == 0. )
for( i = 0; i < leny; i++ )
y[i] = 0.;
else
for( i = 0; i < leny; i++ )
y[i] *= beta;
}
else
{
if( beta == 0. )
for( i = 0; i < leny; i++ )
y[i*incy] = 0.;
else
for( i = 0; i < leny; i++ )
y[i*incy] *= beta;
}
}
if( alpha == 0. )
;
else if( trans == 'N' )
{
if( incy == 1 )
{
for( i = 0; i < n; i++, a += lda )
{
doublereal s = x[i*incx];
if( s == 0. )
continue;
s *= alpha;
for( j = 0; j <= m - 2; j += 2 )
{
doublereal t0 = y[j] + s*a[j];
doublereal t1 = y[j+1] + s*a[j+1];
y[j] = t0; y[j+1] = t1;
}
for( ; j < m; j++ )
y[j] += s*a[j];
}
}
else
{
for( i = 0; i < n; i++, a += lda )
{
doublereal s = x[i*incx];
if( s == 0. )
continue;
s *= alpha;
for( j = 0; j < m; j++ )
y[j*incy] += s*a[j];
}
}
}
else
{
if( incx == 1 )
{
for( i = 0; i < n; i++, a += lda )
{
doublereal s = 0;
for( j = 0; j <= m - 2; j += 2 )
s += x[j]*a[j] + x[j+1]*a[j+1];
for( ; j < m; j++ )
s += x[j]*a[j];
y[i*incy] += alpha*s;
}
}
else
{
for( i = 0; i < n; i++, a += lda )
{
doublereal s = 0;
for( j = 0; j < m; j++ )
s += x[j*incx]*a[j];
y[i*incy] += alpha*s;
}
}
}
return 0;
/* End of DGEMV . */
} /* dgemv_ */

153
3rdparty/lapack/dger.c → 3rdparty/lapack/dger_custom.c vendored
View File

@ -1,29 +1,10 @@
/* dger.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
/* Subroutine */ int dger_(integer *m, integer *n, doublereal *alpha,
doublereal *x, integer *incx, doublereal *y, integer *incy,
doublereal *a, integer *lda)
/* Subroutine */ int dger_(integer *_m, integer *_n, doublereal *_alpha,
doublereal *x, integer *_incx, doublereal *y, integer *_incy,
doublereal *a, integer *_lda)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
integer i__, j, ix, jy, kx, info;
doublereal temp;
extern /* Subroutine */ int xerbla_(char *, integer *);
/* .. Scalar Arguments .. */
/* .. */
@ -111,80 +92,70 @@
/* Test the input parameters. */
/* Parameter adjustments */
--x;
--y;
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
info = 0;
if (*m < 0) {
info = 1;
} else if (*n < 0) {
info = 2;
} else if (*incx == 0) {
info = 5;
} else if (*incy == 0) {
info = 7;
} else if (*lda < max(1,*m)) {
info = 9;
}
if (info != 0) {
xerbla_("DGER ", &info);
return 0;
integer i, j, m = *_m, n = *_n, incx = *_incx, incy = *_incy, lda = *_lda;
doublereal alpha = *_alpha;
integer info = 0;
if (m < 0)
info = 1;
else if (n < 0)
info = 2;
else if (incx == 0)
info = 5;
else if (incy == 0)
info = 7;
else if (lda < max(1,m))
info = 9;
if (info != 0)
{
xerbla_("DGER ", &info);
return 0;
}
/* Quick return if possible. */
if (incx < 0)
x -= (m-1)*incx;
if (incy < 0)
y -= (n-1)*incy;
if (*m == 0 || *n == 0 || *alpha == 0.) {
return 0;
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if( alpha == 0 )
;
else if( incx == 1 )
{
for( j = 0; j < n; j++, a += lda )
{
doublereal s = y[j*incy];
if( s == 0 )
continue;
s *= alpha;
for( i = 0; i <= m - 2; i += 2 )
{
doublereal t0 = a[i] + x[i]*s;
doublereal t1 = a[i+1] + x[i+1]*s;
a[i] = t0; a[i+1] = t1;
}
for( ; i < m; i++ )
a[i] += x[i]*s;
}
}
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if (*incy > 0) {
jy = 1;
} else {
jy = 1 - (*n - 1) * *incy;
}
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (y[jy] != 0.) {
temp = *alpha * y[jy];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
a[i__ + j * a_dim1] += x[i__] * temp;
/* L10: */
}
}
jy += *incy;
/* L20: */
}
} else {
if (*incx > 0) {
kx = 1;
} else {
kx = 1 - (*m - 1) * *incx;
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (y[jy] != 0.) {
temp = *alpha * y[jy];
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
a[i__ + j * a_dim1] += x[ix] * temp;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
else
{
for( j = 0; j < n; j++, a += lda )
{
doublereal s = y[j*incy];
if( s == 0 )
continue;
s *= alpha;
for( i = 0; i < m; i++ )
a[i] += x[i*incx]*s;
}
}
return 0;

1047
3rdparty/lapack/dlamch.c vendored
View File

File diff suppressed because it is too large Load Diff

58
3rdparty/lapack/dlamch_custom.c vendored Normal file
View File

@ -0,0 +1,58 @@
#include "clapack.h"
#include <float.h>
#include <stdio.h>
/* *********************************************************************** */
doublereal dlamc3_(doublereal *a, doublereal *b)
{
/* System generated locals */
doublereal ret_val;
/* -- LAPACK auxiliary routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLAMC3 is intended to force A and B to be stored prior to doing */
/* the addition of A and B , for use in situations where optimizers */
/* might hold one of these in a register. */
/* Arguments */
/* ========= */
/* A (input) DOUBLE PRECISION */
/* B (input) DOUBLE PRECISION */
/* The values A and B. */
/* ===================================================================== */
/* .. Executable Statements .. */
ret_val = *a + *b;
return ret_val;
/* End of DLAMC3 */
} /* dlamc3_ */
/* simpler version of dlamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
#ifndef DBL_DIGITS
#define DBL_DIGITS 53
#endif
const doublereal lapack_dlamch_tab[] =
{
0, FLT_RADIX, DBL_EPSILON, DBL_MAX_EXP, DBL_MIN_EXP, DBL_DIGITS, DBL_MAX,
DBL_EPSILON*FLT_RADIX, 1, DBL_MIN*(1 + DBL_EPSILON), DBL_MIN
};

42
3rdparty/lapack/dlartg.c → 3rdparty/lapack/dlartg_custom.c vendored
View File

@ -1,15 +1,3 @@
/* dlartg.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
@ -20,17 +8,14 @@
integer i__1;
doublereal d__1, d__2;
/* Builtin functions */
double log(doublereal), pow_di(doublereal *, integer *), sqrt(doublereal);
/* Local variables */
integer i__;
doublereal f1, g1, eps, scale;
integer count;
doublereal safmn2, safmx2;
extern doublereal dlamch_(char *);
doublereal safmin;
static doublereal safmn2, safmx2;
static doublereal safmin;
static volatile logical FIRST = TRUE_;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
@ -97,15 +82,16 @@
/* .. */
/* .. Executable Statements .. */
/* IF( FIRST ) THEN */
safmin = dlamch_("S");
eps = dlamch_("E");
d__1 = dlamch_("B");
i__1 = (integer) (log(safmin / eps) / log(dlamch_("B")) / 2.);
safmn2 = pow_di(&d__1, &i__1);
safmx2 = 1. / safmn2;
/* FIRST = .FALSE. */
/* END IF */
if( FIRST )
{
safmin = dlamch_("S");
eps = dlamch_("E");
d__1 = dlamch_("B");
i__1 = (integer) (log(safmin / eps) / log(dlamch_("B")) / 2.);
safmn2 = pow_di(&d__1, &i__1);
safmx2 = 1. / safmn2;
FIRST = FALSE_;
}
if (*g == 0.) {
*cs = 1.;
*sn = 0.;

0
3rdparty/lapack/dlasr.c → 3rdparty/lapack/dlasr_custom.c vendored
View File

654
3rdparty/lapack/ilaenv.c vendored
View File

@ -1,654 +0,0 @@
/* ilaenv.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
#include "string.h"
/* Table of constant values */
static integer c__1 = 1;
static real c_b163 = 0.f;
static real c_b164 = 1.f;
static integer c__0 = 0;
integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4)
{
/* System generated locals */
integer ret_val;
/* Builtin functions */
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
integer s_cmp(char *, char *, ftnlen, ftnlen);
/* Local variables */
integer i__;
char c1[1], c2[2], c3[3], c4[2];
integer ic, nb, iz, nx;
logical cname;
integer nbmin;
logical sname;
extern integer ieeeck_(integer *, real *, real *);
char subnam[6];
extern integer iparmq_(integer *, char *, char *, integer *, integer *,
integer *, integer *);
ftnlen name_len, opts_len;
name_len = (ftnlen)strlen (name__);
opts_len = (ftnlen)strlen (opts);
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* January 2007 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ILAENV is called from the LAPACK routines to choose problem-dependent */
/* parameters for the local environment. See ISPEC for a description of */
/* the parameters. */
/* ILAENV returns an INTEGER */
/* if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC */
/* if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value. */
/* This version provides a set of parameters which should give good, */
/* but not optimal, performance on many of the currently available */
/* computers. Users are encouraged to modify this subroutine to set */
/* the tuning parameters for their particular machine using the option */
/* and problem size information in the arguments. */
/* This routine will not function correctly if it is converted to all */
/* lower case. Converting it to all upper case is allowed. */
/* Arguments */
/* ========= */
/* ISPEC (input) INTEGER */
/* Specifies the parameter to be returned as the value of */
/* ILAENV. */
/* = 1: the optimal blocksize; if this value is 1, an unblocked */
/* algorithm will give the best performance. */
/* = 2: the minimum block size for which the block routine */
/* should be used; if the usable block size is less than */
/* this value, an unblocked routine should be used. */
/* = 3: the crossover point (in a block routine, for N less */
/* than this value, an unblocked routine should be used) */
/* = 4: the number of shifts, used in the nonsymmetric */
/* eigenvalue routines (DEPRECATED) */
/* = 5: the minimum column dimension for blocking to be used; */
/* rectangular blocks must have dimension at least k by m, */
/* where k is given by ILAENV(2,...) and m by ILAENV(5,...) */
/* = 6: the crossover point for the SVD (when reducing an m by n */
/* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds */
/* this value, a QR factorization is used first to reduce */
/* the matrix to a triangular form.) */
/* = 7: the number of processors */
/* = 8: the crossover point for the multishift QR method */
/* for nonsymmetric eigenvalue problems (DEPRECATED) */
/* = 9: maximum size of the subproblems at the bottom of the */
/* computation tree in the divide-and-conquer algorithm */
/* (used by xGELSD and xGESDD) */
/* =10: ieee NaN arithmetic can be trusted not to trap */
/* =11: infinity arithmetic can be trusted not to trap */
/* 12 <= ISPEC <= 16: */
/* xHSEQR or one of its subroutines, */
/* see IPARMQ for detailed explanation */
/* NAME (input) CHARACTER*(*) */
/* The name of the calling subroutine, in either upper case or */
/* lower case. */
/* OPTS (input) CHARACTER*(*) */
/* The character options to the subroutine NAME, concatenated */
/* into a single character string. For example, UPLO = 'U', */
/* TRANS = 'T', and DIAG = 'N' for a triangular routine would */
/* be specified as OPTS = 'UTN'. */
/* N1 (input) INTEGER */
/* N2 (input) INTEGER */
/* N3 (input) INTEGER */
/* N4 (input) INTEGER */
/* Problem dimensions for the subroutine NAME; these may not all */
/* be required. */
/* Further Details */
/* =============== */
/* The following conventions have been used when calling ILAENV from the */
/* LAPACK routines: */
/* 1) OPTS is a concatenation of all of the character options to */
/* subroutine NAME, in the same order that they appear in the */
/* argument list for NAME, even if they are not used in determining */
/* the value of the parameter specified by ISPEC. */
/* 2) The problem dimensions N1, N2, N3, N4 are specified in the order */
/* that they appear in the argument list for NAME. N1 is used */
/* first, N2 second, and so on, and unused problem dimensions are */
/* passed a value of -1. */
/* 3) The parameter value returned by ILAENV is checked for validity in */
/* the calling subroutine. For example, ILAENV is used to retrieve */
/* the optimal blocksize for STRTRI as follows: */
/* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) */
/* IF( NB.LE.1 ) NB = MAX( 1, N ) */
/* ===================================================================== */
/* .. Local Scalars .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Executable Statements .. */
switch (*ispec) {
case 1: goto L10;
case 2: goto L10;
case 3: goto L10;
case 4: goto L80;
case 5: goto L90;
case 6: goto L100;
case 7: goto L110;
case 8: goto L120;
case 9: goto L130;
case 10: goto L140;
case 11: goto L150;
case 12: goto L160;
case 13: goto L160;
case 14: goto L160;
case 15: goto L160;
case 16: goto L160;
}
/* Invalid value for ISPEC */
ret_val = -1;
return ret_val;
L10:
/* Convert NAME to upper case if the first character is lower case. */
ret_val = 1;
s_copy(subnam, name__, (ftnlen)1, name_len);
ic = *(unsigned char *)subnam;
iz = 'Z';
if (iz == 90 || iz == 122) {
/* ASCII character set */
if (ic >= 97 && ic <= 122) {
*(unsigned char *)subnam = (char) (ic - 32);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 97 && ic <= 122) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
}
/* L20: */
}
}
} else if (iz == 233 || iz == 169) {
/* EBCDIC character set */
if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 &&
ic <= 169) {
*(unsigned char *)subnam = (char) (ic + 64);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >=
162 && ic <= 169) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64);
}
/* L30: */
}
}
} else if (iz == 218 || iz == 250) {
/* Prime machines: ASCII+128 */
if (ic >= 225 && ic <= 250) {
*(unsigned char *)subnam = (char) (ic - 32);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 225 && ic <= 250) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
}
/* L40: */
}
}
}
*(unsigned char *)c1 = *(unsigned char *)subnam;
sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D';
cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z';
if (! (cname || sname)) {
return ret_val;
}
s_copy(c2, subnam + 1, (ftnlen)1, (ftnlen)2);
s_copy(c3, subnam + 3, (ftnlen)1, (ftnlen)3);
s_copy(c4, c3 + 1, (ftnlen)1, (ftnlen)2);
switch (*ispec) {
case 1: goto L50;
case 2: goto L60;
case 3: goto L70;
}
L50:
/* ISPEC = 1: block size */
/* In these examples, separate code is provided for setting NB for */
/* real and complex. We assume that NB will take the same value in */
/* single or double precision. */
nb = 1;
if (s_cmp(c2, "GE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
} else if (s_cmp(c3, "QRF", (ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3,
"RQF", (ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)
1, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)1, (ftnlen)3)
== 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "HRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "BRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "TRI", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "PO", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
} else if (sname && s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nb = 32;
} else if (sname && s_cmp(c3, "GST", (ftnlen)1, (ftnlen)3) == 0) {
nb = 64;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
nb = 64;
} else if (s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nb = 32;
} else if (s_cmp(c3, "GST", (ftnlen)1, (ftnlen)3) == 0) {
nb = 64;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nb = 32;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nb = 32;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nb = 32;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nb = 32;
}
}
} else if (s_cmp(c2, "GB", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
if (*n4 <= 64) {
nb = 1;
} else {
nb = 32;
}
} else {
if (*n4 <= 64) {
nb = 1;
} else {
nb = 32;
}
}
}
} else if (s_cmp(c2, "PB", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
if (*n2 <= 64) {
nb = 1;
} else {
nb = 32;
}
} else {
if (*n2 <= 64) {
nb = 1;
} else {
nb = 32;
}
}
}
} else if (s_cmp(c2, "TR", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRI", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "LA", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "UUM", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (sname && s_cmp(c2, "ST", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "EBZ", (ftnlen)1, (ftnlen)3) == 0) {
nb = 1;
}
}
ret_val = nb;
return ret_val;
L60:
/* ISPEC = 2: minimum block size */
nbmin = 2;
if (s_cmp(c2, "GE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "QRF", (ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)1, (
ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)1, (ftnlen)3) == 0)
{
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "HRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "BRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "TRI", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nbmin = 8;
} else {
nbmin = 8;
}
} else if (sname && s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nbmin = 2;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nbmin = 2;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nbmin = 2;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nbmin = 2;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nbmin = 2;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nbmin = 2;
}
}
}
ret_val = nbmin;
return ret_val;
L70:
/* ISPEC = 3: crossover point */
nx = 0;
if (s_cmp(c2, "GE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "QRF", (ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
ftnlen)1, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)1, (
ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)1, (ftnlen)3) == 0)
{
if (sname) {
nx = 128;
} else {
nx = 128;
}
} else if (s_cmp(c3, "HRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nx = 128;
} else {
nx = 128;
}
} else if (s_cmp(c3, "BRD", (ftnlen)1, (ftnlen)3) == 0) {
if (sname) {
nx = 128;
} else {
nx = 128;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)1, (ftnlen)2) == 0) {
if (sname && s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nx = 32;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)1, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRD", (ftnlen)1, (ftnlen)3) == 0) {
nx = 32;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nx = 128;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)1, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)1, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)1, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)1, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)1, (ftnlen)2) == 0) {
nx = 128;
}
}
}
ret_val = nx;
return ret_val;
L80:
/* ISPEC = 4: number of shifts (used by xHSEQR) */
ret_val = 6;
return ret_val;
L90:
/* ISPEC = 5: minimum column dimension (not used) */
ret_val = 2;
return ret_val;
L100:
/* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) */
ret_val = (integer) ((real) min(*n1,*n2) * 1.6f);
return ret_val;
L110:
/* ISPEC = 7: number of processors (not used) */
ret_val = 1;
return ret_val;
L120:
/* ISPEC = 8: crossover point for multishift (used by xHSEQR) */
ret_val = 50;
return ret_val;
L130:
/* ISPEC = 9: maximum size of the subproblems at the bottom of the */
/* computation tree in the divide-and-conquer algorithm */
/* (used by xGELSD and xGESDD) */
ret_val = 25;
return ret_val;
L140:
/* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap */
/* ILAENV = 0 */
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__1, &c_b163, &c_b164);
}
return ret_val;
L150:
/* ISPEC = 11: infinity arithmetic can be trusted not to trap */
/* ILAENV = 0 */
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__0, &c_b163, &c_b164);
}
return ret_val;
L160:
/* 12 <= ISPEC <= 16: xHSEQR or one of its subroutines. */
ret_val = iparmq_(ispec, name__, opts, n1, n2, n3, n4)
;
return ret_val;
/* End of ILAENV */
} /* ilaenv_ */

191
3rdparty/lapack/ilaenv_custom.c vendored Normal file
View File

@ -0,0 +1,191 @@
/* ilaenv.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
#include "string.h"
/* Table of constant values */
static integer c__1 = 1;
static real c_b163 = 0.f;
static real c_b164 = 1.f;
static integer c__0 = 0;
integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4)
{
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* January 2007 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ILAENV is called from the LAPACK routines to choose problem-dependent */
/* parameters for the local environment. See ISPEC for a description of */
/* the parameters. */
/* ILAENV returns an INTEGER */
/* if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC */
/* if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value. */
/* This version provides a set of parameters which should give good, */
/* but not optimal, performance on many of the currently available */
/* computers. Users are encouraged to modify this subroutine to set */
/* the tuning parameters for their particular machine using the option */
/* and problem size information in the arguments. */
/* This routine will not function correctly if it is converted to all */
/* lower case. Converting it to all upper case is allowed. */
/* Arguments */
/* ========= */
/* ISPEC (input) INTEGER */
/* Specifies the parameter to be returned as the value of */
/* ILAENV. */
/* = 1: the optimal blocksize; if this value is 1, an unblocked */
/* algorithm will give the best performance. */
/* = 2: the minimum block size for which the block routine */
/* should be used; if the usable block size is less than */
/* this value, an unblocked routine should be used. */
/* = 3: the crossover point (in a block routine, for N less */
/* than this value, an unblocked routine should be used) */
/* = 4: the number of shifts, used in the nonsymmetric */
/* eigenvalue routines (DEPRECATED) */
/* = 5: the minimum column dimension for blocking to be used; */
/* rectangular blocks must have dimension at least k by m, */
/* where k is given by ILAENV(2,...) and m by ILAENV(5,...) */
/* = 6: the crossover point for the SVD (when reducing an m by n */
/* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds */
/* this value, a QR factorization is used first to reduce */
/* the matrix to a triangular form.) */
/* = 7: the number of processors */
/* = 8: the crossover point for the multishift QR method */
/* for nonsymmetric eigenvalue problems (DEPRECATED) */
/* = 9: maximum size of the subproblems at the bottom of the */
/* computation tree in the divide-and-conquer algorithm */
/* (used by xGELSD and xGESDD) */
/* =10: ieee NaN arithmetic can be trusted not to trap */
/* =11: infinity arithmetic can be trusted not to trap */
/* 12 <= ISPEC <= 16: */
/* xHSEQR or one of its subroutines, */
/* see IPARMQ for detailed explanation */
/* NAME (input) CHARACTER*(*) */
/* The name of the calling subroutine, in either upper case or */
/* lower case. */
/* OPTS (input) CHARACTER*(*) */
/* The character options to the subroutine NAME, concatenated */
/* into a single character string. For example, UPLO = 'U', */
/* TRANS = 'T', and DIAG = 'N' for a triangular routine would */
/* be specified as OPTS = 'UTN'. */
/* N1 (input) INTEGER */
/* N2 (input) INTEGER */
/* N3 (input) INTEGER */
/* N4 (input) INTEGER */
/* Problem dimensions for the subroutine NAME; these may not all */
/* be required. */
/* Further Details */
/* =============== */
/* The following conventions have been used when calling ILAENV from the */
/* LAPACK routines: */
/* 1) OPTS is a concatenation of all of the character options to */
/* subroutine NAME, in the same order that they appear in the */
/* argument list for NAME, even if they are not used in determining */
/* the value of the parameter specified by ISPEC. */
/* 2) The problem dimensions N1, N2, N3, N4 are specified in the order */
/* that they appear in the argument list for NAME. N1 is used */
/* first, N2 second, and so on, and unused problem dimensions are */
/* passed a value of -1. */
/* 3) The parameter value returned by ILAENV is checked for validity in */
/* the calling subroutine. For example, ILAENV is used to retrieve */
/* the optimal blocksize for STRTRI as follows: */
/* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) */
/* IF( NB.LE.1 ) NB = MAX( 1, N ) */
/* ===================================================================== */
/* .. Local Scalars .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Executable Statements .. */
switch (*ispec) {
case 1:
/* ISPEC = 1: block size */
/* In these examples, separate code is provided for setting NB for */
/* real and complex. We assume that NB will take the same value in */
/* single or double precision. */
return 1;
case 2:
/* ISPEC = 2: minimum block size */
return 2;
case 3:
/* ISPEC = 3: crossover point */
return 3;
case 4:
/* ISPEC = 4: number of shifts (used by xHSEQR) */
return 6;
case 5:
/* ISPEC = 5: minimum column dimension (not used) */
return 2;
case 6:
/* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) */
return (integer) ((real) min(*n1,*n2) * 1.6f);
case 7:
/* ISPEC = 7: number of processors (not used) */
return 1;
case 8:
/* ISPEC = 8: crossover point for multishift (used by xHSEQR) */
return 50;
case 9:
/* ISPEC = 9: maximum size of the subproblems at the bottom of the */
/* computation tree in the divide-and-conquer algorithm */
/* (used by xGELSD and xGESDD) */
return 25;
case 10:
/* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap */
return ieeeck_(&c__1, &c_b163, &c_b164);
case 11:
/* ISPEC = 11: infinity arithmetic can be trusted not to trap */
return ieeeck_(&c__0, &c_b163, &c_b164);
case 12:
case 13:
case 14:
case 15:
case 16:
/* 12 <= ISPEC <= 16: xHSEQR or one of its subroutines. */
return iparmq_(ispec, name__, opts, n1, n2, n3, n4);
default:
break;
}
/* Invalid value for ISPEC */
return -1;
/* End of ILAENV */
} /* ilaenv_ */

312
3rdparty/lapack/sgemv.c vendored
View File

@ -1,312 +0,0 @@
/* sgemv.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
/* Subroutine */ int sgemv_(char *trans, integer *m, integer *n, real *alpha,
real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
integer *incy)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
integer i__, j, ix, iy, jx, jy, kx, ky, info;
real temp;
integer lenx, leny;
extern logical lsame_(char *, char *);
extern /* Subroutine */ int xerbla_(char *, integer *);
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* SGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Arguments */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - REAL . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - REAL array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared */
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - REAL array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - REAL . */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - REAL array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the */
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- Written on 22-October-1986. */
/* Jack Dongarra, Argonne National Lab. */
/* Jeremy Du Croz, Nag Central Office. */
/* Sven Hammarling, Nag Central Office. */
/* Richard Hanson, Sandia National Labs. */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--x;
--y;
/* Function Body */
info = 0;
if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
) {
info = 1;
} else if (*m < 0) {
info = 2;
} else if (*n < 0) {
info = 3;
} else if (*lda < max(1,*m)) {
info = 6;
} else if (*incx == 0) {
info = 8;
} else if (*incy == 0) {
info = 11;
}
if (info != 0) {
xerbla_("SGEMV ", &info);
return 0;
}
/* Quick return if possible. */
if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
return 0;
}
/* Set LENX and LENY, the lengths of the vectors x and y, and set */
/* up the start points in X and Y. */
if (lsame_(trans, "N")) {
lenx = *n;
leny = *m;
} else {
lenx = *m;
leny = *n;
}
if (*incx > 0) {
kx = 1;
} else {
kx = 1 - (lenx - 1) * *incx;
}
if (*incy > 0) {
ky = 1;
} else {
ky = 1 - (leny - 1) * *incy;
}
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
/* First form y := beta*y. */
if (*beta != 1.f) {
if (*incy == 1) {
if (*beta == 0.f) {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[i__] = 0.f;
/* L10: */
}
} else {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[i__] = *beta * y[i__];
/* L20: */
}
}
} else {
iy = ky;
if (*beta == 0.f) {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[iy] = 0.f;
iy += *incy;
/* L30: */
}
} else {
i__1 = leny;
for (i__ = 1; i__ <= i__1; ++i__) {
y[iy] = *beta * y[iy];
iy += *incy;
/* L40: */
}
}
}
}
if (*alpha == 0.f) {
return 0;
}
if (lsame_(trans, "N")) {
/* Form y := alpha*A*x + y. */
jx = kx;
if (*incy == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (x[jx] != 0.f) {
temp = *alpha * x[jx];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
y[i__] += temp * a[i__ + j * a_dim1];
/* L50: */
}
}
jx += *incx;
/* L60: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (x[jx] != 0.f) {
temp = *alpha * x[jx];
iy = ky;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
y[iy] += temp * a[i__ + j * a_dim1];
iy += *incy;
/* L70: */
}
}
jx += *incx;
/* L80: */
}
}
} else {
/* Form y := alpha*A'*x + y. */
jy = ky;
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp = 0.f;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp += a[i__ + j * a_dim1] * x[i__];
/* L90: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L100: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
temp = 0.f;
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
temp += a[i__ + j * a_dim1] * x[ix];
ix += *incx;
/* L110: */
}
y[jy] += *alpha * temp;
jy += *incy;
/* L120: */
}
}
}
return 0;
/* End of SGEMV . */
} /* sgemv_ */

204
3rdparty/lapack/sgemv_custom.c vendored Normal file
View File

@ -0,0 +1,204 @@
#include "clapack.h"
#include <assert.h>
/* Subroutine */ int sgemv_(char *_trans, integer *_m, integer *_n, real *_alpha,
real *a, integer *_lda, real *x, integer *_incx, real *_beta, real *y,
integer *_incy)
{
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* SGEMV performs one of the matrix-vector operations */
/* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
/* where alpha and beta are scalars, x and y are vectors and A is an */
/* m by n matrix. */
/* Arguments */
/* ========== */
/* TRANS - CHARACTER*1. */
/* On entry, TRANS specifies the operation to be performed as */
/* follows: */
/* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
/* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
/* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
/* Unchanged on exit. */
/* M - INTEGER. */
/* On entry, M specifies the number of rows of the matrix A. */
/* M must be at least zero. */
/* Unchanged on exit. */
/* N - INTEGER. */
/* On entry, N specifies the number of columns of the matrix A. */
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - REAL . */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* A - REAL array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. */
/* Unchanged on exit. */
/* LDA - INTEGER. */
/* On entry, LDA specifies the first dimension of A as declared */
/* in the calling (sub) program. LDA must be at least */
/* max( 1, m ). */
/* Unchanged on exit. */
/* X - REAL array of DIMENSION at least */
/* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
/* Before entry, the incremented array X must contain the */
/* vector x. */
/* Unchanged on exit. */
/* INCX - INTEGER. */
/* On entry, INCX specifies the increment for the elements of */
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* BETA - REAL . */
/* On entry, BETA specifies the scalar beta. When BETA is */
/* supplied as zero then Y need not be set on input. */
/* Unchanged on exit. */
/* Y - REAL array of DIMENSION at least */
/* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
/* and at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
/* Before entry with BETA non-zero, the incremented array Y */
/* must contain the vector y. On exit, Y is overwritten by the */
/* updated vector y. */
/* INCY - INTEGER. */
/* On entry, INCY specifies the increment for the elements of */
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* Level 2 Blas routine. */
/* -- Written on 22-October-1986. */
/* Jack Dongarra, Argonne National Lab. */
/* Jeremy Du Croz, Nag Central Office. */
/* Sven Hammarling, Nag Central Office. */
/* Richard Hanson, Sandia National Labs. */
/* Test the input parameters. */
/* Function Body */
char trans = lapack_toupper(_trans[0]);
integer i, j, m = *_m, n = *_n, lda = *_lda, incx = *_incx, incy = *_incy;
integer leny = trans == 'N' ? m : n, lenx = trans == 'N' ? n : m;
real alpha = *_alpha, beta = *_beta;
integer info = 0;
if (trans != 'N' && trans != 'T' && trans != 'C')
info = 1;
else if (m < 0)
info = 2;
else if (n < 0)
info = 3;
else if (lda < max(1,m))
info = 6;
else if (incx == 0)
info = 8;
else if (incy == 0)
info = 11;
if (info != 0)
{
xerbla_("SGEMV ", &info);
return 0;
}
if( incy < 0 )
y -= incy*(leny - 1);
if( incx < 0 )
x -= incx*(lenx - 1);
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if( beta != 1.f )
{
if( incy == 1 )
{
if( beta == 0.f )
for( i = 0; i < leny; i++ )
y[i] = 0.f;
else
for( i = 0; i < leny; i++ )
y[i] *= beta;
}
else
{
if( beta == 0.f )
for( i = 0; i < leny; i++ )
y[i*incy] = 0.f;
else
for( i = 0; i < leny; i++ )
y[i*incy] *= beta;
}
}
if( alpha == 0.f )
;
else if( trans == 'N' )
{
for( i = 0; i < n; i++, a += lda )
{
real s = x[i*incx];
if( s == 0.f )
continue;
s *= alpha;
for( j = 0; j <= m - 4; j += 4 )
{
real t0 = y[j] + s*a[j];
real t1 = y[j+1] + s*a[j+1];
y[j] = t0; y[j+1] = t1;
t0 = y[j+2] + s*a[j+2];
t1 = y[j+3] + s*a[j+3];
y[j+2] = t0; y[j+3] = t1;
}
for( ; j < m; j++ )
y[j] += s*a[j];
}
}
else
{
for( i = 0; i < n; i++, a += lda )
{
real s = 0;
for( j = 0; j <= m - 4; j += 4 )
s += x[j]*a[j] + x[j+1]*a[j+1] + x[j+2]*a[j+2] + x[j+3]*a[j+3];
for( ; j < m; j++ )
s += x[j]*a[j];
y[i*incy] += alpha*s;
}
}
return 0;
/* End of SGEMV . */
} /* sgemv_ */

161
3rdparty/lapack/sger.c → 3rdparty/lapack/sger_custom.c vendored
View File

@ -1,28 +1,9 @@
/* sger.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
/* Subroutine */ int sger_(integer *m, integer *n, real *alpha, real *x,
integer *incx, real *y, integer *incy, real *a, integer *lda)
/* Subroutine */ int sger_(integer *_m, integer *_n, real *_alpha,
real *x, integer *_incx, real *y, integer *_incy,
real *a, integer *_lda)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
/* Local variables */
integer i__, j, ix, jy, kx, info;
real temp;
extern /* Subroutine */ int xerbla_(char *, integer *);
/* .. Scalar Arguments .. */
/* .. */
@ -52,11 +33,11 @@
/* N must be at least zero. */
/* Unchanged on exit. */
/* ALPHA - REAL . */
/* ALPHA - SINGLE PRECISION. */
/* On entry, ALPHA specifies the scalar alpha. */
/* Unchanged on exit. */
/* X - REAL array of dimension at least */
/* X - SINGLE PRECISION array of dimension at least */
/* ( 1 + ( m - 1 )*abs( INCX ) ). */
/* Before entry, the incremented array X must contain the m */
/* element vector x. */
@ -67,7 +48,7 @@
/* X. INCX must not be zero. */
/* Unchanged on exit. */
/* Y - REAL array of dimension at least */
/* Y - SINGLE PRECISION array of dimension at least */
/* ( 1 + ( n - 1 )*abs( INCY ) ). */
/* Before entry, the incremented array Y must contain the n */
/* element vector y. */
@ -78,7 +59,7 @@
/* Y. INCY must not be zero. */
/* Unchanged on exit. */
/* A - REAL array of DIMENSION ( LDA, n ). */
/* A - SINGLE PRECISION array of DIMENSION ( LDA, n ). */
/* Before entry, the leading m by n part of the array A must */
/* contain the matrix of coefficients. On exit, A is */
/* overwritten by the updated matrix. */
@ -110,80 +91,70 @@
/* Test the input parameters. */
/* Parameter adjustments */
--x;
--y;
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
info = 0;
if (*m < 0) {
info = 1;
} else if (*n < 0) {
info = 2;
} else if (*incx == 0) {
info = 5;
} else if (*incy == 0) {
info = 7;
} else if (*lda < max(1,*m)) {
info = 9;
}
if (info != 0) {
xerbla_("SGER ", &info);
return 0;
integer i, j, m = *_m, n = *_n, incx = *_incx, incy = *_incy, lda = *_lda;
real alpha = *_alpha;
integer info = 0;
if (m < 0)
info = 1;
else if (n < 0)
info = 2;
else if (incx == 0)
info = 5;
else if (incy == 0)
info = 7;
else if (lda < max(1,m))
info = 9;
if (info != 0)
{
xerbla_("SGER ", &info);
return 0;
}
/* Quick return if possible. */
if (incx < 0)
x -= (m-1)*incx;
if (incy < 0)
y -= (n-1)*incy;
if (*m == 0 || *n == 0 || *alpha == 0.f) {
return 0;
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if( alpha == 0 )
;
else if( incx == 1 )
{
for( j = 0; j < n; j++, a += lda )
{
real s = y[j*incy];
if( s == 0 )
continue;
s *= alpha;
for( i = 0; i <= m - 2; i += 2 )
{
real t0 = a[i] + x[i]*s;
real t1 = a[i+1] + x[i+1]*s;
a[i] = t0; a[i+1] = t1;
}
for( ; i < m; i++ )
a[i] += x[i]*s;
}
}
/* Start the operations. In this version the elements of A are */
/* accessed sequentially with one pass through A. */
if (*incy > 0) {
jy = 1;
} else {
jy = 1 - (*n - 1) * *incy;
}
if (*incx == 1) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (y[jy] != 0.f) {
temp = *alpha * y[jy];
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
a[i__ + j * a_dim1] += x[i__] * temp;
/* L10: */
}
}
jy += *incy;
/* L20: */
}
} else {
if (*incx > 0) {
kx = 1;
} else {
kx = 1 - (*m - 1) * *incx;
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (y[jy] != 0.f) {
temp = *alpha * y[jy];
ix = kx;
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
a[i__ + j * a_dim1] += x[ix] * temp;
ix += *incx;
/* L30: */
}
}
jy += *incy;
/* L40: */
}
else
{
for( j = 0; j < n; j++, a += lda )
{
real s = y[j*incy];
if( s == 0 )
continue;
s *= alpha;
for( i = 0; i < m; i++ )
a[i] += x[i*incx]*s;
}
}
return 0;

1045
3rdparty/lapack/slamch.c vendored
View File

File diff suppressed because it is too large Load Diff

88
3rdparty/lapack/slamch_custom.c vendored Normal file
View File

@ -0,0 +1,88 @@
#include "clapack.h"
#include <float.h>
#include <stdio.h>
/* *********************************************************************** */
doublereal slamc3_(real *a, real *b)
{
/* System generated locals */
real ret_val;
/* -- LAPACK auxiliary routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* SLAMC3 is intended to force A and B to be stored prior to doing */
/* the addition of A and B , for use in situations where optimizers */
/* might hold one of these in a register. */
/* Arguments */
/* ========= */
/* A (input) REAL */
/* B (input) REAL */
/* The values A and B. */
/* ===================================================================== */
/* .. Executable Statements .. */
ret_val = *a + *b;
return ret_val;
/* End of SLAMC3 */
} /* slamc3_ */
const unsigned char lapack_toupper_tab[] =
{
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
90, 91, 92, 93, 94, 95, 96, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 123, 124, 125, 126, 127, 128, 129, 130, 131,
132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167,
168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185,
186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203,
204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221,
222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255
};
/* simpler version of dlamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
#ifndef FLT_DIGITS
#define FLT_DIGITS 24
#endif
const unsigned char lapack_lamch_tab[] =
{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9, 0, 10, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9,
0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};
const doublereal lapack_slamch_tab[] =
{
0, FLT_RADIX, FLT_EPSILON, FLT_MAX_EXP, FLT_MIN_EXP, FLT_DIGITS, FLT_MAX,
FLT_EPSILON*FLT_RADIX, 1, FLT_MIN*(1 + FLT_EPSILON), FLT_MIN
};

41
3rdparty/lapack/slartg.c → 3rdparty/lapack/slartg_custom.c vendored
View File

@ -1,15 +1,3 @@
/* slartg.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "clapack.h"
@ -19,17 +7,13 @@
integer i__1;
real r__1, r__2;
/* Builtin functions */
double log(doublereal), pow_ri(real *, integer *), sqrt(doublereal);
/* Local variables */
integer i__;
real f1, g1, eps, scale;
integer count;
real safmn2, safmx2;
extern doublereal slamch_(char *);
real safmin;
static real safmn2, safmx2;
static real safmin;
static logical FIRST = TRUE_;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
@ -96,15 +80,16 @@
/* .. */
/* .. Executable Statements .. */
/* IF( FIRST ) THEN */
safmin = slamch_("S");
eps = slamch_("E");
r__1 = slamch_("B");
i__1 = (integer) (log(safmin / eps) / log(slamch_("B")) / 2.f);
safmn2 = pow_ri(&r__1, &i__1);
safmx2 = 1.f / safmn2;
/* FIRST = .FALSE. */
/* END IF */
if(FIRST)
{
safmin = slamch_("S");
eps = slamch_("E");
r__1 = slamch_("B");
i__1 = (integer) (log(safmin / eps) / log(slamch_("B")) / 2.f);
safmn2 = pow_ri(&r__1, &i__1);
safmx2 = 1.f / safmn2;
FIRST = FALSE_;
}
if (*g == 0.f) {
*cs = 1.f;
*sn = 0.f;

0
3rdparty/lapack/slasr.c → 3rdparty/lapack/slasr_custom.c vendored
View File