new release candidate 1.2.91

new release candidate 1.2.91

[Brian Gough]

I've upgraded libtool, fixed configure for detection of log1p and
improved the test messages.

http://www.getafile.com/cgi-bin/merlot/get/network-theory/gsl/gsl-1.2.91.tar.gz

Please retest and report the results to
gsl-discuss@sources.redhat.com.  Thanks.

Brian

Re: new release candidate 1.2.91

[Gert Van den Eynde]

Hi,

tested on 

Linux pc1727 2.4.19-4GB #1 Fri Sep 13 13:14:56 UTC 2002 i686 unknown

on SuSE 8.1 using gcc version 3.2.1


All tests passed.


Bye,

Gert

Re: new release candidate 1.2.91

[Philip Kendall]

On Sun, Dec 15, 2002 at 08:24:42PM +0000, Brian Gough wrote:
> I've upgraded libtool, fixed configure for detection of log1p and
> improved the test messages.
> 
> http://www.getafile.com/cgi-bin/merlot/get/network-theory/gsl/gsl-1.2.91.tar.gz

Same problem as in my previous message as with 1.2.90 on Solaris 8/gcc
3.1.

On a more positive note, using Sun's cc:

cass30:pak:~/data2/gsl-1.2.91$ cc -V
cc: Sun WorkShop 6 update 2 C 5.3 Patch 111679-08 2002/05/09

