%perlcode %{

use strict;
use Carp qw/croak/;
use Scalar::Util 'blessed';
use Math::Complex;
use Data::Dumper;

use Math::GSL           qw/:all/;
use Math::GSL::Errno    qw/:all/;
use Math::GSL::Eigen    qw/:all/;
use Math::GSL::Test     qw/is_similar/;
use Math::GSL::BLAS     qw/gsl_blas_zgemm/;
use Math::GSL::CBLAS    qw/$CblasNoTrans/;
use Math::GSL::Complex  qw/:all/;
use Math::GSL::Permutation;
use Math::GSL::VectorComplex qw/:all/;

use Math::GSL::Linalg qw/
    gsl_linalg_complex_LU_decomp
    gsl_linalg_complex_LU_det
    gsl_linalg_complex_LU_lndet
    gsl_linalg_complex_LU_invert
/;

use overload
    '*'      => \&_multiplication,
    '+'      => \&_addition,
    '-'      => \&_subtract,
    '=='     => \&_equal,
    '!='     => \&_not_equal,
#    'abs'    => \&_abs,
    fallback => 1;

our @EXPORT_OK = qw/
                gsl_matrix_complex_alloc
                gsl_matrix_complex_calloc
                gsl_matrix_complex_alloc_from_block
                gsl_matrix_complex_alloc_from_matrix
                gsl_vector_complex_alloc_row_from_matrix
                gsl_vector_complex_alloc_col_from_matrix
                gsl_matrix_complex_free
                gsl_matrix_complex_submatrix
                gsl_matrix_complex_row
                gsl_matrix_complex_column
                gsl_matrix_complex_diagonal
                gsl_matrix_complex_subdiagonal
                gsl_matrix_complex_superdiagonal
                gsl_matrix_complex_subrow
                gsl_matrix_complex_subcolumn
                gsl_matrix_complex_view_array
                gsl_matrix_complex_view_array_with_tda
                gsl_matrix_complex_view_vector
                gsl_matrix_complex_view_vector_with_tda
                gsl_matrix_complex_const_submatrix
                gsl_matrix_complex_const_row
                gsl_matrix_complex_const_column
                gsl_matrix_complex_const_diagonal
                gsl_matrix_complex_const_subdiagonal
                gsl_matrix_complex_const_superdiagonal
                gsl_matrix_complex_const_subrow
                gsl_matrix_complex_const_subcolumn
                gsl_matrix_complex_const_view_array
                gsl_matrix_complex_const_view_array_with_tda
                gsl_matrix_complex_const_view_vector
                gsl_matrix_complex_const_view_vector_with_tda
                gsl_matrix_complex_get
                gsl_matrix_complex_set
                gsl_matrix_complex_ptr
                gsl_matrix_complex_const_ptr
                gsl_matrix_complex_set_zero
                gsl_matrix_complex_set_identity
                gsl_matrix_complex_set_all
                gsl_matrix_complex_fread
                gsl_matrix_complex_fwrite
                gsl_matrix_complex_fscanf
                gsl_matrix_complex_fprintf
                gsl_matrix_complex_memcpy
                gsl_matrix_complex_swap
                gsl_matrix_complex_swap_rows
                gsl_matrix_complex_swap_columns
                gsl_matrix_complex_swap_rowcol
                gsl_matrix_complex_transpose
                gsl_matrix_complex_transpose_memcpy
                gsl_matrix_complex_isnull
                gsl_matrix_complex_ispos
                gsl_matrix_complex_isneg
                gsl_matrix_complex_add
                gsl_matrix_complex_sub
                gsl_matrix_complex_mul_elements
                gsl_matrix_complex_div_elements
                gsl_matrix_complex_scale
                gsl_matrix_complex_add_constant
                gsl_matrix_complex_add_diagonal
                gsl_matrix_complex_get_row
                gsl_matrix_complex_get_col
                gsl_matrix_complex_set_row
                gsl_matrix_complex_set_col
/;

sub _multiplication {
    my ($left,$right) = @_;
    my $lcopy = $left->copy;

    if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') ) {
        if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
            return _dot_product($left,$right);
            #gsl_matrix_complex_mul_elements($lcopy->raw, $right->raw);
        } else {
            croak "Math::GSL::MatrixComplex - multiplication of elements of matrices must be called with two objects matrices and must have the same number of columns and rows";
        }
    } else {
        gsl_matrix_complex_scale($lcopy->raw, $right);
    }
    return $lcopy;
}

