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