everything works (bar lots of warnings about 'end-of-loop code not
reached' while compiling) without optimization (the default).
However, if I explicitly add CFLAGS=-fast, `make check' fails early:

FAIL: gsl_isinf(inf) (0 observed vs 1 expected)
FAIL: gsl_isinf(-inf) (0 observed vs -1 expected)
FAIL: gsl_finite(inf) (1 observed vs 0 expected)
FAIL: gsl_finite(nan) (1 observed vs 0 expected)

Phil

-- 
"Buffy want beer." "You can't have beer."
"Want. Beer." "Giles, don't make caveslayer unhappy."
                                        Buffy, Giles and Xander: Beer Bad

Re: new release candidate 1.2.91

[Slaven Peles]

On Sunday 15 December 2002 03:24 pm, Brian Gough wrote:
> I've upgraded libtool, fixed configure for detection of log1p and
> improved the test messages.
>
> http://www.getafile.com/cgi-bin/merlot/get/network-theory/gsl/gsl-1.2.91.ta
>r.gz
>
> Please retest and report the results to
> gsl-discuss@sources.redhat.com.  Thanks.
>
> Brian

Compiles, passes all tests on Compaq Alpha

Compaq C V6.3-025 on Compaq Tru64 UNIX V5.1 (Rev. 732)
Compiler Driver V6.3-026 (sys) cc Driver

Slaven

Re: new release candidate 1.2.91

[Grant Ingram]

Brian Gough wrote:

> Please retest and report the results to
> gsl-discuss@sources.redhat.com.  Thanks.

o.k so this is actually testing rather than re-testing but...

Compiling on a cygwin installation gives the following failure:-

Output from "make check":-

make[2]: Nothing to be done for `test_gsl_histogram.sh'.
make[2]: Leaving directory `/home/grant/gsl-1.2.91'
make  check-TESTS
make[2]: Entering directory `/home/grant/gsl-1.2.91'
test.exp.1.tmp test.obs.1.tmp differ: char 6, line 1
FAIL: test_gsl_histogram.sh
===================
1 of 1 tests failed
===================
make[2]: Leaving directory `/home/grant/gsl-1.2.91'
make[1]: Leaving directory `/home/grant/gsl-1.2.91'

I'd previously installed (and actually used) the gsl version 1.2 on the same system
without problems.  Though I don't recall if I ran make check.  The only hint I
could find in the INSTALL file was to try installing without the shared libraries
"./configure --disable-shared" which I tried and "make check" produced exactly the
same error message.  The INSTALL file also mentions that the test should fail with
gcc v2.95/2.96 in the eigen directory - but that doesn't seem to have happened -
perhaps I have a recent enough version of gcc?

Anyway I hope this is of some use.

Cheers,

Grant.

******************
System Details:-
$ uname --all
CYGWIN_98-4.10 ENG9023 1.3.12(0.54/3/2) 2002-07-06 02:16 i686 unknown

$ gcc --version
2.95.3-5

Re: new release candidate 1.2.91

[Brian Gough]

Grant Ingram writes:
 > make[2]: Entering directory `/home/grant/gsl-1.2.91'
 > test.exp.1.tmp test.obs.1.tmp differ: char 6, line 1
 > FAIL: test_gsl_histogram.sh

Thanks.  If you look in the shell script test_gsl_histogram.sh you
could try the commands manually and see what the problem is.

Brian

Re: new release candidate 1.2.91

[Lukas Dobrek]

On Sun, Dec 15, 2002 at 08:24:42PM +0000, Brian Gough wrote:
> I've upgraded libtool, fixed configure for detection of log1p and
> improved the test messages.
> 
> http://www.getafile.com/cgi-bin/merlot/get/network-theory/gsl/gsl-1.2.91.tar.gz
> 
> Please retest and report the results to
> gsl-discuss@sources.redhat.com.  Thanks.
On Axp Linux the situation is like before. 
gcc-3.2 all test passed. 
gcc-2.95.4 everythings fine appart of eigen/herm but this is compiler bug.

Take Care

Lukasz

PS. Brian sorry for sending to you directly it was mutt misconfiguration, 
    want happend again. 




-- 
£ukasz Dobrek
   An optimist believes that we live in the best of all possible worlds.
   A pessimist is sure that this must be so.
http://www.pld-linux.org

Re: new release candidate 1.2.91

[Brian Gough]

Thanks to everyone for the testing. I'll go with 1.2.91.

Brian

Re: new release candidate 1.2.91

[John Marraffino]

| Grant Ingram writes:
|  > make[2]: Entering directory `/home/grant/gsl-1.2.91'
|  > test.exp.1.tmp test.obs.1.tmp differ: char 6, line 1
|  > FAIL: test_gsl_histogram.sh
| 
| Thanks.  If you look in the shell script test_gsl_histogram.sh you
| could try the commands manually and see what the problem is.
| 
| Brian

I noticed the same thing and did the same experiment you suggest. The
problem is that, on CYGWIN systems, the file text.exp.1.tmp, produced
via the cat command, somehow acquires carriage returns as well as
newlines. The file test.obs.1.tmp has only newlines and, consequently,
cmp declares them different.

Cheers
John
-
=======================================================================
John Marraffino                                       marafino@fnal.gov
Fermi National Accelerator Lab                     PHONE: (630)840-4483
Computing Division/Computational Physics Department  FAX: (630)840-2783
=======================================================================
"The surest sign that intelligent life exists elsewhere in the universe
is that it has never tried to contact us."
				     Calvin and Hobbes (Bill Watterson)

Re: new release candidate 1.2.91

[Brian Gough]

John Marraffino writes:
 > I noticed the same thing and did the same experiment you suggest. The
 > problem is that, on CYGWIN systems, the file text.exp.1.tmp, produced
 > via the cat command, somehow acquires carriage returns as well as
 > newlines. The file test.obs.1.tmp has only newlines and, consequently,
 > cmp declares them different.

If any CYGWIN users would like to suggest a patch to the test I will
add it (maybe to use 'echo' instead of 'cat' ??).

LU "bug" and fix

[Slaven Peles]

[-- Attachment #1: Type: text/plain, Size: 368 bytes --]

LU decomposition does not give an error message in case that input matrix is 
singular, and I guess it should... I changed functions gsl_linalg_LU_decomp 
and gsl_linalg_complex_LU_decomp (in files linalg/lu.c and linalg/luc.c) so 
they cry foul whenever user puts a singular matrix as the input. Perhaps this 
should be added to gsl cvs at some point?

Cheers,
Slaven

[-- Attachment #2: luc.c --]
[-- Type: text/x-c, Size: 8454 bytes --]

/* linalg/luc.c
 * 
 * Copyright (C) 2001 Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex_math.h>

#include <gsl/gsl_linalg.h>

/* Factorise a general N x N complex matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.  
 * 
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation. 
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (p->size != A->size1)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i, j, k;

      *signum = 1;
      gsl_permutation_init (p);

      for (j = 0; j < N - 1; j++)
	{
	  /* Find maximum in the j-th column */

	  gsl_complex ajj = gsl_matrix_complex_get (A, j, j);
          double max = gsl_complex_abs (ajj);
	  size_t i_pivot = j;

	  for (i = j + 1; i < N; i++)
	    {
	      gsl_complex aij = gsl_matrix_complex_get (A, i, j);
              double ai = gsl_complex_abs (aij);

	      if (ai > max)
		{
		  max = ai;
		  i_pivot = i;
		}
	    }

	  if (i_pivot != j)
	    {
	      gsl_matrix_complex_swap_rows (A, j, i_pivot);
	      gsl_permutation_swap (p, j, i_pivot);
	      *signum = -(*signum);
	    }

	  ajj = gsl_matrix_complex_get (A, j, j);

	  if (!(GSL_REAL(ajj) == 0.0 && GSL_IMAG(ajj) == 0.0))
	    {
	      for (i = j + 1; i < N; i++)
		{
		  gsl_complex aij_orig = gsl_matrix_complex_get (A, i, j);
                  gsl_complex aij = gsl_complex_div (aij_orig, ajj);
		  gsl_matrix_complex_set (A, i, j, aij);

		  for (k = j + 1; k < N; k++)
		    {
		      gsl_complex aik = gsl_matrix_complex_get (A, i, k);
		      gsl_complex ajk = gsl_matrix_complex_get (A, j, k);
                      
                      /* aik = aik - aij * ajk */

                      gsl_complex aijajk = gsl_complex_mul (aij, ajk);
                      gsl_complex aik_new = gsl_complex_sub (aik, aijajk);

		      gsl_matrix_complex_set (A, i, k, aik_new);
		    }
		}
	    }
	  else /* If LU matrix is singular exit with error */
	    {
	      GSL_ERROR ("LU matrix must not be singular", GSL_EINVAL);
	    }
	}
      
      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_complex_memcpy (x, b);

      /* Solve for x */

      gsl_linalg_complex_LU_svx (LU, p, x);

      return GSL_SUCCESS;
    }
}