sub _dot_product {
    my ($left,$right) = @_;

    croak "Math::GSL::Matrix - incompatible matrices in multiplication"
        unless ( blessed $right && $right->isa('Math::GSL::MatrixComplex') and $left->rows == $right->cols );
    my $C = Math::GSL::MatrixComplex->new($left->rows, $right->cols);
    my $complex = gsl_complex_rect(1,0);
    gsl_blas_zgemm($CblasNoTrans, $CblasNoTrans, $complex, $left->raw, $right->raw, $complex, $C->raw);
    return $C;
}

sub _addition {
    my ($left, $right) = @_;

    my $lcopy = $left->copy;

    if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') && blessed $left && $left->isa('Math::GSL::MatrixComplex') ) {
        if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
            gsl_matrix_complex_add($lcopy->raw, $right->raw);
        } else {
            croak "Math::GSL - addition of matrices must be called with two objects matrices and must have the same number of columns and rows";
        }
    } else {
        gsl_matrix_complex_add_constant($lcopy->raw, $right);
    }
    return $lcopy;
}

sub _subtract {
    my ($left, $right) = @_;

    my $lcopy = $left->copy;

    if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') && blessed $left && $left->isa('Math::GSL::MatrixComplex') ) {
        if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
            gsl_matrix_complex_sub($lcopy->raw, $right->raw);
        } else {
            croak "Math::GSL - subtraction of matrices must be called with two objects matrices and must have the same number of columns and rows";
        }
    } else {
        gsl_matrix_complex_add_constant($lcopy->raw, $right*-1);
    }
    return $lcopy;
}

sub _equal {
    my ($left, $right) = @_;
    return is_similar( [ map { Re $_ } $left->as_list ],
                       [ map { Re $_ } $right->as_list ]) &&
           is_similar( [ map { Im $_ } $left->as_list ],
                       [ map { Im $_ } $right->as_list ]);
}

sub _not_equal {
    my ($left, $right) = @_;
    return !_equal($left,$right);
}

sub copy {
    my $self = shift;
    my $copy = Math::GSL::MatrixComplex->new( $self->rows, $self->cols );
    if ( gsl_matrix_complex_memcpy($copy->raw, $self->raw) != $GSL_SUCCESS ) {
        croak "Math::GSL - error copying memory, aborting";
    }
    return $copy;
}
our %EXPORT_TAGS = ( all => \@EXPORT_OK  );

=encoding utf8

=head1 NAME

Math::GSL::MatrixComplex - Complex Matrices

=head1 SYNOPSIS

    use Math::GSL::MatrixComplex qw/:all/;
    my $matrix1 = Math::GSL::MatrixComplex->new(5,5);  # OO interface
    my $matrix3 = $matrix1 + $matrix1;
    my $matrix4 = $matrix1 - $matrix1;
    if($matrix1 == $matrix4) ...
    if($matrix1 != $matrix3) ...

    my $matrix2 = gsl_matrix_complex_alloc(5,5);        # standard interface

=head1 Objected Oriented Interface to GSL Math::GSL::MatrixComplex

=head2 new()

Creates a new MatrixComplex object of the given size.

    my $matrix = Math::GSL::MatrixComplex->new(10,10);
=cut

sub new
{
    my ($class, $rows, $cols) = @_;
    my $this = {};
    my $matrix;
    if ( defined $rows       && defined $cols &&
        $rows > 0            && $cols > 0     &&
        (int $rows == $rows) && (int $cols == $cols)){

        $matrix  = gsl_matrix_complex_alloc($rows,$cols);
    } else {
        croak( __PACKAGE__.'::new($x,$y) - $x and $y must be positive integers');
    }
    gsl_matrix_complex_set_zero($matrix);
    $this->{_matrix} = $matrix;
    ($this->{_rows}, $this->{_cols}) = ($rows,$cols);
    bless $this, $class;
}
=head2 raw()

