Comment out some more free() calls, which need to be turned on per-subsystem to weed...
[Math-GSL.git] / pod / VectorComplex.pod
blob103bd9eb95d630ddb90a8bf23a11b188bd23702d
1 %perlcode %{
2 use Scalar::Util 'blessed';
3 use Data::Dumper;
4 use Carp qw/croak/;
5 use Math::GSL::Errno qw/:all/;
6 use Math::GSL::BLAS qw/gsl_blas_ddot/;
7 use Math::GSL::Complex qw/:all/;
8 use Math::GSL::Test qw/is_status_ok/;
9 use Math::Complex;
10 use overload
11     '*'      => \&_multiplication,
12     '+'      => \&_addition,
13     '-'      => \&_subtract,
14     fallback => 1,
17 @EXPORT_all  = qw/
18                  gsl_vector_complex_alloc gsl_vector_complex_calloc gsl_vector_complex_alloc_from_block gsl_vector_complex_alloc_from_vector
19                  gsl_vector_complex_free gsl_vector_complex_view_array gsl_vector_complex_view_array_with_stride gsl_vector_complex_const_view_array
20                  gsl_vector_complex_const_view_array_with_stride gsl_vector_complex_subvector gsl_vector_complex_subvector_with_stride
21                  gsl_vector_complex_const_subvector gsl_vector_complex_const_subvector_with_stride gsl_vector_complex_real gsl_vector_complex_imag
22                  gsl_vector_complex_const_real gsl_vector_complex_const_imag gsl_vector_complex_get gsl_vector_complex_set
23                  gsl_vector_complex_ptr gsl_vector_complex_const_ptr gsl_vector_complex_set_zero gsl_vector_complex_set_all
24                  gsl_vector_complex_set_basis gsl_vector_complex_fread gsl_vector_complex_fwrite gsl_vector_complex_fscanf
25                  gsl_vector_complex_fprintf gsl_vector_complex_memcpy gsl_vector_complex_reverse gsl_vector_complex_swap
26                  gsl_vector_complex_swap_elements gsl_vector_complex_isnull gsl_vector_complex_ispos gsl_vector_complex_isneg
28 @EXPORT_OK = (@EXPORT_all);
29 %EXPORT_TAGS = ( all => \@EXPORT_all );
31 =head1 NAME
33 Math::GSL::VectorComplex - Complex Vectors
35 =head1 SYNOPSIS
37     use Math::GSL::VectorComplex qw/:all/;
38     my $vec1 = Math::GSL::VectorComplex->new([1 + 2*i, 7*i, 5, -3 ]);
39     my $vec2 = $vec1 * 5;
40     my $vec3 = Math::GSL::Vector>new(10);   # 10 element zero vector 
41     my $vec4 = $vec1 + $vec2;
43     # set the element at index 1 to -i
44     # and the element at index 3 to i
45     $vec3->set([ 1, -i ], [ 9, i ]);
47     my @vec = $vec2->as_list;               # return elements as Perl list
49     my $dot_product = $vec1 * $vec2;
50     my $length      = $vec2->length;
51     my $first       = $vec1->get(0);
54 =cut
56 =head1 Objected Oriented Interface to GSL Math::GSL::VectorComplex
58 =head2 new()
60 Creates a new Vector of the given size.
62     my $vector = Math::GSL::VectorComplex->new(3);
64 You can also create and set directly the values of the vector like this :
66    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
68 =cut
70 sub new {
71     my ($class, $values) = @_;
72     my $length  = $#$values;
73     my $this = {};
74     my $vector;
76     # we expect $values to have Math::Complex objects
77     @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values;
79     if ( ref $values eq 'ARRAY' ){
80         die __PACKAGE__.'::new($x) - $x must be a nonempty array reference' if $length == -1;
81         $vector  = gsl_vector_complex_alloc($length+1);
82         map { gsl_vector_complex_set($vector, $_, $values->[$_] ) }  (0 .. $length);
83         $this->{_length} = $length+1;
84     } elsif ( (int($values) == $values) && ($values > 0)) {
85         $vector  = gsl_vector_complex_alloc($values);
86         gsl_vector_complex_set_zero($vector);
87         $this->{_length} = $values;
88     } else {
89         die __PACKAGE__.'::new($x) - $x must be an int or array reference';
90     }
91     $this->{_vector} = $vector;
92     bless $this, $class;
96 =head2 raw()
98 Get the underlying GSL vector object created by SWIG, useful for using gsl_vector_* functions which do not have an OO counterpart.
100     my $vector    = Math::GSL::VectorComplex->new(3);
101     my $gsl_vector = $vector->raw;
102     my $stuff      = gsl_vector_get($gsl_vector, 1);
104 =cut
106 sub raw { 
107     my $self = shift;
108     return $self->{_vector};
111 =head2 min()
113 Returns the minimum value contained in the vector.
115    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
116    my $minimum = $vector->min;
118 =cut 
120 sub min {
121     my $self=shift;
122     return gsl_vector_min($self->raw);
125 =head2 max()
127 Returns the minimum value contained in the vector.
129    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
130    my $maximum = $vector->max;
132 =cut 
134 sub max {
135     my $self=shift;
136     return gsl_vector_max($self->raw);
139 =head2 length()
141 Returns the number of elements contained in the vector.
143    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
144    my $length = $vector->length;
146 =cut 
148 sub length { my $self=shift; $self->{_length} }
150 =head2  as_list() 
152 Gets the content of a Math::GSL::Vector object as a Perl list.
154     my $vector = Math::GSL::VectorComplex->new(3);
155     ...
156     my @values = $vector->as_list;
157 =cut
159 sub as_list {
160     my $self=shift;
161     # this is wrong
162     return map { cplxe( gsl_complex_abs($_), gsl_complex_arg($_) ) } $self->get( [ 0 .. $self->length - 1  ] );
165 =head2  get()
167 Gets the value of an of a Math::GSL::Vector object.
169     my $vector = Math::GSL::VectorComplex->new(3);
170     ...
171     my @values = $vector->get(2);
173 You can also enter an array of indices to receive their corresponding values:
175     my $vector = Math::GSL::VectorComplex->new(3);
176     ...
177     my @values = $vector->get([0,2]);
179 =cut
181 sub get {
182     my ($self, $indices) = @_;
183     return  map {  gsl_vector_complex_get($self->raw, $_ ) } @$indices ;
186 =head2 reverse()
188 Returns the a vector with the elements in reversed order.
190     use Math::Complex;
191     my $v1 = Math::GSL::VectorComplex->new([ 1, 2, 3*i]);
192     my $v2 = $v1->reverse;
194 =cut
196 sub reverse {
197     my $self = shift;
198     my $copy = $self->copy();
199     unless ( is_status_ok( gsl_vector_complex_reverse( $copy->raw )) ) {
200         die( __PACKAGE__.": error reversing vector " . gsl_strerror($status) );
201     }
202     return $copy;
205 =head2  set() 
207 Sets values of an of a Math::GSL::Vector object.
209     my $vector = Math::GSL::VectorComplex->new(3);
210     $vector->set([1,2], [8,23]);
212 This sets the second and third value to 8 and 23.
214 =cut
216 sub set {
217     my ($self, $indices, $values) = @_;
218     die (__PACKAGE__.'::set($indices, $values) - $indices and $values must be array references of the same length')
219         unless ( ref $indices eq 'ARRAY' && ref $values eq 'ARRAY' &&  $#$indices == $#$values );
220     eval {
221         map {  gsl_vector_complex_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
222     };
223     # better error handling?
224     warn $@ if $@;
225     return;
228 =head2 copy()
230 Returns a copy of the vector, which has the same length and values but resides at a different location in memory.
232     my $vector = Math::GSL::VectorComplex->new([10 .. 20]);
233     my $copy   = $vector->copy;
235 =cut
237 sub copy {
238     my $self = shift;
239     my $copy = Math::GSL::VectorComplex->new( $self->length );
240     my $status = gsl_vector_complex_memcpy($copy->raw, $self->raw);
241     if ( $status != $GSL_SUCCESS ) {
242         croak "Math::GSL - error copying memory, aborting. $! status=$status";
243     }
244     return $copy;
247 =head2 swap()
249 Exchanges the values in the vectors $v with $w by copying.
251     my $v = Math::GSL::VectorComplex->new([1..5]);
252     my $w = Math::GSL::VectorComplex->new([3..7]);
253     $v->swap( $w );
255 =cut
257 sub swap() {
258     my ($self,$other) = @_;
259     croak "Math::GSL::VectorComplex : \$v->swap(\$w) - \$w must be a Math::GSL::VectorComplex"
260         unless ref $other eq 'Math::GSL::VectorComplex';
261     gsl_vector_complex_swap( $self->raw, $other->raw );
262     return $self;
265 sub _multiplication {
266     my ($left,$right) = @_;
267     my $lcopy = $left->copy;
269     if ( blessed $right && $right->isa(__PACKAGE__) ) {
270         return $lcopy->dot_product($right);
271     } else {
272         # will be in upcoming gsl 1.12
273         # gsl_vector_complex_scale($lcopy->raw, $right);
274     }
275     return $lcopy;
278 sub _subtract {
279     my ($left, $right, $flip) = @_;
281     if ($flip) {
282         my $lcopy = $left->copy;
283         # will be in upcoming gsl 1.12
284         # gsl_vector_complex_scale($lcopy->raw, -1 );
285         gsl_vector_add_constant($lcopy->raw, $right);
286         return $lcopy;
287     } else {
288         return _addition($left, -1.0*$right);
289     }
292 sub _addition {
293     my ($left, $right, $flip) = @_;
294     my $lcopy = $left->copy;
296     if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
297         if ( $left->length == $right->length ) {
298             gsl_vector_complex_add($lcopy->raw, $right->raw);
299         } else {
300             croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
301         }
302     } else {
303         gsl_vector_complex_add_constant($lcopy->raw, $right);
304     }
305     return $lcopy;
308 sub dot_product_pp {
309     my ($left,$right) = @_;
310     my $sum=0;
311     if ( blessed $right && $right->isa('Math::GSL::Vector') && $left->length == $right->length ) {
312          my @l = $left->as_list;
313          my @r = $right->as_list;
314          map { $sum += $l[$_] * $r[$_] } (0..$#l);
315         return $sum;
316     } else {
317         croak "dot_product() must be called with two vectors";
318     }
321 sub dot_product {
322     my ($left,$right) = @_;
324     my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
325     croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
326     return $product;
329 =head1 AUTHORS
331 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
333 =head1 COPYRIGHT AND LICENSE
335 Copyright (C) 2008-2009 Jonathan Leto and Thierry Moisan
337 This program is free software; you can redistribute it and/or modify it
338 under the same terms as Perl itself.
340 =cut