int
gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* Apply permutation to RHS */

      gsl_permute_vector_complex (p, x);

      /* Solve for c using forward-substitution, L c = P b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}


int
gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Compute residual, residual = (A * x  - b) */

      gsl_vector_complex_memcpy (residual, b);

      {
        gsl_complex one = GSL_COMPLEX_ONE;
        gsl_complex negone = GSL_COMPLEX_NEGONE;
        gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, residual);
      }

      /* Find correction, delta = - (A^-1) * residual, and apply it */

      gsl_linalg_complex_LU_svx (LU, p, residual);

      {
        gsl_complex negone= GSL_COMPLEX_NEGONE;
        gsl_blas_zaxpy (negone, residual, x);
      }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
{
  size_t i, n = LU->size1;

  int status = GSL_SUCCESS;

  gsl_matrix_complex_set_identity (inverse);

  for (i = 0; i < n; i++)
    {
      gsl_vector_complex_view c = gsl_matrix_complex_column (inverse, i);
      int status_i = gsl_linalg_complex_LU_svx (LU, p, &(c.vector));

      if (status_i)
	status = status_i;
    }

  return status;
}

gsl_complex
gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex det = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
      det = gsl_complex_mul (det, zi);
    }

  return det;
}


double
gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);
      lndet += log (gsl_complex_abs (z));
    }

  return lndet;
}


gsl_complex
gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex phase = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);
      
      double r = gsl_complex_abs(z);

      if (r == 0)
	{
	  phase = gsl_complex_rect(0.0, 0.0);
	  break;
	}
      else
        {
          z = gsl_complex_div_real(z, r);
          phase = gsl_complex_mul(phase, z);
        }
    }

  return phase;
}

