/ 
PDL-2.058
River stage three • 100 direct dependents • 108 total dependents
/ Libtmp/GSL/LINALG/gsl_linalg.pd

NAME

PDL::GSL::LINALG - PDL interface to linear algebra routines in GSL

SYNOPSIS

use PDL::LiteF;
use PDL::MatrixOps; # for 'x'
use PDL::GSL::LINALG;
my $A = pdl [
  [0.18, 0.60, 0.57, 0.96],
  [0.41, 0.24, 0.99, 0.58],
  [0.14, 0.30, 0.97, 0.66],
  [0.51, 0.13, 0.19, 0.85],
];
my $B = sequence(2,4); # column vectors
LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null);
# transpose so first dim means is vector, higher dims thread
LU_solve($lu, $p, $B->transpose, my $x=null);
$x = $x->inplace->transpose; # now can be matrix-multiplied

DESCRIPTION

This is an interface to the linear algebra package present in the GNU Scientific Library. Functions are named as in GSL, but with the initial gsl_linalg_ removed. They are provided in both real and complex double precision.

Currently only LU decomposition interfaces here. Pull requests welcome!

EOD

pp_addpm({At=>Bot},<<'EOD'); =head1 SEE ALSO

PDL

The GSL documentation for linear algebra is online at https://www.gnu.org/software/gsl/doc/html/linalg.html

LU decomposition of the given (real or complex) matrix.

EOF );

pp_def('LU_solve', HandleBad => 0, Pars => '[phys]LU(n,m); indx [phys]ipiv(p); [phys]B(n); [o,phys]x(n)', GenericTypes => [qw(C D)], Code => <<'EOF', gsl_matrix$TDC(,_complex) m; gsl_permutation p; gsl_vector$TDC(,_complex) b, x; int s; MATRIX_SETUP(m, $SIZE(m), $SIZE(n), $SIZE(n), $P(LU)) PERM_SETUP(p, $SIZE(p), $P(ipiv)) VECTOR_SETUP(b, $SIZE(n), $P(B)) VECTOR_SETUP(x, $SIZE(n), $P(x)) gsl_linalg$TDC(,_complex)_LU_solve(&m, &p, &b, &x); PERM_RETURN(p, $SIZE(p), $P(ipiv)) EOF Doc => <<'EOF', =for ref

Solve A x = B using the LU and permutation from "LU_decomp", real or complex.

EOF );

pp_def('LU_det', HandleBad => 0, Pars => '[phys]LU(n,m); int [phys]signum(); [o]det()', GenericTypes => [qw(C D)], Code => <<'EOF', gsl_matrix$TDC(,_complex) m; MATRIX_SETUP(m, $SIZE(m), $SIZE(n), $SIZE(n), $P(LU)) types (D) %{ $det() = gsl_linalg_LU_det(&m, $signum()); %} types (C) %{ gsl_complex z = gsl_linalg_complex_LU_det(&m, $signum()); $det() = GSL_REAL(z) + I*GSL_IMAG(z); %} EOF Doc => <<'EOF', =for ref

Find the determinant from the LU decomp.

EOF );

pp_def('solve_tridiag', HandleBad => 0, Pars => '[phys]diag(n); [phys]superdiag(n); [phys]subdiag(n); [phys]B(n); [o,phys]x(n)', GenericTypes => [qw(D)], Code => <<'EOF', gsl_vector d, sup, sub, b, x; VECTOR_SETUP(d, $SIZE(n), $P(diag)) VECTOR_SETUP(sup, $SIZE(n)-1, $P(superdiag)) VECTOR_SETUP(sub, $SIZE(n)-1, $P(subdiag)) VECTOR_SETUP(b, $SIZE(n), $P(B)) VECTOR_SETUP(x, $SIZE(n), $P(x)) #define CONST_VEC (const gsl_vector *) gsl_linalg_solve_tridiag( CONST_VEC &d, CONST_VEC &sup, CONST_VEC &sub, CONST_VEC &b, &x ); #undef CONST_VEC EOF Doc => <<'EOF', =for ref

Solve A x = B where A is a tridiagonal system. Real only, because GSL does not have a complex function.

EOF );

pp_done();