Package: lapack3-dev
Version: 3.0.20000531a-6
Severity: normal
Hi, I'm a little hesitant to file this bug report, since lapack is
such a well-tested package.
BACKGROUND:
The lapack function dgesdd_(), used for computing the SVD of a general
matrix, requires an externally supplied workspace. One way to know
how big this workspace must be is to call dgesdd_() with parameter 12,
LWORK, set to -1. In this case the usual SVD commputation is not
carried out. Instead, the optimal workspace size is returned through
parameter 11, WORK. A correctly sized workspace can then be allocated
and used in subsequent calls to dgesdd_().
PROBLEM:
After upgrading to lapack3-dev, I'm observing that the initial call to
dgesdd_() no longer returns a valid workspace size. Attempting to use
the returned value results in an error and, apparently, a call to
exit().
System info is at the end of this email.
EXAMPLE:
I've included a code snippet below which exercises the bug. You can
compile it by saving as dgesddExample.c:
> gcc dgesddExample.c -o dgesddExample -llapack
You can then run it with no arguments, in which case the erronious
return value from dgesdd_() will be corrected:
> ./dgesddExample
Warning: corrected LWORK from 1 to 27.
Both calls to dgesdd() completed successfully.
Or you can run it with an argument, in which case the erronious return
value from dgesdd_() will not be corrected, and the second call to
dgesdd_() will raise an error:
> ./dgesddExample crash
** On entry to DGEBRD parameter number 10 had an illegal value
Thanks for your help and hard work,
David
--------------------- Code starts here ---------------------
/**
* This file exercises a bug in the lapack routine dgesdd(), in which
* the reported workspace requirement is insufficient.
*/
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
/**
* This is a declaration for the LAPACK routine dgesdd(), which
* computes the singular value decomposition of a matrix using a
* fast divide & conquer algorithm.
*/
void dgesdd_(char* JOBZ, long int* M, long int* N, double* A,
long int* LDA, double* S, double* U, long int* LDU,
double* VT, long int* LDVT, double* WORK, long int* LWORK,
long int* IWORK, long int* INFO);
/**
* The main() function here simply computes the SVD of a 4x3 matrix:
*
*
* |1.0, 5.0. 9.0|
* A = |2.0, 6.0, 10.0|
* |3.0, 7.0, 11.0|
* |4.0, 8.0, 12.0|
*
* When called with no arguments, the program checks the LWORK size
* returned from dgesdd_(), and corrects it. If an input argument is
* supplied, then the value of LWORK is not corrected, and the program
* crashes.
*/
int main(int argc, char* argv[])
{
const int A_ROWS = 4;
const int A_COLUMNS = 3;
const int A_SIZE = A_ROWS * A_COLUMNS;
const int MINIMUM_DIMENSION = A_ROWS > A_COLUMNS ? A_COLUMNS : A_ROWS;
const int MAXIMUM_DIMENSION = A_ROWS > A_COLUMNS ? A_ROWS : A_COLUMNS;
/* We define a couple of constants to find the minimum acceptable */
/* size of the WORK array which will be passed to dgesdd() */
const int tempVal0 = (MAXIMUM_DIMENSION > 6 * MINIMUM_DIMENSION ?
MAXIMUM_DIMENSION : 6 * MINIMUM_DIMENSION);
const long int minimumLWORK = 3 * MINIMUM_DIMENSION + tempVal0;
/* Arguments for lapack call. */
char JOBZ = 'N'; /* Compute no columns of U and no rows of VT. */
double A[A_SIZE]; /* We'll have a 4x3 matrix. */
double S[MINIMUM_DIMENSION]; /* We'll compute 3 singular values. */
long int IWORK[8 * MINIMUM_DIMENSION]; /* Workspace required by lapack. */
long int M = A_ROWS;
long int N = A_COLUMNS;
long int LDA = M;
double U; /* Since JOBZ == 'N', U is not referenced. */
long int LDU = 1;
double VT; /* Since JOBZ == 'N', VT is not referenced. */
long int LDVT = 1;
double* WORK = (double*)malloc(sizeof(double));
long int LWORK = -1; /* Request info on optimal workspace size. */
long int INFO;
/* Initialize the matrix A. */
int aIndex;
for(aIndex = 0; aIndex < A_SIZE; ++aIndex) {
A[aIndex] = aIndex + 1;
}
/* Call once to request optimal workspace size. */
dgesdd_(&JOBZ, &M, &N, A, &LDA, S, &U, &LDU, &VT, &LDVT,
WORK, &LWORK, IWORK, &INFO);
/* Check for errors. */
if(INFO != 0L) {
fprintf(stderr, "Error: first call to dgesdd_ returns &ld.\n",
INFO);
exit(1);
}
/* Recover optimal workspace size. */
LWORK = (long int)(WORK[0]);
/* Sanity check. Is the returned value of LWORK legal? */
if(LWORK < minimumLWORK) {
if(argc == 1) {
fprintf(stderr, "Warning: changing LWORK from %ld to %ld.\n",
LWORK, minimumLWORK);
LWORK = minimumLWORK;
}
} else {
fprintf(stderr,
"Submitter was mistaken: the returned value for LWORK "
"is perfectly legal");
}
/* Resize workspace. */
free(WORK);
WORK = (double*)malloc(LWORK * sizeof(double));
/* Call dgesdd() again to do the SVD. */
dgesdd_(&JOBZ, &M, &N, A, &LDA, S, &U, &LDU, &VT, &LDVT,
WORK, &LWORK, IWORK, &INFO);
/* Check for errors. */
if(INFO != 0L) {
fprintf(stderr, "Error: secondd call to dgesdd_ returns &ld.\n",
INFO);
exit(3);
}
// Phew! Made it.
fprintf(stderr, "Both calls to dgesdd() completed successfully.\n");
exit(0);
}
--------------------- Code ends here ---------------------
-- System Information:
Debian Release: testing/unstable
APT prefers unstable
APT policy: (500, 'unstable')
Architecture: i386 (i686)
Shell: /bin/sh linked to /bin/bash
Kernel: Linux 2.6.9-1-686-smp
Locale: LANG=C, LC_CTYPE=C (charmap=ANSI_X3.4-1968)
Versions of packages lapack3-dev depends on:
ii g77 4:3.4.4-5 The GNU Fortran 77 compiler
ii lapack3 3.0.20000531a-6 library of linear algebra routines
ii libc6-dev 2.3.5-6 GNU C Library: Development Librari
ii refblas3-dev [libblas-3. 1.2-8 Basic Linear Algebra Subroutines 3
lapack3-dev recommends no packages.
-- no debconf information
--
To UNSUBSCRIBE, email to [EMAIL PROTECTED]
with a subject of "unsubscribe". Trouble? Contact [EMAIL PROTECTED]
