This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: gsl interface to LAPACK?


José Luis Gómez Dans writes:
 > 	I guess my questions is whether someone has written a driver for
 > LAPACK, for which you pass a matrix and the software decides how to call
 > the appropriate LAPACK routine.
 > 

I don't know of anyone who has written a driver.

I would like to include something like this in the library, but there
is a licensing issue .... it is not clear to me whether LAPACK has a
license which is compatible with the GPL.

See my posting about the Debian lapack package for details,
http://bugs.debian.org/101683 and on the octave-maintainers list,
http://www.octave.org/mailing-lists/octave-maintainers/2001.

I've attached an example file showing how to call the fortran lapack
with gsl_vector/gsl_matrix arguments.

regards
Brian Gough

----------------------------------------------------------------------
/* With gcc: add an underscore to the end of the function name and
   convert all arguments into pointers.  The Fortran routines use the
   transpose of matrices as they are stored in C, so you can either
   transpose the matrix before input or remember that the output is
   transposed.  I use the liblapack.a library from the Debian lapack-dev
   package, so I did not have to compile lapack myself, although I could
   have done that with g77.  */

#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>

int
main ()
{
  size_t M = 5, N = 5;
  size_t L = GSL_MIN (M, N);

  gsl_matrix_complex *m = gsl_matrix_complex_alloc (M, N);
  gsl_vector *s = gsl_vector_alloc (L);
  gsl_matrix_complex *u = gsl_matrix_complex_alloc (M, N);
  gsl_matrix_complex *v = gsl_matrix_complex_alloc (M, N);

  int lwork = 2 * L + GSL_MAX (M, N);
  gsl_vector_complex *work = gsl_vector_complex_alloc (lwork);
  gsl_vector_complex *rwork = gsl_vector_complex_alloc (5 * L);

  int status;
  size_t i, j;

  for (i = 0; i < M; i++)
    for (j = 0; j < N; j++)
      {
	gsl_complex z = gsl_complex_rect (1.0 / (i + j + 1.0), 0.0);
	gsl_matrix_complex_set (m, i, j, z);
      }

  zgesvd_ ("A", "A",		/*  these are char * pointer arguments */
	   (int *) &(m->size1), (int *) &(m->size2),
	   m->data, (int *) &(m->tda),
	   s->data,
	   u->data, (int *) &(u->tda),
	   v->data, (int *) &(v->tda),
	   work->data, &lwork, rwork->data, &status);

  gsl_matrix_complex_fprintf (stdout, m, "%g");
  printf ("\n");
  gsl_matrix_complex_fprintf (stdout, u, "%g");
  printf ("\n");
  gsl_matrix_complex_fprintf (stdout, v, "%g");
  printf ("\n");

  gsl_vector_fprintf (stdout, s, "%g");
  printf ("\n");
}


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]