Package: liblapack3gf
Version: 3.0.20000531a-1.3
Severity: important
For a few years now, I have built R with R's own lapack routines. I tried to
switch back yesterday only to get a decisive bug report (#464833) to convince
me otherwise :)
I have spent a good part of today trying to figure out what triggers the bug.
What happens when you call 'example(kappa)' in R is that a few invocations to
svd() occur. The killer is
> svd(h9,0,0)
Error in La.svd(x, nu, nv) :
BLAS/LAPACK routine 'DGEBRD' gave error code -10
where h9 is a 9x9 matrix (see below). The second+third arguments are key,
they basically signal to R and then Lapack to call DGESDD with 'JOBZ' equal
to 'N' -- ie perform a 'small' svd where no U or VT are returned.
What then happens is that svd() calls some other code that eventually sets
two calls to DGESS. The first is meant to provide the size of the main
workspace -- and this fails. This is the bug. We get a return of 1 when we
should get 603.
The bug is not triggered in the
> svd(h9)
case as for full U and V returns, are different path is taken.
I attach an example program below. It runs as follows
[EMAIL PROTECTED]:~/src/progs/C> ./lapack_dgesdd
info 0 tmp 1.000000 != recd 603
1.725883,0.321633,0.031039,0.001979,0.000088,0.000003,0.000000,0.000000,0.000000
Kappa == max(s)/min(s) == 4.861940e+11
[EMAIL PROTECTED]:~/src/progs/C>
and works because I override the faulty output from the first call to
DGESDD. Once you comment that out, it terminates with a similar error as R:
[EMAIL PROTECTED]:~/src/progs/C> ./lapack_dgesdd
info 0 tmp 1.000000 != recd 603
** On entry to DGEBRD parameter number 10 had an illegal value
The program is below. The key overrride line is
lwork = recsize;
which 'fixes' things, ie suppresses the bug. Else it rears its head.
I haven't looked at the Debian sources of lapack to see if our Fortran code
is different from what R ships -- my Fortran knowledge is probably too
rudimentary for this.
Hope this helps, I'd be glad to help with follow-ups.
Dirk
Source code below. I build via
gcc -O3 -Wall -o lapack_dgesdd lapack_dgesdd.c -llapack
-----------------------------------------------------------------------------
/*
* Demonstrate a bug in Lapack's DGESDD
*
* This bug was triggered when R switched from its own Lapack routine to
* Debian's liblapack-dev [ versions: R 2.6.2-1, liblapack 3.0.20000531a-1.3 ]
*
* Simply running
* > example(kappa)
* in R triggers it. The offending code is
* > kappa(h9, exact = TRUE)
* which calls
* if (exact) {
* s <- svd(z, nu = 0, nv = 0)$d
* max(s)/min(s[s > 0])
* }
* and nu=0, nv=0 is key here.
*
* The example creates a 9x9 matrix which is then passed
* via kappa() [ an R function in src/library/base/R/kappa.R ]
* to svd() [ an R function in src/library/base/R/svd.R ]
* to La.svd() [ an R function in src/library/base/R/LAPACK.R ]
* to La_svd [ a C interface routine in src/modules/lapack/Lapack.c)
* to modLa_svd [ a C routine called via a function pointer redirect,
* both src/modules/lapack/Lapack.c ]
* to dgessdd [ a Fortran routine in src/modules/lapack/dlapack1.f ]
*
* Dirk Eddelbuettel, 08 Feb 2009
*
*/
#include <stdio.h>
#include <stdlib.h>
/* defined in /usr/share/R/include/R_ext/Lapack.h, but also see RS.h in the
same directory re underscores */
void dgesdd_(const char *jobz,
const int *m, const int *n,
double *a, const int *lda, double *s,
double *u, const int *ldu,
double *vt, const int *ldvt,
double *work, const int *lwork, int *iwork, int *info);
static int min(a,b) { return a < b ? a : b; };
static int max(a,b) { return a > b ? a : b; };
int main (void) {
char jobu = 'N'; /* 'jobz' job type */
int n = 9; /* 'm' rows */
int p = 9; /* 'n' cols */
double xvals[] = { /* 'a' input/output array */
1.000000000000, 0.500000000000, 0.3333333333333, 0.2500000000000,
0.2000000000000, 0.1666666666667, 0.1428571428571, 0.1250000000000,
0.1111111111111,
0.500000000000, 0.333333333333, 0.2500000000000, 0.2000000000000,
0.1666666666667, 0.1428571428571, 0.1250000000000, 0.1111111111111,
0.1000000000000,
0.333333333333, 0.250000000000, 0.2000000000000, 0.1666666666667,
0.1428571428571, 0.1250000000000, 0.1111111111111, 0.1000000000000,
0.0909090909091,
0.250000000000, 0.200000000000, 0.1666666666667, 0.1428571428571,
0.1250000000000, 0.1111111111111, 0.1000000000000, 0.0909090909091,
0.0833333333333,
0.200000000000, 0.166666666667, 0.1428571428571, 0.1250000000000,
0.1111111111111, 0.1000000000000, 0.0909090909091, 0.0833333333333,
0.0769230769231,
0.166666666667, 0.142857142857, 0.1250000000000, 0.1111111111111,
0.1000000000000, 0.0909090909091, 0.0833333333333, 0.0769230769231,
0.0714285714286,
0.142857142857, 0.125000000000, 0.1111111111111, 0.1000000000000,
0.0909090909091, 0.0833333333333, 0.0769230769231, 0.0714285714286,
0.0666666666667,
0.125000000000, 0.111111111111, 0.1000000000000, 0.0909090909091,
0.0833333333333, 0.0769230769231, 0.0714285714286, 0.0666666666667,
0.0625000000000,
0.111111111111, 0.100000000000, 0.0909090909091, 0.0833333333333,
0.0769230769231, 0.0714285714286, 0.0666666666667, 0.0625000000000,
0.0588235294118
};
int lda = n; /* 'lda' R just uses n again */
double *s, /* 's' output of svd values */
*u, /* 'u' */
*v, /* 'v' */
*work; /* 'tmp */
int ldu = 1; /* 'ldu' */
int ldvt = 1; /* 'ldvt' */
int lwork = 8*9; /* 'lwork' */
int *iwork, info, i;
int recsize;
s = calloc(p, sizeof(double));
u = calloc(ldu, sizeof(double));
v = calloc(ldvt, sizeof(double));
work = calloc(1, sizeof(double));
iwork = calloc(8*n, sizeof(double));
/* ask for optimal size of work array */
/* according to dlapck1.f in R's sources, dgesdd can also provide the
optimal size
* If LWORK = -1 but other input arguments are legal, WORK(1)
* returns the optimal LWORK.
* however this fails here
*/
lwork = -1;
dgesdd_(&jobu, &n, &p, xvals, &lda, s, u, &ldu, v, &ldvt, work, &lwork,
iwork, &info);
if (info != 0) {
fprintf(stderr, "error code %d from Lapack routine '%s'", info,
"dgesdd");
return -1;
}
lwork = work[0];
/* resize is what we should get, but we don't */
recsize = 3*min(n,p)*min(n,p) +
max(max(n,p),4*min(n,p)*min(n,p)+4*min(n,p));
printf("info %d tmp %f != recd %d\n", info, work[0], recsize);
/* so we overrride this */
lwork = recsize;
work = calloc(lwork, sizeof(double));
dgesdd_(&jobu, &n, &p, xvals, &lda, s, u, &ldu, v, &ldvt, work, &lwork,
iwork, &info);
if (info == 0) {
for (i=0; i<9; i++) {
printf("%f%s", s[i], i<9-1 ? "," : "\n");
}
printf("Kappa == max(s)/min(s) == %e\n", s[0]/s[8]);
}
return 0;
}
-----------------------------------------------------------------------------
-- System Information:
Debian Release: lenny/sid
APT prefers testing
APT policy: (500, 'testing')
Architecture: i386 (i686)
Kernel: Linux 2.6.21-2-686 (SMP w/2 CPU cores)
Locale: LANG=en_US.UTF-8, LC_CTYPE=en_US.UTF-8 (charmap=UTF-8)
Shell: /bin/sh linked to /bin/bash
Versions of packages liblapack3gf depends on:
ii debconf [debconf-2.0] 1.5.18 Debian configuration management sy
ii libblas3gf [libblas.so. 1.2-1.3 Basic Linear Algebra Subroutines 3
ii libc6 2.7-6 GNU C Library: Shared libraries
ii libgcc1 1:4.3-20080116-1 GCC support library
ii libgfortran3 4.3-20080116-1 Runtime library for GNU Fortran ap
liblapack3gf recommends no packages.
-- debconf information:
liblapack3gf/minor:
lapack2/ttr:
liblapack3gf/sig:
liblapack3gf/crit:
--
Three out of two people have difficulties with fractions.
--
To UNSUBSCRIBE, email to [EMAIL PROTECTED]
with a subject of "unsubscribe". Trouble? Contact [EMAIL PROTECTED]
