Here is an example: please try it and see what you get.
Included are files svdd-demo.c and array.h and test.dat
To compile, use
% gcc svdd-demo.c -ansi -Wall -pedantic -lm -llapack
% a.out
reading in matrix A from test.dat:
2.117079 2.992091
2.106365 1.970807
2.131384 3.130390
2.344235 2.667219
1.838804 3.243629
** On entry to DGESDD parameter number 12 had an illegal value
STOP 0
If you change m to <=4 it works fine. So there is an error for m "substantially
larger" than n (=2 here). Works fine on SuSE and RHEL3/4.
Thanks!
__________________________________________________
Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around
http://mail.yahoo.com #include <stdio.h>
#include <stdlib.h>
#include "array.h"
#define MIN(i,j) ((i)<(j) ? (i) : (j))
#define MAX(i,j) ((i)>(j) ? (i) : (j))
void dgesdd_( char *jobz,
int *m,
int *n,
double *a,
int *lda,
double *s,
double *u,
int *ldu,
double *vt,
int *ldvt,
double *work,
int *lwork,
int *iwork,
int *info);
int svdd(double **a, int m, int n, double *d, double **u, double **v)
{
double *A, *U, *VT;
int lwork = 3*MIN(m,n)*MIN(m,n)
+ MAX(MAX(m,n),4*MIN(m,n)*MIN(m,n)+4*MIN(m,n));
int liwork = 8*MIN(m,n);
char jobz = 'A';
double *work;
int *iwork;
int i, j, k, info;
MAKE_VECTOR(A, m*n);
for (j=0, k=0; j<n; j++)
for (i=0; i<m; i++)
A[k++] = a[i][j];
MAKE_VECTOR(U, m*m);
MAKE_VECTOR(VT, n*n);
MAKE_VECTOR(work, lwork);
MAKE_VECTOR(iwork, liwork);
dgesdd_(&jobz, &m, &n, A, &m, d, U, &m, VT, &n,
work, &lwork, iwork, &info);
for (j=0, k=0; j<m; j++)
for (i=0; i<m; i++)
u[i][j] = U[k++];
/* VT, as calculated by dgesdd_(), is the transpose of the right
* multiplier. Here we undo the transpose so that the matrix
* v[][] returned by this function is not transposed anymore.
*/
for (j=0, k=0; j<n; j++)
for (i=0; i<n; i++)
v[j][i] = VT[k++];
FREE_VECTOR(A);
FREE_VECTOR(U);
FREE_VECTOR(VT);
FREE_VECTOR(work);
FREE_VECTOR(iwork);
return info;
}
int main(void)
{
FILE *fp;
double **a, **u, **v, *d;
int m=5, n=2;
int ld = MIN(m,n);
int i, j;
int result;
fp = fopen("test.dat", "r");
if (fp == NULL) {
fprintf(stderr, "Error opening file test.dat\n");
return EXIT_FAILURE;
}
MAKE_MATRIX(a, m, n);
MAKE_MATRIX(u, m, m);
MAKE_MATRIX(v, n, n);
MAKE_VECTOR(d, ld);
puts("reading in matrix A from test.dat:");
for(i=0; i<m; i++) {
for(j=0; j<n; j++) {
fscanf(fp, "%lf", &a[i][j]);
printf("%12.6f ", a[i][j]);
}
putchar('\n');
}
putchar('\n');
putchar('\n');
fclose(fp);
result = svdd(a, m, n, d, u, v);
printf("svdd returned %d\n", result);
if (result!=0) {
fprintf(stderr, "exiting!\n");
return EXIT_FAILURE;
}
putchar('\n');
printf("singular values are:\n");
for (i=0; i<ld; i++)
printf("%12.6f ", d[i]);
putchar('\n');
putchar('\n');
puts("The left multiplier matrix, u, is:");
for(i=0; i<m; i++) {
for(j=0; j<m; j++)
printf("%12.6f ", u[i][j]);
putchar('\n');
}
putchar('\n');
puts("The right multiplier matrix, v, is:");
for(i=0; i<n; i++) {
for(j=0; j<n; j++)
printf("%12.6f ", v[i][j]);
putchar('\n');
}
putchar('\n');
FREE_MATRIX(a);
FREE_MATRIX(u);
FREE_MATRIX(v);
FREE_VECTOR(d);
return EXIT_SUCCESS;
}
/* array.h
o This file defines the following macros:
MAKE_1ARRAY(a,n) make a 1D array of length "n"
MAKE_2ARRAY(a,m,n) make a 2D array of dimensions "m x n"
MAKE_3ARRAY(a,l,m,n) make a 3D array of dimensions "l x m x n"
MAKE_4ARRAY(a,k,l,m,n) make a 4D array of dimensions "k x l x m x n"
FREE_1ARRAY(a) free memory allocated by MAKE_1ARRAY()
FREE_2ARRAY(a) free memory allocated by MAKE_2ARRAY()
FREE_3ARRAY(a) free memory allocated by MAKE_3ARRAY()
FREE_4ARRAY(a) free memory allocated by MAKE_4ARRAY()
o Additionally, it defines the following convenience macros as synonyms:
MAKE_VECTOR(a,n) same as MAKE_1ARRAY(a,n)
FREE_VECTOR(a) same as FREE_1ARRAY(a)
MAKE_MATRIX(a,m,n) same as MAKE_2ARRAY(a,m,n)
FREE_MATRIX(a) same as FREE_2ARRAY(a)
o Additionally, it declares and uses the identifiers:
ARRAY_H2RESERVED
ARRAY_H3RESERVED
ARRAY_H4RESERVED
within local blocks within the macros. THESE IDENTIFIERS
ARE RESERVED. The user should not use identifiers with this
names within the scope of these macros.
o vector/matrix/array elements can be of _any_ type.
o If malloc() fails during the execution of any of the MAKE_* macros:
. An "out of memory" message is printed to stderr.
. The macro's first argument is set to NULL.
o After a call to any of the FREE_*() macros, the macro's argument
is set to NULL.
o The FREE_*() macros can be applied to previously FREE_*()-ed
arrays with no ill effect.
o Note that macro arguments are evaluated more than once.
o Customization
When malloc() returns NULL, MAKE_1ARRAY prints a message on stderr
and the program continues. The user can alter this behavior by
redefining the MAKE_1ARRAY macro. For instance, to call exit()
whenever malloc() fails, the user can do:
#undef MAKE_1ARRAY
#define MAKE_1ARRAY(a,n) do { \
(a) = malloc((n) * sizeof *(a)); \
if ((a)==NULL) { \
fprintf(stderr, "*** in file %s, function %s(), line %d: " \
"out of memory!\n", __FILE__, __func__, __LINE__); \
exit(EXIT_FAILURE); \
} \
} while (0)
Since only MAKE_1ARRAY calls malloc() explicitly, this change affects
the behavior of not only MAKE_1ARRAY but all other MAKE_* macros as well.
---- SAMPLE USAGE -------------
#include "array.h"
int main(void)
{
float ***a; // can use any other type instead of "float"
size_t p=3, q=4, r=5; // will make a 3x4x5 3-D array of float
MAKE_3ARRAY(a, p, q, r);
if (a==NULL)
return EXIT_FAILURE;
a[2][0][1] = 3.14;
printf("%g \n", a[2][0][1]);
FREE_3ARRAY(a);
return EXIT_SUCCESS;
}
---- END OF SAMPLE USAGE -------
Author: Rouben Rostamian, <rostamian@xxxxxxxx>
June 2000
Revised Jul 2000
Revised Sep 2002
Revised Jun 2003 - macros don't call exit anymore
changed macro names to all-caps
Revised Aug 2003 - changed reserved names to all caps
Revised Dec 2003 - changed array index types to size_t (were int)
- added an "out of memory" message, printed to stderr
- most FREE* macros now come before MAKE* macros,
for possible improved efficiency in preprocessing
*/
#ifndef H_ARRAY_
#define H_ARRAY_
#include <stdio.h>
#include <stdlib.h>
/* ---------- 1D arrays ---------------------- */
#define MAKE_1ARRAY(a,n) do { \
(a) = malloc((n) * sizeof *(a)); \
if ((a)==NULL) \
fprintf(stderr, "*** in file %s, function %s(), line %d: " \
"out of memory!\n", __FILE__, __func__, __LINE__); \
} while (0)
#define FREE_1ARRAY(a) do { \
free(a); \
a = NULL; \
} while (0)
/* ---------- 2D arrays ---------------------- */
#define FREE_2ARRAY(a) do { \
size_t ARRAY_H2RESERVED; \
if (a==NULL) \
break; \
for (ARRAY_H2RESERVED=0; (a)[ARRAY_H2RESERVED]!=NULL; ARRAY_H2RESERVED++)\
FREE_1ARRAY((a)[ARRAY_H2RESERVED]); \
FREE_1ARRAY(a); \
} while (0)
/* note: parenthesize first arg because it may be given as `*a' */
#define MAKE_2ARRAY(a,m,n) do { \
size_t ARRAY_H2RESERVED; \
MAKE_1ARRAY(a,(m)+1); \
if (a==NULL) \
break; \
(a)[m] = NULL; \
for (ARRAY_H2RESERVED=0; ARRAY_H2RESERVED<(m); ARRAY_H2RESERVED++) { \
MAKE_1ARRAY((a)[ARRAY_H2RESERVED],(n)); \
if ((a)[ARRAY_H2RESERVED]==NULL) { \
FREE_2ARRAY(a); \
break; \
} \
} \
} while (0)
/* ---------- 3D arrays ---------------------- */
#define FREE_3ARRAY(a) do { \
size_t ARRAY_H3RESERVED; \
if (a==NULL) \
break; \
for (ARRAY_H3RESERVED=0; (a)[ARRAY_H3RESERVED]!=NULL; ARRAY_H3RESERVED++)\
FREE_2ARRAY((a)[ARRAY_H3RESERVED]); \
FREE_1ARRAY(a); \
} while (0)
#define MAKE_3ARRAY(a,p,q,r) do { \
size_t ARRAY_H3RESERVED; \
MAKE_1ARRAY(a,(p)+1); \
if (a==NULL) \
break; \
(a)[p] = NULL; \
for (ARRAY_H3RESERVED=0; ARRAY_H3RESERVED<(p); ARRAY_H3RESERVED++) { \
MAKE_2ARRAY((a)[ARRAY_H3RESERVED],(q),(r)); \
if ((a)[ARRAY_H3RESERVED]==NULL) { \
FREE_3ARRAY(a); \
break; \
} \
} \
} while (0)
/* ---------- 3D arrays ---------------------- */
#define FREE_4ARRAY(a) do { \
size_t ARRAY_H4RESERVED; \
if (a==NULL) \
break; \
for (ARRAY_H4RESERVED=0; (a)[ARRAY_H4RESERVED]!=NULL; ARRAY_H4RESERVED++)\
FREE_3ARRAY((a)[ARRAY_H4RESERVED]); \
FREE_1ARRAY(a); \
} while (0)
#define MAKE_4ARRAY(a,p,q,r,s) do { \
size_t ARRAY_H4RESERVED; \
MAKE_1ARRAY(a,(p)+1); \
if (a==NULL) \
break; \
(a)[p] = NULL; \
for (ARRAY_H4RESERVED=0; ARRAY_H4RESERVED<(p); ARRAY_H4RESERVED++) { \
MAKE_3ARRAY((a)[ARRAY_H4RESERVED],(q),(r),(s)); \
if ((a)[ARRAY_H4RESERVED]==NULL) { \
FREE_4ARRAY(a); \
break; \
} \
} \
} while (0)
/* ---------- synonyms ---------------------- */
#define MAKE_VECTOR(a,n) MAKE_1ARRAY(a,n)
#define MAKE_MATRIX(a,m,n) MAKE_2ARRAY(a,m,n)
#define FREE_VECTOR(a) FREE_1ARRAY(a)
#define FREE_MATRIX(a) FREE_2ARRAY(a)
#endif /* H_ARRAY_ */
Attachment:
test.dat
Description: 2998025909-test.dat
