[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]

lapack error (example)



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 umbc edu>
   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


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]