[-- Attachment #3: lu.c --]
[-- Type: text/x-c, Size: 7137 bytes --]

/* linalg/lu.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

#define REAL double

/* Factorise a general N x N matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.  
 * 
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation. 
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (p->size != A->size1)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i, j, k;

      *signum = 1;
      gsl_permutation_init (p);

      for (j = 0; j < N - 1; j++)
	{
	  /* Find maximum in the j-th column */

	  REAL ajj, max = fabs (gsl_matrix_get (A, j, j));
	  size_t i_pivot = j;

	  for (i = j + 1; i < N; i++)
	    {
	      REAL aij = fabs (gsl_matrix_get (A, i, j));

	      if (aij > max)
		{
		  max = aij;
		  i_pivot = i;
		}
	    }

	  if (i_pivot != j)
	    {
	      gsl_matrix_swap_rows (A, j, i_pivot);
	      gsl_permutation_swap (p, j, i_pivot);
	      *signum = -(*signum);
	    }

	  ajj = gsl_matrix_get (A, j, j);

	  if (ajj != 0.0)
	    {
	      for (i = j + 1; i < N; i++)
		{
		  REAL aij = gsl_matrix_get (A, i, j) / ajj;
		  gsl_matrix_set (A, i, j, aij);

		  for (k = j + 1; k < N; k++)
		    {
		      REAL aik = gsl_matrix_get (A, i, k);
		      REAL ajk = gsl_matrix_get (A, j, k);
		      gsl_matrix_set (A, i, k, aik - aij * ajk);
		    }
		}
	    }
	  else /* If LU matrix is singular exit with error */
	    {
	      GSL_ERROR ("LU matrix must not be singular", GSL_EINVAL);
	    }
	}
      
      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve for x */

      gsl_linalg_LU_svx (LU, p, x);

      return GSL_SUCCESS;
    }
}


int
gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* Apply permutation to RHS */

      gsl_permute_vector (p, x);

      /* Solve for c using forward-substitution, L c = P b */

      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}


int
gsl_linalg_LU_refine (const gsl_matrix * A, const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Compute residual, residual = (A * x  - b) */

      gsl_vector_memcpy (residual, b);
      gsl_blas_dgemv (CblasNoTrans, 1.0, A, x, -1.0, residual);

      /* Find correction, delta = - (A^-1) * residual, and apply it */

      gsl_linalg_LU_svx (LU, p, residual);
      gsl_blas_daxpy (-1.0, residual, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)
{
  size_t i, n = LU->size1;

  int status = GSL_SUCCESS;

  gsl_matrix_set_identity (inverse);

  for (i = 0; i < n; i++)
    {
      gsl_vector_view c = gsl_matrix_column (inverse, i);
      int status_i = gsl_linalg_LU_svx (LU, p, &(c.vector));

      if (status_i)
	status = status_i;
    }

  return status;
}

double
gsl_linalg_LU_det (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  double det = (double) signum;

  for (i = 0; i < n; i++)
    {
      det *= gsl_matrix_get (LU, i, i);
    }

  return det;
}


double
gsl_linalg_LU_lndet (gsl_matrix * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      lndet += log (fabs (gsl_matrix_get (LU, i, i)));
    }

  return lndet;
}


int
gsl_linalg_LU_sgndet (gsl_matrix * LU, int signum)
{
  size_t i, n = LU->size1;

  int s = signum;

  for (i = 0; i < n; i++)
    {
      double u = gsl_matrix_get (LU, i, i);

      if (u < 0)
	{
	  s *= -1;
	}
      else if (u == 0)
	{
	  s = 0;
	  break;
	}
    }

  return s;
}

Re: LU "bug" and fix

[Brian Gough]

Slaven Peles writes:
 > LU decomposition does not give an error message in case that input
 > matrix is singular

A singular matrix should still have a valid LU decomposition
(L*U=original matrix) so it's not an error in that case.

Brian