Get the underlying GSL matrix object created by SWIG, useful for using gsl_matrix_* functions which do not have an OO counterpart.

    my $matrix     = Math::GSL::MatrixComplex->new(3,3);
    my $gsl_matrix = $matrix->raw;
    my $stuff      = gsl_matrix_complex_get($gsl_matrix, 1, 2);

=cut
sub raw  { (shift)->{_matrix} }
=head2 rows()

Returns the number of rows in the matrix.

    my $rows = $matrix->rows;
=cut

sub rows { (shift)->{_rows}   }

=head2 cols()

Returns the number of columns in the matrix.

    my $cols = $matrix->cols;
=cut

sub cols { (shift)->{_cols}   }

=head2 as_vector()

Returns a 1xN or Nx1 matrix as a Math::GSL::VectorComplex object. Dies if called on a matrix that is not a single row or column. Useful for turning the output of C<col()> or C<row()> into a vector, like so:

    my $vector1 = $matrix->col(0)->as_vector;
    my $vector2 = $matrix->row(1)->as_vector;

=cut

sub as_vector($)
{
    my ($self) = @_;
    croak(__PACKAGE__.'::as_vector() - must be a single row or column matrix') unless ($self->cols == 1 or $self->rows == 1);

    # TODO: there is a faster way to do this
    return Math::GSL::VectorComplex->new([ $self->as_list ] );

}

=head2  as_list()

Get the contents of a Math::GSL::Matrix object as a Perl list.

    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    ...
    my @matrix = $matrix->as_list;
=cut

sub as_list($)
{
    my $self = shift;
    my @matrix;
    for my $row ( 0 .. $self->rows-1) {
       push @matrix, map { gsl_matrix_complex_get($self->raw, $row, $_) } (0 .. $self->cols-1 );
    }
    return map {
            Math::Complex->make(
                gsl_real($_),
                gsl_imag($_))
            } @matrix;
}

=head2 row()

Returns a row matrix of the row you enter.

    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    ...
    my $matrix_row = $matrix->row(0);

=cut

sub row
{
    my ($self, $row) = @_;
    croak (__PACKAGE__.'::$matrix->row($row) - invalid $row value')
        unless (($row < $self->rows) and $row >= 0);

   my $rowmat = Math::GSL::MatrixComplex->new(1,$self->cols);

   for my $n (0 .. $self->cols-1) {
        gsl_matrix_complex_set( $rowmat->raw, 0, $n,
            gsl_matrix_complex_get($self->raw, $row, $n)
        );
   }

    return $rowmat;
}

=head2 col()

Returns a col matrix of the column you enter.

    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    ...
    my $matrix_col = $matrix->col(0);

=cut

sub col
{
    my ($self, $col) = @_;
    croak (__PACKAGE__."::\$matrix->col(\$col) - $col not a valid column")
        unless (defined $col && $col < $self->cols and $col >= 0);

    my $colmat = Math::GSL::MatrixComplex->new($self->rows, 1);

    map { gsl_matrix_complex_set($colmat->raw, $_, 0,
            gsl_matrix_complex_get($self->raw, $_, $col),
            )
        } (0 .. $self->rows-1);

    return $colmat;
}

=head2 set_row()

Sets a the values of a row with the elements of an array.

    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    $matrix->set_row(0, [8, 6, 2]);

You can also set multiple rows at once with chained calls:
    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    $matrix->set_row(0, [8, 6, 2])
           ->set_row(1, [2, 4, 1]);
    ...

=cut

sub set_row {
    my ($self, $row, $values) = @_;
    my $length = $#$values;
    die __PACKAGE__.'::set_row($x, $values) - $values must be a nonempty array reference' if $length == -1;
    die __PACKAGE__.'::set_row($x, $values) - $x must be a valid row number' if ($row < 0 || $row >= $self->rows);
    die __PACKAGE__.'::set_row($x, $values) - $values must contains the same number of elements as there is columns in the matrix' if($length != $self->cols-1);
    # $values may have Math::Complex objects
    @$values = map { Math::GSL::Complex::gsl_complex_rect(Re($_), Im($_)) } @$values;
    # warn Dumper [ @$values ];
    map { gsl_matrix_complex_set($self->raw, $row, $_, $values->[$_]) } (0..$length);
    return $self;
}

=head2 set_col()

Sets a the values of a column with the elements of an array.

    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    $matrix->set_col(0, [8, 6, 2]);

