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