Re: gsl interface to LAPACK?

Re: gsl interface to LAPACK?

[Brian Gough]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 2585 bytes --]

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");
}

gsl interface to LAPACK?

[José Luis Gómez Dans]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 1116 bytes --]

Hi,
	I am beginning to use GSL, and I am both a newbie in both GSL
and C. So I apologise if I ask something really stupid.

	I was playing with GSL to convert some code I have in Fortran
95. Mostly, I got it to work. The problem is that I need to solve a
generalised eiganvalue/vector problem (complex, non-Hermitian).
According to the manuals, you need to use LAPACK for this. I understand
that (I have used LAPACK in F77), but I now have all my matrices and
vectors as gsl_complex_matrix. If I want to use LAPACK, I have to get
the data into a CLAPACK friendly format, which is rather tedious and
error-prone (there might be many an instance when this might apply).
Moreover, the matrices might be Hermitian or not...

	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.

	Thank you for your time,
	José
-- 
José L Gómez Dans			PhD student
Tel: +44 114 222 5582			Radar & Communications Group
FAX; +44 870 132 2990			Department of Electronic Engineering
					University of Sheffield UK