You can also set multiple columns at once with chained calls:
    my $matrix = Math::GSL::MatrixComplex->new(3,3);
    $matrix->set_col(0, [8, 6, 2])
           ->set_col(1, [2, 4, 1]);
    ...

=cut

sub set_col {
    my ($self, $col, $values) = @_;
    my $length = $#$values;
    die __PACKAGE__.'::set_col($x, $values) - $values must be a nonempty array reference' if $length == -1;
    die __PACKAGE__.'::set_col($x, $values) - $x must be a valid column number' if ($col < 0 || $col >= $self->cols);
    die __PACKAGE__.'::set_col($x, $values) - $values must contains the same number of elements as there is rowss in the matrix' if($length != $self->rows-1);
    # $values may have Math::Complex objects
    @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values;
    map { gsl_matrix_complex_set($self->raw, $_, $col, $values->[$_]) } (0..$length);
    return $self;
}

=head2 is_square()

Returns true if a matrix is square, i.e. it has the same number of rows as columns, false otherwise.

=cut

sub is_square($)
{
    my $self = shift;
    return ($self->rows == $self->cols) ? 1 : 0 ;
}

=head2 det()

Returns the determinant of a matrix (computed by LU decomposition) or dies if called on a non-square matrix.

    my $det = $matrix->det();

=cut

sub det($)
{
    my $self = shift;
    croak(__PACKAGE__."- determinant only exists for square matrices") unless $self->is_square;

    my $p  = Math::GSL::Permutation->new( $self->rows );
    my $LU = $self->copy;

    my $s = gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);

    # $z is a gsl_complex
    my $z = gsl_linalg_complex_LU_det($LU->raw, $s );
    return Math::Complex->make( gsl_real($z), gsl_imag($z) );
}

=head2 zero()

Set a matrix to the zero matrix.

    $matrix->zero;

=cut

sub zero # brrr!
{
    my $self=shift;
    gsl_matrix_complex_set_zero($self->raw);
    return $self;
}

=head2 identity()

Set a matrix to the identity matrix, i.e. one on the diagonal and zero elsewhere.

    my $I = $matrix->identity;

=cut

sub identity
{
    my $self=shift;
    gsl_matrix_complex_set_identity($self->raw);
    return $self;
}

=head2 inverse()

Returns the inverse of a matrix or dies when called on a non-square matrix.

    my $inverse = $matrix->inverse;

=cut

sub inverse($)
{
    my $self = shift;
    croak(__PACKAGE__."- inverse only exists for square matrices") unless $self->is_square;

    my $p  = Math::GSL::Permutation->new( $self->rows );
    my $LU = $self->copy;
    my $inverse = $self->copy;

    # should check return status
    gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);
    gsl_linalg_complex_LU_invert($LU->raw, $p->raw,$inverse->raw);

    return $inverse;

}

=head2 is_hermitian()

Returns true if the matrix is hermitian, false otherwise

    my $test = $matrix->is_hermitian;

=cut

sub is_hermitian()
{
    my ($self) = @_;
    my $test = $self->is_square;
    my $transpose = $self->copy;
    gsl_matrix_complex_transpose($transpose->raw);

    if($test == 1) {
        for my $row (0..$self->rows - 1) {
            map { gsl_matrix_complex_set($transpose->raw, $row, $_, gsl_complex_conjugate(gsl_matrix_complex_get($transpose->raw, $row, $_))) } (0..$self->cols - 1);
        }
        if($self != $transpose){
            $test = 0;
        }
    }
    return $test;
}
=head2 eigenvalues()

=cut


sub eigenvalues($)
{
    my $self=shift;
    my ($r,$c) = ($self->rows,$self->cols);
    croak "Math::GSL::MatrixComplex : \$matrix->eigenvalues - \$matrix must be square" unless ($r == $c);
    my $vector = Math::GSL::Vector->new($r);
    my $eigen  = gsl_eigen_herm_alloc($r);


    # GSL has no generalized complex matrix routines
    # can only continue if the matrix is hermitian
    croak (__PACKAGE__."::eigenvalues : non-hermitian matrices are not currently supported")  unless $self->is_hermitian;
    gsl_eigen_herm($self->raw, $vector->raw, $eigen);

    return $vector->as_list;
}

=head2 lndet()

