5 use Scalar::Util 'blessed';
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
27 '*' => \&_multiplication,
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);
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";
113 gsl_matrix_complex_scale($lcopy->raw, $right);
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);
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);
138 croak "Math::GSL - addition of matrices must be called with two objects matrices and must have the same number of columns and rows";
141 gsl_matrix_complex_add_constant($lcopy->raw, $right);
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);
155 croak "Math::GSL - subtraction of matrices must be called with two objects matrices and must have the same number of columns and rows";
158 gsl_matrix_complex_add_constant($lcopy->raw, $right*-1);
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 ]);
172 my ($left, $right) = @_;
173 return !_equal($left,$right);
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";
184 our %EXPORT_TAGS = ( all => \@EXPORT_OK );
188 Math::GSL::MatrixComplex - Complex Matrices
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
205 Creates a new MatrixComplex object of the given size.
207 my $matrix = Math::GSL::MatrixComplex->new(10,10);
212 my ($class, $rows, $cols) = @_;
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);
221 croak( __PACKAGE__.'::new($x,$y) - $x and $y must be positive integers');
223 gsl_matrix_complex_set_zero($matrix);
224 $this->{_matrix} = $matrix;
225 ($this->{_rows}, $this->{_cols}) = ($rows,$cols);
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);
237 sub raw { (shift)->{_matrix} }
240 Returns the number of rows in the matrix.
242 my $rows = $matrix->rows;
245 sub rows { (shift)->{_rows} }
249 Returns the number of columns in the matrix.
251 my $cols = $matrix->cols;
254 sub cols { (shift)->{_cols} }
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;
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 ] );
277 Get the contents of a Math::GSL::Matrix object as a Perl list.
279 my $matrix = Math::GSL::MatrixComplex->new(3,3);
281 my @matrix = $matrix->as_list;
288 for my $row ( 0 .. $self->rows-1) {
289 push @matrix, map { gsl_matrix_complex_get($self->raw, $row, $_) } (0 .. $self->cols-1 );
300 Returns a row matrix of the row you enter.
302 my $matrix = Math::GSL::MatrixComplex->new(3,3);
304 my $matrix_row = $matrix->row(0);
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)
327 Returns a col matrix of the column you enter.
329 my $matrix = Math::GSL::MatrixComplex->new(3,3);
331 my $matrix_col = $matrix->col(0);
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),
346 } (0 .. $self->rows-1);
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]);
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);
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]);
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);
408 Returns true if a matrix is square, i.e. it has the same number of rows as columns, false otherwise.
415 return ($self->rows == $self->cols) ? 1 : 0 ;
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();
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) );
443 Set a matrix to the zero matrix.
452 gsl_matrix_complex_set_zero($self->raw);
458 Set a matrix to the identity matrix, i.e. one on the diagonal and zero elsewhere.
460 my $I = $matrix->identity;
467 gsl_matrix_complex_set_identity($self->raw);
473 Returns the inverse of a matrix or dies when called on a non-square matrix.
475 my $inverse = $matrix->inverse;
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);
496 =head2 is_hermitian()
498 Returns true if the matrix is hermitian, false otherwise
500 my $test = $matrix->is_hermitian;
507 my $test = $self->is_square;
508 my $transpose = $self->copy;
509 gsl_matrix_complex_transpose($transpose->raw);
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);
515 if($self != $transpose){
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;
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();
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);
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>
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/>
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.