1%perlcode %{ 2 3use strict; 4use Carp qw/croak/; 5use Scalar::Util 'blessed'; 6use Math::Complex; 7use Data::Dumper; 8 9use Math::GSL qw/:all/; 10use Math::GSL::Errno qw/:all/; 11use Math::GSL::Eigen qw/:all/; 12use Math::GSL::Test qw/is_similar/; 13use Math::GSL::BLAS qw/gsl_blas_zgemm/; 14use Math::GSL::CBLAS qw/$CblasNoTrans/; 15use Math::GSL::Complex qw/:all/; 16use Math::GSL::Permutation; 17use Math::GSL::VectorComplex qw/:all/; 18 19use 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 24/; 25 26use overload 27 '*' => \&_multiplication, 28 '+' => \&_addition, 29 '-' => \&_subtract, 30 '==' => \&_equal, 31 '!=' => \&_not_equal, 32# 'abs' => \&_abs, 33 fallback => 1; 34 35our @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 99/; 100 101sub _multiplication { 102 my ($left,$right) = @_; 103 my $lcopy = $left->copy; 104 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; 116} 117 118sub _dot_product { 119 my ($left,$right) = @_; 120 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; 127} 128 129sub _addition { 130 my ($left, $right) = @_; 131 132 my $lcopy = $left->copy; 133 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; 144} 145 146sub _subtract { 147 my ($left, $right) = @_; 148 149 my $lcopy = $left->copy; 150 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; 161} 162 163sub _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 ]); 169} 170 171sub _not_equal { 172 my ($left, $right) = @_; 173 return !_equal($left,$right); 174} 175 176sub 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; 183} 184our %EXPORT_TAGS = ( all => \@EXPORT_OK ); 185 186=encoding utf8 187 188=head1 NAME 189 190Math::GSL::MatrixComplex - Complex Matrices 191 192=head1 SYNOPSIS 193 194 use Math::GSL::MatrixComplex qw/:all/; 195 my $matrix1 = Math::GSL::MatrixComplex->new(5,5); # OO interface 196 my $matrix3 = $matrix1 + $matrix1; 197 my $matrix4 = $matrix1 - $matrix1; 198 if($matrix1 == $matrix4) ... 199 if($matrix1 != $matrix3) ... 200 201 my $matrix2 = gsl_matrix_complex_alloc(5,5); # standard interface 202 203=head1 Objected Oriented Interface to GSL Math::GSL::MatrixComplex 204 205=head2 new() 206 207Creates a new MatrixComplex object of the given size. 208 209 my $matrix = Math::GSL::MatrixComplex->new(10,10); 210=cut 211 212sub new 213{ 214 my ($class, $rows, $cols) = @_; 215 my $this = {}; 216 my $matrix; 217 if ( defined $rows && defined $cols && 218 $rows > 0 && $cols > 0 && 219 (int $rows == $rows) && (int $cols == $cols)){ 220 221 $matrix = gsl_matrix_complex_alloc($rows,$cols); 222 } else { 223 croak( __PACKAGE__.'::new($x,$y) - $x and $y must be positive integers'); 224 } 225 gsl_matrix_complex_set_zero($matrix); 226 $this->{_matrix} = $matrix; 227 ($this->{_rows}, $this->{_cols}) = ($rows,$cols); 228 bless $this, $class; 229} 230=head2 raw() 231 232Get the underlying GSL matrix object created by SWIG, useful for using gsl_matrix_* functions which do not have an OO counterpart. 233 234 my $matrix = Math::GSL::MatrixComplex->new(3,3); 235 my $gsl_matrix = $matrix->raw; 236 my $stuff = gsl_matrix_complex_get($gsl_matrix, 1, 2); 237 238=cut 239sub raw { (shift)->{_matrix} } 240=head2 rows() 241 242Returns the number of rows in the matrix. 243 244 my $rows = $matrix->rows; 245=cut 246 247sub rows { (shift)->{_rows} } 248 249=head2 cols() 250 251Returns the number of columns in the matrix. 252 253 my $cols = $matrix->cols; 254=cut 255 256sub cols { (shift)->{_cols} } 257 258=head2 as_vector() 259 260Returns 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: 261 262 my $vector1 = $matrix->col(0)->as_vector; 263 my $vector2 = $matrix->row(1)->as_vector; 264 265=cut 266 267sub as_vector($) 268{ 269 my ($self) = @_; 270 croak(__PACKAGE__.'::as_vector() - must be a single row or column matrix') unless ($self->cols == 1 or $self->rows == 1); 271 272 # TODO: there is a faster way to do this 273 return Math::GSL::VectorComplex->new([ $self->as_list ] ); 274 275} 276 277=head2 as_list() 278 279Get the contents of a Math::GSL::Matrix object as a Perl list. 280 281 my $matrix = Math::GSL::MatrixComplex->new(3,3); 282 ... 283 my @matrix = $matrix->as_list; 284=cut 285 286sub as_list($) 287{ 288 my $self = shift; 289 my @matrix; 290 for my $row ( 0 .. $self->rows-1) { 291 push @matrix, map { gsl_matrix_complex_get($self->raw, $row, $_) } (0 .. $self->cols-1 ); 292 } 293 return map { 294 Math::Complex->make( 295 gsl_real($_), 296 gsl_imag($_)) 297 } @matrix; 298} 299 300=head2 row() 301 302Returns a row matrix of the row you enter. 303 304 my $matrix = Math::GSL::MatrixComplex->new(3,3); 305 ... 306 my $matrix_row = $matrix->row(0); 307 308=cut 309 310sub row 311{ 312 my ($self, $row) = @_; 313 croak (__PACKAGE__.'::$matrix->row($row) - invalid $row value') 314 unless (($row < $self->rows) and $row >= 0); 315 316 my $rowmat = Math::GSL::MatrixComplex->new(1,$self->cols); 317 318 for my $n (0 .. $self->cols-1) { 319 gsl_matrix_complex_set( $rowmat->raw, 0, $n, 320 gsl_matrix_complex_get($self->raw, $row, $n) 321 ); 322 } 323 324 return $rowmat; 325} 326 327=head2 col() 328 329Returns a col matrix of the column you enter. 330 331 my $matrix = Math::GSL::MatrixComplex->new(3,3); 332 ... 333 my $matrix_col = $matrix->col(0); 334 335=cut 336 337sub col 338{ 339 my ($self, $col) = @_; 340 croak (__PACKAGE__."::\$matrix->col(\$col) - $col not a valid column") 341 unless (defined $col && $col < $self->cols and $col >= 0); 342 343 my $colmat = Math::GSL::MatrixComplex->new($self->rows, 1); 344 345 map { gsl_matrix_complex_set($colmat->raw, $_, 0, 346 gsl_matrix_complex_get($self->raw, $_, $col), 347 ) 348 } (0 .. $self->rows-1); 349 350 return $colmat; 351} 352 353=head2 set_row() 354 355Sets a the values of a row with the elements of an array. 356 357 my $matrix = Math::GSL::MatrixComplex->new(3,3); 358 $matrix->set_row(0, [8, 6, 2]); 359 360You can also set multiple rows at once with chained calls: 361 my $matrix = Math::GSL::MatrixComplex->new(3,3); 362 $matrix->set_row(0, [8, 6, 2]) 363 ->set_row(1, [2, 4, 1]); 364 ... 365 366=cut 367 368sub set_row { 369 my ($self, $row, $values) = @_; 370 my $length = $#$values; 371 die __PACKAGE__.'::set_row($x, $values) - $values must be a nonempty array reference' if $length == -1; 372 die __PACKAGE__.'::set_row($x, $values) - $x must be a valid row number' if ($row < 0 || $row >= $self->rows); 373 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); 374 # $values may have Math::Complex objects 375 @$values = map { Math::GSL::Complex::gsl_complex_rect(Re($_), Im($_)) } @$values; 376 # warn Dumper [ @$values ]; 377 map { gsl_matrix_complex_set($self->raw, $row, $_, $values->[$_]) } (0..$length); 378 return $self; 379} 380 381=head2 set_col() 382 383Sets a the values of a column with the elements of an array. 384 385 my $matrix = Math::GSL::MatrixComplex->new(3,3); 386 $matrix->set_col(0, [8, 6, 2]); 387 388You can also set multiple columns at once with chained calls: 389 my $matrix = Math::GSL::MatrixComplex->new(3,3); 390 $matrix->set_col(0, [8, 6, 2]) 391 ->set_col(1, [2, 4, 1]); 392 ... 393 394=cut 395 396sub set_col { 397 my ($self, $col, $values) = @_; 398 my $length = $#$values; 399 die __PACKAGE__.'::set_col($x, $values) - $values must be a nonempty array reference' if $length == -1; 400 die __PACKAGE__.'::set_col($x, $values) - $x must be a valid column number' if ($col < 0 || $col >= $self->cols); 401 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); 402 # $values may have Math::Complex objects 403 @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values; 404 map { gsl_matrix_complex_set($self->raw, $_, $col, $values->[$_]) } (0..$length); 405 return $self; 406} 407 408=head2 is_square() 409 410Returns true if a matrix is square, i.e. it has the same number of rows as columns, false otherwise. 411 412=cut 413 414sub is_square($) 415{ 416 my $self = shift; 417 return ($self->rows == $self->cols) ? 1 : 0 ; 418} 419 420=head2 det() 421 422Returns the determinant of a matrix (computed by LU decomposition) or dies if called on a non-square matrix. 423 424 my $det = $matrix->det(); 425 426=cut 427 428sub det($) 429{ 430 my $self = shift; 431 croak(__PACKAGE__."- determinant only exists for square matrices") unless $self->is_square; 432 433 my $p = Math::GSL::Permutation->new( $self->rows ); 434 my $LU = $self->copy; 435 436 my $s = gsl_linalg_complex_LU_decomp($LU->raw, $p->raw); 437 438 # $z is a gsl_complex 439 my $z = gsl_linalg_complex_LU_det($LU->raw, $s ); 440 return Math::Complex->make( gsl_real($z), gsl_imag($z) ); 441} 442 443=head2 zero() 444 445Set a matrix to the zero matrix. 446 447 $matrix->zero; 448 449=cut 450 451sub zero # brrr! 452{ 453 my $self=shift; 454 gsl_matrix_complex_set_zero($self->raw); 455 return $self; 456} 457 458=head2 identity() 459 460Set a matrix to the identity matrix, i.e. one on the diagonal and zero elsewhere. 461 462 my $I = $matrix->identity; 463 464=cut 465 466sub identity 467{ 468 my $self=shift; 469 gsl_matrix_complex_set_identity($self->raw); 470 return $self; 471} 472 473=head2 inverse() 474 475Returns the inverse of a matrix or dies when called on a non-square matrix. 476 477 my $inverse = $matrix->inverse; 478 479=cut 480 481sub inverse($) 482{ 483 my $self = shift; 484 croak(__PACKAGE__."- inverse only exists for square matrices") unless $self->is_square; 485 486 my $p = Math::GSL::Permutation->new( $self->rows ); 487 my $LU = $self->copy; 488 my $inverse = $self->copy; 489 490 # should check return status 491 gsl_linalg_complex_LU_decomp($LU->raw, $p->raw); 492 gsl_linalg_complex_LU_invert($LU->raw, $p->raw,$inverse->raw); 493 494 return $inverse; 495 496} 497 498=head2 is_hermitian() 499 500Returns true if the matrix is hermitian, false otherwise 501 502 my $test = $matrix->is_hermitian; 503 504=cut 505 506sub is_hermitian() 507{ 508 my ($self) = @_; 509 my $test = $self->is_square; 510 my $transpose = $self->copy; 511 gsl_matrix_complex_transpose($transpose->raw); 512 513 if($test == 1) { 514 for my $row (0..$self->rows - 1) { 515 map { gsl_matrix_complex_set($transpose->raw, $row, $_, gsl_complex_conjugate(gsl_matrix_complex_get($transpose->raw, $row, $_))) } (0..$self->cols - 1); 516 } 517 if($self != $transpose){ 518 $test = 0; 519 } 520 } 521 return $test; 522} 523=head2 eigenvalues() 524 525=cut 526 527 528sub eigenvalues($) 529{ 530 my $self=shift; 531 my ($r,$c) = ($self->rows,$self->cols); 532 croak "Math::GSL::MatrixComplex : \$matrix->eigenvalues - \$matrix must be square" unless ($r == $c); 533 my $vector = Math::GSL::Vector->new($r); 534 my $eigen = gsl_eigen_herm_alloc($r); 535 536 537 # GSL has no generalized complex matrix routines 538 # can only continue if the matrix is hermitian 539 croak (__PACKAGE__."::eigenvalues : non-hermitian matrices are not currently supported") unless $self->is_hermitian; 540 gsl_eigen_herm($self->raw, $vector->raw, $eigen); 541 542 return $vector->as_list; 543} 544 545=head2 lndet() 546 547Returns 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. 548 549 my $lndet = $matrix->lndet(); 550 551=cut 552 553sub lndet($) 554{ 555 my $self = shift; 556 croak(__PACKAGE__."- log determinant only exists for square matrices") unless $self->is_square; 557 558 my $p = Math::GSL::Permutation->new( $self->rows ); 559 my $LU = $self->copy; 560 561 gsl_linalg_complex_LU_decomp($LU->raw, $p->raw); 562 return gsl_linalg_complex_LU_lndet($LU->raw); 563} 564=head1 DESCRIPTION 565 566=over 1 567 568=item C<gsl_matrix_complex_alloc > 569 570=item C<gsl_matrix_complex_calloc > 571 572=item C<gsl_matrix_complex_alloc_from_block > 573 574=item C<gsl_matrix_complex_alloc_from_matrix > 575 576=item C<gsl_vector_complex_alloc_row_from_matrix > 577 578=item C<gsl_vector_complex_alloc_col_from_matrix > 579 580=item C<gsl_matrix_complex_free > 581 582=item C<gsl_matrix_complex_submatrix > 583 584=item C<gsl_matrix_complex_row > 585 586=item C<gsl_matrix_complex_column > 587 588=item C<gsl_matrix_complex_diagonal > 589 590=item C<gsl_matrix_complex_subdiagonal > 591 592=item C<gsl_matrix_complex_superdiagonal > 593 594=item C<gsl_matrix_complex_subrow > 595 596=item C<gsl_matrix_complex_subcolumn > 597 598=item C<gsl_matrix_complex_view_array > 599 600=item C<gsl_matrix_complex_view_array_with_tda > 601 602=item C<gsl_matrix_complex_view_vector > 603 604=item C<gsl_matrix_complex_view_vector_with_tda > 605 606=item C<gsl_matrix_complex_const_submatrix > 607 608=item C<gsl_matrix_complex_const_row > 609 610=item C<gsl_matrix_complex_const_column > 611 612=item C<gsl_matrix_complex_const_diagonal > 613 614=item C<gsl_matrix_complex_const_subdiagonal > 615 616=item C<gsl_matrix_complex_const_superdiagonal > 617 618=item C<gsl_matrix_complex_const_subrow > 619 620=item C<gsl_matrix_complex_const_subcolumn > 621 622=item C<gsl_matrix_complex_const_view_array > 623 624=item C<gsl_matrix_complex_const_view_array_with_tda > 625 626=item C<gsl_matrix_complex_const_view_vector > 627 628=item C<gsl_matrix_complex_const_view_vector_with_tda > 629 630=item C<gsl_matrix_complex_get> 631 632=item C<gsl_matrix_complex_set> 633 634=item C<gsl_matrix_complex_ptr> 635 636=item C<gsl_matrix_complex_const_ptr> 637 638=item C<gsl_matrix_complex_set_zero > 639 640=item C<gsl_matrix_complex_set_identity > 641 642=item C<gsl_matrix_complex_set_all > 643 644=item C<gsl_matrix_complex_fread > 645 646=item C<gsl_matrix_complex_fwrite > 647 648=item C<gsl_matrix_complex_fscanf > 649 650=item C<gsl_matrix_complex_fprintf > 651 652=item C<gsl_matrix_complex_memcpy> 653 654=item C<gsl_matrix_complex_swap> 655 656=item C<gsl_matrix_complex_swap_rows> 657 658=item C<gsl_matrix_complex_swap_columns> 659 660=item C<gsl_matrix_complex_swap_rowcol> 661 662=item C<gsl_matrix_complex_transpose > 663 664=item C<gsl_matrix_complex_transpose_memcpy > 665 666=item C<gsl_matrix_complex_isnull > 667 668=item C<gsl_matrix_complex_ispos > 669 670=item C<gsl_matrix_complex_isneg > 671 672=item C<gsl_matrix_complex_add > 673 674=item C<gsl_matrix_complex_sub > 675 676=item C<gsl_matrix_complex_mul_elements > 677 678=item C<gsl_matrix_complex_div_elements > 679 680=item C<gsl_matrix_complex_scale > 681 682=item C<gsl_matrix_complex_add_constant > 683 684=item C<gsl_matrix_complex_add_diagonal > 685 686=item C<gsl_matrix_complex_get_row> 687 688=item C<gsl_matrix_complex_get_col> 689 690=item C<gsl_matrix_complex_set_row> 691 692=item C<gsl_matrix_complex_set_col> 693 694=back 695 696For more informations on the functions, we refer you to the GSL official documentation 697L<http://www.gnu.org/software/gsl/manual/html_node/> 698 699 700=head1 AUTHORS 701 702Jonathan "Duke" Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com> 703 704=head1 COPYRIGHT AND LICENSE 705 706Copyright (C) 2008-2021 Jonathan "Duke" Leto and Thierry Moisan 707 708This program is free software; you can redistribute it and/or modify it 709under the same terms as Perl itself. 710 711=cut 712 713%} 714 715