Returns the natural log of the absolute value of the determinant of a matrix (computed by LU decomposition) or dies if called on a non-square matrix.

    my $lndet = $matrix->lndet();

=cut

sub lndet($)
{
    my $self = shift;
    croak(__PACKAGE__."- log determinant only exists for square matrices") unless $self->is_square;

    my $p  = Math::GSL::Permutation->new( $self->rows );
    my $LU = $self->copy;

    gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);
    return gsl_linalg_complex_LU_lndet($LU->raw);
}
=head1 DESCRIPTION

=over 1

=item C<gsl_matrix_complex_alloc >

=item C<gsl_matrix_complex_calloc >

=item C<gsl_matrix_complex_alloc_from_block >

=item C<gsl_matrix_complex_alloc_from_matrix >

=item C<gsl_vector_complex_alloc_row_from_matrix >

=item C<gsl_vector_complex_alloc_col_from_matrix >

=item C<gsl_matrix_complex_free >

=item C<gsl_matrix_complex_submatrix >

=item C<gsl_matrix_complex_row >

=item C<gsl_matrix_complex_column >

=item C<gsl_matrix_complex_diagonal >

=item C<gsl_matrix_complex_subdiagonal >

=item C<gsl_matrix_complex_superdiagonal >

=item C<gsl_matrix_complex_subrow >

=item C<gsl_matrix_complex_subcolumn >

=item C<gsl_matrix_complex_view_array >

=item C<gsl_matrix_complex_view_array_with_tda >

=item C<gsl_matrix_complex_view_vector >

=item C<gsl_matrix_complex_view_vector_with_tda >

=item C<gsl_matrix_complex_const_submatrix >

=item C<gsl_matrix_complex_const_row >

=item C<gsl_matrix_complex_const_column >

=item C<gsl_matrix_complex_const_diagonal >

=item C<gsl_matrix_complex_const_subdiagonal >

=item C<gsl_matrix_complex_const_superdiagonal >

=item C<gsl_matrix_complex_const_subrow >

=item C<gsl_matrix_complex_const_subcolumn >

=item C<gsl_matrix_complex_const_view_array >

=item C<gsl_matrix_complex_const_view_array_with_tda >

=item C<gsl_matrix_complex_const_view_vector >

=item C<gsl_matrix_complex_const_view_vector_with_tda >

=item C<gsl_matrix_complex_get>

=item C<gsl_matrix_complex_set>

=item C<gsl_matrix_complex_ptr>

=item C<gsl_matrix_complex_const_ptr>

=item C<gsl_matrix_complex_set_zero >

=item C<gsl_matrix_complex_set_identity >

=item C<gsl_matrix_complex_set_all >

=item C<gsl_matrix_complex_fread >

=item C<gsl_matrix_complex_fwrite >

=item C<gsl_matrix_complex_fscanf >

=item C<gsl_matrix_complex_fprintf >

=item C<gsl_matrix_complex_memcpy>

=item C<gsl_matrix_complex_swap>

=item C<gsl_matrix_complex_swap_rows>

=item C<gsl_matrix_complex_swap_columns>

=item C<gsl_matrix_complex_swap_rowcol>

=item C<gsl_matrix_complex_transpose >

=item C<gsl_matrix_complex_transpose_memcpy >

=item C<gsl_matrix_complex_isnull >

=item C<gsl_matrix_complex_ispos >

=item C<gsl_matrix_complex_isneg >

=item C<gsl_matrix_complex_add >

=item C<gsl_matrix_complex_sub >

=item C<gsl_matrix_complex_mul_elements >

=item C<gsl_matrix_complex_div_elements >

=item C<gsl_matrix_complex_scale >

=item C<gsl_matrix_complex_add_constant >

=item C<gsl_matrix_complex_add_diagonal >

=item C<gsl_matrix_complex_get_row>

=item C<gsl_matrix_complex_get_col>

=item C<gsl_matrix_complex_set_row>

=item C<gsl_matrix_complex_set_col>

=back

For more information on the functions, we refer you to the GSL official documentation
L<http://www.gnu.org/software/gsl/manual/html_node/>


=head1 AUTHORS

Jonathan "Duke" Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2008-2024 Jonathan "Duke" Leto and Thierry Moisan

This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=cut

%}