This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl interface to LAPACK?
- To: j dot l dot gomez-dans at sheffield dot ac dot uk
- Subject: Re: gsl interface to LAPACK?
- From: Brian Gough <bjg at network-theory dot co dot uk>
- Date: Fri, 20 Jul 2001 13:53:27 +0100 (BST)
- Cc: gsl-discuss at sourceware dot cygnus dot com
- References: <20010720124828.A2446@shef.ac.uk>
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");
}