Update Copyright years
[Math-GSL.git] / pod / MatrixComplex.pod
blobcaefad9d17968aeaa5a20a20cb33ead1804a8a7c
1 %perlcode %{
3 use strict;
4 use Carp qw/croak/;
5 use Scalar::Util 'blessed';
6 use Math::Complex;
7 use Data::Dumper;
9 use Math::GSL           qw/:all/;
10 use Math::GSL::Errno    qw/:all/;
11 use Math::GSL::Eigen    qw/:all/;
12 use Math::GSL::Test     qw/is_similar/;
13 use Math::GSL::BLAS     qw/gsl_blas_zgemm/;
14 use Math::GSL::CBLAS    qw/$CblasNoTrans/;
15 use Math::GSL::Complex  qw/:all/;
16 use Math::GSL::Permutation;
17 use Math::GSL::VectorComplex qw/:all/;
19 use Math::GSL::Linalg qw/
20     gsl_linalg_complex_LU_decomp
21     gsl_linalg_complex_LU_det
22     gsl_linalg_complex_LU_lndet
23     gsl_linalg_complex_LU_invert
26 use overload
27     '*'      => \&_multiplication,
28     '+'      => \&_addition,
29     '-'      => \&_subtract,
30     '=='     => \&_equal,
31     '!='     => \&_not_equal,
32 #    'abs'    => \&_abs,
33     fallback => 1;
35 our @EXPORT_OK = qw/
36                 gsl_matrix_complex_alloc
37                 gsl_matrix_complex_calloc
38                 gsl_matrix_complex_alloc_from_block
39                 gsl_matrix_complex_alloc_from_matrix
40                 gsl_vector_complex_alloc_row_from_matrix
41                 gsl_vector_complex_alloc_col_from_matrix
42                 gsl_matrix_complex_free
43                 gsl_matrix_complex_submatrix
44                 gsl_matrix_complex_row
45                 gsl_matrix_complex_column
46                 gsl_matrix_complex_diagonal
47                 gsl_matrix_complex_subdiagonal
48                 gsl_matrix_complex_superdiagonal
49                 gsl_matrix_complex_subrow
50                 gsl_matrix_complex_subcolumn
51                 gsl_matrix_complex_view_array
52                 gsl_matrix_complex_view_array_with_tda
53                 gsl_matrix_complex_view_vector
54                 gsl_matrix_complex_view_vector_with_tda
55                 gsl_matrix_complex_const_submatrix
56                 gsl_matrix_complex_const_row
57                 gsl_matrix_complex_const_column
58                 gsl_matrix_complex_const_diagonal
59                 gsl_matrix_complex_const_subdiagonal
60                 gsl_matrix_complex_const_superdiagonal
61                 gsl_matrix_complex_const_subrow
62                 gsl_matrix_complex_const_subcolumn 
63                 gsl_matrix_complex_const_view_array
64                 gsl_matrix_complex_const_view_array_with_tda
65                 gsl_matrix_complex_const_view_vector
66                 gsl_matrix_complex_const_view_vector_with_tda
67                 gsl_matrix_complex_get
68                 gsl_matrix_complex_set
69                 gsl_matrix_complex_ptr
70                 gsl_matrix_complex_const_ptr
71                 gsl_matrix_complex_set_zero
72                 gsl_matrix_complex_set_identity
73                 gsl_matrix_complex_set_all
74                 gsl_matrix_complex_fread
75                 gsl_matrix_complex_fwrite
76                 gsl_matrix_complex_fscanf
77                 gsl_matrix_complex_fprintf
78                 gsl_matrix_complex_memcpy
79                 gsl_matrix_complex_swap
80                 gsl_matrix_complex_swap_rows
81                 gsl_matrix_complex_swap_columns
82                 gsl_matrix_complex_swap_rowcol
83                 gsl_matrix_complex_transpose
84                 gsl_matrix_complex_transpose_memcpy
85                 gsl_matrix_complex_isnull
86                 gsl_matrix_complex_ispos
87                 gsl_matrix_complex_isneg
88                 gsl_matrix_complex_add
89                 gsl_matrix_complex_sub
90                 gsl_matrix_complex_mul_elements
91                 gsl_matrix_complex_div_elements
92                 gsl_matrix_complex_scale
93                 gsl_matrix_complex_add_constant
94                 gsl_matrix_complex_add_diagonal
95                 gsl_matrix_complex_get_row
96                 gsl_matrix_complex_get_col
97                 gsl_matrix_complex_set_row
98                 gsl_matrix_complex_set_col
101 sub _multiplication {
102     my ($left,$right) = @_;
103     my $lcopy = $left->copy;
105     if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') ) {
106         if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
107             return _dot_product($left,$right);
108             #gsl_matrix_complex_mul_elements($lcopy->raw, $right->raw);
109         } else {
110             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";
111         } 
112     } else {
113         gsl_matrix_complex_scale($lcopy->raw, $right);
114     }
115     return $lcopy;
118 sub _dot_product {
119     my ($left,$right) = @_;
121     croak "Math::GSL::Matrix - incompatible matrices in multiplication"
122         unless ( blessed $right && $right->isa('Math::GSL::MatrixComplex') and $left->rows == $right->cols );
123     my $C = Math::GSL::MatrixComplex->new($left->rows, $right->cols);
124     my $complex = gsl_complex_rect(1,0);
125     gsl_blas_zgemm($CblasNoTrans, $CblasNoTrans, $complex, $left->raw, $right->raw, $complex, $C->raw);
126     return $C;
129 sub _addition {
130     my ($left, $right) = @_;
132     my $lcopy = $left->copy;
134     if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') && blessed $left && $left->isa('Math::GSL::MatrixComplex') ) {
135         if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
136             gsl_matrix_complex_add($lcopy->raw, $right->raw);
137         } else {
138             croak "Math::GSL - addition of matrices must be called with two objects matrices and must have the same number of columns and rows";
139         }
140     } else {
141         gsl_matrix_complex_add_constant($lcopy->raw, $right);
142     }
143     return $lcopy;
146 sub _subtract {
147     my ($left, $right) = @_;
149     my $lcopy = $left->copy;
151     if ( blessed $right && $right->isa('Math::GSL::MatrixComplex') && blessed $left && $left->isa('Math::GSL::MatrixComplex') ) {
152         if ( $left->cols == $right->cols && $left->rows == $right->rows ) {
153             gsl_matrix_complex_sub($lcopy->raw, $right->raw);
154         } else {
155             croak "Math::GSL - subtraction of matrices must be called with two objects matrices and must have the same number of columns and rows";
156         }
157     } else {
158         gsl_matrix_complex_add_constant($lcopy->raw, $right*-1);
159     }
160     return $lcopy;
163 sub _equal {
164     my ($left, $right) = @_;
165     return is_similar( [ map { Re $_ } $left->as_list ],
166                        [ map { Re $_ } $right->as_list ]) &&
167            is_similar( [ map { Im $_ } $left->as_list ],
168                        [ map { Im $_ } $right->as_list ]);
171 sub _not_equal {
172     my ($left, $right) = @_;
173     return !_equal($left,$right);
176 sub copy {
177     my $self = shift;
178     my $copy = Math::GSL::MatrixComplex->new( $self->rows, $self->cols );
179     if ( gsl_matrix_complex_memcpy($copy->raw, $self->raw) != $GSL_SUCCESS ) {
180         croak "Math::GSL - error copying memory, aborting";
181     }
182     return $copy;
184 our %EXPORT_TAGS = ( all => \@EXPORT_OK  );
186 =head1 NAME
188 Math::GSL::MatrixComplex - Complex Matrices
190 =head1 SYNOPSIS
192     use Math::GSL::MatrixComplex qw/:all/;
193     my $matrix1 = Math::GSL::MatrixComplex->new(5,5);  # OO interface
194     my $matrix3 = $matrix1 + $matrix1;
195     my $matrix4 = $matrix1 - $matrix1;
196     if($matrix1 == $matrix4) ...
197     if($matrix1 != $matrix3) ...
199     my $matrix2 = gsl_matrix_complex_alloc(5,5);        # standard interface
201 =head1 Objected Oriented Interface to GSL Math::GSL::MatrixComplex
203 =head2 new()
205 Creates a new MatrixComplex object of the given size.
207     my $matrix = Math::GSL::MatrixComplex->new(10,10);
208 =cut
210 sub new 
212     my ($class, $rows, $cols) = @_;
213     my $this = {}; 
214     my $matrix;
215     if ( defined $rows       && defined $cols && 
216         $rows > 0            && $cols > 0     && 
217         (int $rows == $rows) && (int $cols == $cols)){
219         $matrix  = gsl_matrix_complex_alloc($rows,$cols);
220     } else {
221         croak( __PACKAGE__.'::new($x,$y) - $x and $y must be positive integers');
222     }
223     gsl_matrix_complex_set_zero($matrix);
224     $this->{_matrix} = $matrix; 
225     ($this->{_rows}, $this->{_cols}) = ($rows,$cols);
226     bless $this, $class;
228 =head2 raw()
230 Get the underlying GSL matrix object created by SWIG, useful for using gsl_matrix_* functions which do not have an OO counterpart.
232     my $matrix     = Math::GSL::MatrixComplex->new(3,3);
233     my $gsl_matrix = $matrix->raw;
234     my $stuff      = gsl_matrix_complex_get($gsl_matrix, 1, 2);
236 =cut
237 sub raw  { (shift)->{_matrix} }
238 =head2 rows()
240 Returns the number of rows in the matrix.
242     my $rows = $matrix->rows;
243 =cut
245 sub rows { (shift)->{_rows}   }
247 =head2 cols()
249 Returns the number of columns in the matrix.
251     my $cols = $matrix->cols;
252 =cut
254 sub cols { (shift)->{_cols}   }
256 =head2 as_vector()
258 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:
260     my $vector1 = $matrix->col(0)->as_vector;
261     my $vector2 = $matrix->row(1)->as_vector;
263 =cut
265 sub as_vector($)
267     my ($self) = @_;
268     croak(__PACKAGE__.'::as_vector() - must be a single row or column matrix') unless ($self->cols == 1 or $self->rows == 1);
270     # TODO: there is a faster way to do this
271     return Math::GSL::VectorComplex->new([ $self->as_list ] );
275 =head2  as_list() 
277 Get the contents of a Math::GSL::Matrix object as a Perl list.
279     my $matrix = Math::GSL::MatrixComplex->new(3,3);
280     ...
281     my @matrix = $matrix->as_list;
282 =cut
284 sub as_list($)
286     my $self = shift;
287     my @matrix;
288     for my $row ( 0 .. $self->rows-1) {
289        push @matrix, map { gsl_matrix_complex_get($self->raw, $row, $_) } (0 .. $self->cols-1 );
290     }
291     return map {
292             Math::Complex->make(
293                 gsl_real($_),
294                 gsl_imag($_))
295             } @matrix;
298 =head2 row()
300 Returns a row matrix of the row you enter.
302     my $matrix = Math::GSL::MatrixComplex->new(3,3);
303     ...
304     my $matrix_row = $matrix->row(0);
306 =cut
308 sub row
310     my ($self, $row) = @_;
311     croak (__PACKAGE__.'::$matrix->row($row) - invalid $row value') 
312         unless (($row < $self->rows) and $row >= 0);  
314    my $rowmat = Math::GSL::MatrixComplex->new(1,$self->cols);
316    for my $n (0 .. $self->cols-1) {
317         gsl_matrix_complex_set( $rowmat->raw, 0, $n,
318             gsl_matrix_complex_get($self->raw, $row, $n)
319         );
320    }
322     return $rowmat;
325 =head2 col()
327 Returns a col matrix of the column you enter.
329     my $matrix = Math::GSL::MatrixComplex->new(3,3);
330     ...
331     my $matrix_col = $matrix->col(0);
333 =cut
335 sub col
337     my ($self, $col) = @_;
338     croak (__PACKAGE__."::\$matrix->col(\$col) - $col not a valid column")
339         unless (defined $col && $col < $self->cols and $col >= 0);
341     my $colmat = Math::GSL::MatrixComplex->new($self->rows, 1);
343     map { gsl_matrix_complex_set($colmat->raw, $_, 0,
344             gsl_matrix_complex_get($self->raw, $_, $col),
345             )
346         } (0 .. $self->rows-1);
348     return $colmat;
351 =head2 set_row()
353 Sets a the values of a row with the elements of an array.
355     my $matrix = Math::GSL::MatrixComplex->new(3,3);
356     $matrix->set_row(0, [8, 6, 2]);
358 You can also set multiple rows at once with chained calls:
359     my $matrix = Math::GSL::MatrixComplex->new(3,3);
360     $matrix->set_row(0, [8, 6, 2])
361            ->set_row(1, [2, 4, 1]);
362     ...
364 =cut
366 sub set_row {
367     my ($self, $row, $values) = @_;
368     my $length = $#$values;
369     die __PACKAGE__.'::set_row($x, $values) - $values must be a nonempty array reference' if $length == -1;
370     die __PACKAGE__.'::set_row($x, $values) - $x must be a valid row number' if ($row < 0 || $row >= $self->rows);
371     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);
372     # $values may have Math::Complex objects
373     @$values = map { Math::GSL::Complex::gsl_complex_rect(Re($_), Im($_)) } @$values;
374     # warn Dumper [ @$values ];
375     map { gsl_matrix_complex_set($self->raw, $row, $_, $values->[$_]) } (0..$length);
376     return $self;
379 =head2 set_col()
381 Sets a the values of a column with the elements of an array.
383     my $matrix = Math::GSL::MatrixComplex->new(3,3);
384     $matrix->set_col(0, [8, 6, 2]);
386 You can also set multiple columns at once with chained calls:
387     my $matrix = Math::GSL::MatrixComplex->new(3,3);
388     $matrix->set_col(0, [8, 6, 2])
389            ->set_col(1, [2, 4, 1]);
390     ...
392 =cut
394 sub set_col {
395     my ($self, $col, $values) = @_;
396     my $length = $#$values;
397     die __PACKAGE__.'::set_col($x, $values) - $values must be a nonempty array reference' if $length == -1;
398     die __PACKAGE__.'::set_col($x, $values) - $x must be a valid column number' if ($col < 0 || $col >= $self->cols);
399     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);
400     # $values may have Math::Complex objects
401     @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values;
402     map { gsl_matrix_complex_set($self->raw, $_, $col, $values->[$_]) } (0..$length);
403     return $self;
406 =head2 is_square()
408 Returns true if a matrix is square, i.e. it has the same number of rows as columns, false otherwise.
410 =cut
412 sub is_square($)
414     my $self = shift;
415     return ($self->rows == $self->cols) ? 1 : 0 ;
418 =head2 det()
420 Returns the determinant of a matrix (computed by LU decomposition) or dies if called on a non-square matrix.
422     my $det = $matrix->det();
424 =cut
426 sub det($)
428     my $self = shift;
429     croak(__PACKAGE__."- determinant only exists for square matrices") unless $self->is_square;
431     my $p  = Math::GSL::Permutation->new( $self->rows );
432     my $LU = $self->copy;
434     my $s = gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);
436     # $z is a gsl_complex
437     my $z = gsl_linalg_complex_LU_det($LU->raw, $s );
438     return Math::Complex->make( gsl_real($z), gsl_imag($z) );
441 =head2 zero()
443 Set a matrix to the zero matrix.
445     $matrix->zero;
447 =cut
449 sub zero # brrr!
451     my $self=shift;
452     gsl_matrix_complex_set_zero($self->raw);
453     return $self;
456 =head2 identity()
458 Set a matrix to the identity matrix, i.e. one on the diagonal and zero elsewhere.
460     my $I = $matrix->identity;
462 =cut
464 sub identity
466     my $self=shift;
467     gsl_matrix_complex_set_identity($self->raw);
468     return $self;
471 =head2 inverse()
473 Returns the inverse of a matrix or dies when called on a non-square matrix.
475     my $inverse = $matrix->inverse;
477 =cut
479 sub inverse($)
481     my $self = shift;
482     croak(__PACKAGE__."- inverse only exists for square matrices") unless $self->is_square;
484     my $p  = Math::GSL::Permutation->new( $self->rows );
485     my $LU = $self->copy;
486     my $inverse = $self->copy;
488     # should check return status
489     gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);
490     gsl_linalg_complex_LU_invert($LU->raw, $p->raw,$inverse->raw);
492     return $inverse;
496 =head2 is_hermitian()
498 Returns true if the matrix is hermitian, false otherwise
500     my $test = $matrix->is_hermitian;
502 =cut
504 sub is_hermitian()
506     my ($self) = @_;
507     my $test = $self->is_square;
508     my $transpose = $self->copy;
509     gsl_matrix_complex_transpose($transpose->raw);
510     
511     if($test == 1) {
512         for my $row (0..$self->rows - 1) {
513             map { gsl_matrix_complex_set($transpose->raw, $row, $_, gsl_complex_conjugate(gsl_matrix_complex_get($transpose->raw, $row, $_))) } (0..$self->cols - 1);
514         }
515         if($self != $transpose){
516             $test = 0;
517         }
518     }
519     return $test;
521 =head2 eigenvalues()
523 =cut
526 sub eigenvalues($)
528     my $self=shift;
529     my ($r,$c) = ($self->rows,$self->cols);
530     croak "Math::GSL::MatrixComplex : \$matrix->eigenvalues - \$matrix must be square" unless ($r == $c);
531     my $vector = Math::GSL::Vector->new($r);
532     my $eigen  = gsl_eigen_herm_alloc($r);
535     # GSL has no generalized complex matrix routines
536     # can only continue if the matrix is hermitian
537     croak (__PACKAGE__."::eigenvalues : non-hermitian matrices are not currently supported")  unless $self->is_hermitian;
538     gsl_eigen_herm($self->raw, $vector->raw, $eigen);
540     return $vector->as_list;
543 =head2 lndet()
545 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.
547     my $lndet = $matrix->lndet();
549 =cut
551 sub lndet($)
553     my $self = shift;
554     croak(__PACKAGE__."- log determinant only exists for square matrices") unless $self->is_square;
556     my $p  = Math::GSL::Permutation->new( $self->rows );
557     my $LU = $self->copy;
559     gsl_linalg_complex_LU_decomp($LU->raw, $p->raw);
560     return gsl_linalg_complex_LU_lndet($LU->raw);
562 =head1 DESCRIPTION
564 =over 1
566 =item C<gsl_matrix_complex_alloc >
568 =item C<gsl_matrix_complex_calloc >
570 =item C<gsl_matrix_complex_alloc_from_block >
572 =item C<gsl_matrix_complex_alloc_from_matrix >
574 =item C<gsl_vector_complex_alloc_row_from_matrix >
576 =item C<gsl_vector_complex_alloc_col_from_matrix >
578 =item C<gsl_matrix_complex_free >
580 =item C<gsl_matrix_complex_submatrix >
582 =item C<gsl_matrix_complex_row >
584 =item C<gsl_matrix_complex_column >
586 =item C<gsl_matrix_complex_diagonal >
588 =item C<gsl_matrix_complex_subdiagonal >
590 =item C<gsl_matrix_complex_superdiagonal >
592 =item C<gsl_matrix_complex_subrow >
594 =item C<gsl_matrix_complex_subcolumn >
596 =item C<gsl_matrix_complex_view_array >
598 =item C<gsl_matrix_complex_view_array_with_tda >
600 =item C<gsl_matrix_complex_view_vector >
602 =item C<gsl_matrix_complex_view_vector_with_tda >
604 =item C<gsl_matrix_complex_const_submatrix >
606 =item C<gsl_matrix_complex_const_row >
608 =item C<gsl_matrix_complex_const_column >
610 =item C<gsl_matrix_complex_const_diagonal >
612 =item C<gsl_matrix_complex_const_subdiagonal >
614 =item C<gsl_matrix_complex_const_superdiagonal >
616 =item C<gsl_matrix_complex_const_subrow >
618 =item C<gsl_matrix_complex_const_subcolumn >
620 =item C<gsl_matrix_complex_const_view_array >
622 =item C<gsl_matrix_complex_const_view_array_with_tda >
624 =item C<gsl_matrix_complex_const_view_vector >
626 =item C<gsl_matrix_complex_const_view_vector_with_tda >
628 =item C<gsl_matrix_complex_get>
630 =item C<gsl_matrix_complex_set>
632 =item C<gsl_matrix_complex_ptr>
634 =item C<gsl_matrix_complex_const_ptr>
636 =item C<gsl_matrix_complex_set_zero >
638 =item C<gsl_matrix_complex_set_identity >
640 =item C<gsl_matrix_complex_set_all >
642 =item C<gsl_matrix_complex_fread >
644 =item C<gsl_matrix_complex_fwrite >
646 =item C<gsl_matrix_complex_fscanf >
648 =item C<gsl_matrix_complex_fprintf >
650 =item C<gsl_matrix_complex_memcpy>
652 =item C<gsl_matrix_complex_swap>
654 =item C<gsl_matrix_complex_swap_rows>
656 =item C<gsl_matrix_complex_swap_columns>
658 =item C<gsl_matrix_complex_swap_rowcol>
660 =item C<gsl_matrix_complex_transpose >
662 =item C<gsl_matrix_complex_transpose_memcpy >
664 =item C<gsl_matrix_complex_isnull >
666 =item C<gsl_matrix_complex_ispos >
668 =item C<gsl_matrix_complex_isneg >
670 =item C<gsl_matrix_complex_add >
672 =item C<gsl_matrix_complex_sub >
674 =item C<gsl_matrix_complex_mul_elements >
676 =item C<gsl_matrix_complex_div_elements >
678 =item C<gsl_matrix_complex_scale >
680 =item C<gsl_matrix_complex_add_constant >
682 =item C<gsl_matrix_complex_add_diagonal >
684 =item C<gsl_matrix_complex_get_row>
686 =item C<gsl_matrix_complex_get_col>
688 =item C<gsl_matrix_complex_set_row>
690 =item C<gsl_matrix_complex_set_col>
692 =back 
694 For more informations on the functions, we refer you to the GSL offcial documentation
695 L<http://www.gnu.org/software/gsl/manual/html_node/>
698 =head1 AUTHORS
700 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
702 =head1 COPYRIGHT AND LICENSE
704 Copyright (C) 2008-2009 Jonathan Leto and Thierry Moisan
706 This program is free software; you can redistribute it and/or modify it
707 under the same terms as Perl itself.
709 =cut