diff --git a/MANIFEST b/MANIFEST index 257274fb3..f8aca3959 100644 --- a/MANIFEST +++ b/MANIFEST @@ -197,7 +197,6 @@ t/fits-noafh.t t/fits.t t/flexraw-iotypes.t t/flexraw.t -t/flexraw_fortran.t t/func.t t/image2d.t t/imagend.t diff --git a/lib/PDL/Course.pod b/lib/PDL/Course.pod index e218bf690..017bad4d8 100644 --- a/lib/PDL/Course.pod +++ b/lib/PDL/Course.pod @@ -273,7 +273,7 @@ eigenvalues, inverting square matrices, LU-decomposition, and solving a system of linear equations. Though it is not built on L, it should generally work with that module. Also, the methods provided by this module do not depend on external libraries such as -Slatec or GSL. +LAPACK or GSL. =item * L diff --git a/lib/PDL/IO/FlexRaw.pm b/lib/PDL/IO/FlexRaw.pm index b08cf27fd..5d01c54d3 100644 --- a/lib/PDL/IO/FlexRaw.pm +++ b/lib/PDL/IO/FlexRaw.pm @@ -183,6 +183,76 @@ and return the PDLs mapped before the problem arose. This can be dealt with either by reorganizing the data file (large types first helps, as a rule-of-thumb), or more simply by using C. +=head2 Fortran code to create data + +Until PDL 2.099, the test file F compiled a +Fortran program, ran it, then byte-swapped its output, to test this +module's ability to do that. Version 2.099 has dropped external +dependencies, including the use of Fortran. The code it used is +shown here for historical curiosity: + + c Program to test i/o of F77 unformatted files + program rawtest + implicit none + integer i + $f77type a($ndata) + do i = 1, $ndata + a(i) = $val + enddo + open(8,file= + \$'$data' + \$,status='new',form='unformatted') + i = $ndata + write (8) i + write (8) a + close(8) + end + +with this FlexRaw header: + + # FlexRaw file header + f77 + long 1 1 + # Data + $pdltype 1 $ndata + +C<$ndata> was set to 10, C<$val> was C<100.*sin(0.01* i)>, C<$data> +was a filename. C<$f77type> was set to C and C. + +There was also a more complex program: + + c Program to test i/o of F77 unformatted files + program rawtest + implicit none + character a + integer*2 i + integer*4 l + real*4 f + real*8 d + d = 4*atan(1.) + f = d + l = 10**d + i = l + a = ' ' + open(8,file= + \$'$data' + \$,status='new',form='unformatted') + c Choose bad boundaries... + write (8) a,i,l,f,d + close(8) + end + +with this FlexRaw header: + + # FlexRaw file header + byte 1 4 + byte 0 + short 0 + long 0 + float 0 + double 0 + byte 1 4 + =head1 FUNCTIONS =cut diff --git a/lib/PDL/Index.pod b/lib/PDL/Index.pod index 235fc7d2b..de4baa498 100644 --- a/lib/PDL/Index.pod +++ b/lib/PDL/Index.pod @@ -536,10 +536,6 @@ L - a C function for PDL =item * -L - PDL interface to the slatec numerical programming library - -=item * - L - Indexing, slicing, and dicing =item * diff --git a/lib/PDL/Matrix.pm b/lib/PDL/Matrix.pm index 9d58523e8..866745306 100644 --- a/lib/PDL/Matrix.pm +++ b/lib/PDL/Matrix.pm @@ -29,7 +29,7 @@ vectors resp and matrices can be typed in using traditional matrix convention. If you want to know more about matrix operation support in PDL, you -want to read L or L. +want to read L or L. The original pdl class refers to the first index as the first row, the second index as the first column of a matrix. Consider diff --git a/lib/PDL/MatrixOps.pd b/lib/PDL/MatrixOps.pd index d9e178339..f7aad7e44 100644 --- a/lib/PDL/MatrixOps.pd +++ b/lib/PDL/MatrixOps.pd @@ -37,8 +37,8 @@ contains utilities for many common matrix operations: inversion, determinant finding, eigenvalue/vector finding, singular value decomposition, etc. PDL::MatrixOps routines are written in a mixture of Perl and C, so that they are reliably present even when there is no -FORTRAN compiler or external library available (e.g. -L or any of the PDL::GSL family of modules). +external library available (e.g. +L or any of the PDL::GSL family of modules). Matrix manipulation, particularly with large matrices, is a challenging field and no one algorithm is suitable in all cases. The @@ -1030,13 +1030,6 @@ Solve A x = B for matrix A, by back substitution into A's LU decomposition. ($x) = simq($A->dupN(1,1,map +($A_dims[$_]//1)==1?$broadcast[$_]:1, 0..$#broadcast)->copy, $B->transpose, 0); $x = $x->inplace->transpose; - # or with Slatec LINPACK - use PDL::Slatec; - gefa($lu=$A->copy, $ipiv=null, $info=null); - # 1 = do transpose because Fortran's idea of rows vs columns - gesl($lu, $ipiv, $x=$B->transpose->copy, 1); - $x = $x->inplace->transpose; - # or with PDL::LinearAlgebra wrappers of LAPACK $x = msolve($A, $B); diff --git a/lib/PDL/Modules.pod b/lib/PDL/Modules.pod index 9ab9e4d0a..9026b44d9 100644 --- a/lib/PDL/Modules.pod +++ b/lib/PDL/Modules.pod @@ -169,14 +169,9 @@ Linear filtering. Simplex optimization routines. -=item L - -PDL interface to the Slatec library. - =back - =head1 COORDINATE TRANSFORMATIONS =over 5 diff --git a/lib/PDL/PP.pod b/lib/PDL/PP.pod index ea94bdb8e..3656ff8e7 100644 --- a/lib/PDL/PP.pod +++ b/lib/PDL/PP.pod @@ -1316,8 +1316,8 @@ further. Perl can help you again(!) in doing this. In many libraries you have certain calling conventions. This can be exploited. In short, you can write a little parser (which is really not difficult in Perl) that then generates the calls to C from parsed descriptions of the -functions in that library. For an example, please check the I -interface in the C tree of the PDL distribution. If you want to check +functions in that library. For an example, please check the L +interface on CPAN. If you want to check (during debugging) which calls to PP functions your Perl code generated then set C<$::PP_VERBOSE> to a true value just before the relevant C call, or at the top of the file. diff --git a/t/flexraw_fortran.t b/t/flexraw_fortran.t deleted file mode 100644 index 7923c257f..000000000 --- a/t/flexraw_fortran.t +++ /dev/null @@ -1,619 +0,0 @@ -use strict; -use warnings; -use PDL::LiteF; -use PDL::IO::FlexRaw; -use Config; -use Test::More; -use Test::PDL; -use File::Temp qw(tempfile); -use File::Spec; -use File::Which (); - -my $prog_iter; - -(undef, my $data) = tempfile("rawXXXX", SUFFIX=>'_data', TMPDIR=>1); -my $hdr = $data . '.hdr'; -(my $head = $data) =~ s/_data$//; - -$|=1; - -my $ndata = 10; -my $Verbose = 0; -my $DEBUG = 0; -$PDL::Verbose = 0; -$Verbose |= $PDL::Verbose; - -my $null = ' 2>' . File::Spec->devnull; - -my $datalen = length($data); -eval "use PDL::Slatec"; -plan skip_all => "Skipped tests as no Slatec" if $@; -plan skip_all => "temp file path too long for f77 ($datalen chars), skipping all tests" - if $datalen > 70; -eval "use ExtUtils::F77"; -plan skip_all => "Skip all tests as ExtUtils::F77 not found" if $@; - -my $F77; -my $F77flags; - -if ($ExtUtils::F77::VERSION > 1.03) { - $F77 = ExtUtils::F77::compiler(); - $F77flags = ExtUtils::F77::cflags(); -} else { - $F77 = 'f77'; - $F77flags = ''; -} - -sub byte4swap { - my ($file) = @_; - my ($ofile) = $file.'~'; - my ($word); - - open my $ifh, "<", $file or die "Can't open $file to read: $!"; - open my $ofh, ">", $ofile or die "Can't open $ofile to write: $!"; - binmode $ifh; - binmode $ofh; - while ( !eof $ifh ) { - read $ifh, $word, 4; - $word = pack 'c4',reverse unpack 'c4',$word; - print $ofh $word; - } - close $ofh; - close $ifh; - rename $ofile, $file; -} - -sub byte8swap { - my ($file) = @_; - my ($ofile) = $file.'~'; - my ($word); - - open my $ifh, "<", $file or die "Can't open $file to read: $!"; - open my $ofh, ">", $ofile or die "Can't open $ofile to write: $!"; - binmode $ifh; - binmode $ofh; - while ( !$ifh->eof ) { - read $ifh, $word, 8; - $word = pack 'c8',reverse unpack 'c8',$word; - print $ofh $word; - } - $ofh->close; - $ifh->close; - rename $ofile, $file; -} - -# utility to fold long lines preventing problems with 72 -# char limit and long text parameters (e.g. filenames) -sub codefold { - my $oldcode = shift; - my $newcode = ''; - - eval { - # to simplify loop processing, introduces dependence - require IO::String; - - my $in = IO::String->new($oldcode); - my $out = IO::String->new($newcode); - - # find non-comment lines longer than 72 columns and fold - my $line = ''; - while ($line = <$in>) { - - # clean off line-feed stuff - chomp $line; - - # pass comments to output - print $out "$line\n" if $line =~ /^\S/; - - # output code lines (by 72-char chunks if needed) - while ($line ne '') { - - # output first 72 columns of the line - print $out substr($line,0,72) . "\n"; # print first 72 cols - - if ( length($line) > 72 ) { - # make continuation line of the rest of the line - substr($line,0,72) = ' $'; - } else { - $line = ''; - } - - } - } - - # close "files" and return folded code - close($in); - close($out); - }; - - $newcode = $oldcode if $@; - return $newcode; -} - -# createData $head, $code -# -# given a F77 program (in $code), compile and run it. -# It is expected to create a data file called -# "${head}data" -# The executable and code are cleaned up but NOT the -# data file. -# -# Requires the global variables -# $F77 -# $F77flags -# $Verbose -# $DEBUG -# -sub createData { - my $head = shift; - my $code = shift; - - # try and provide a modicum of safety, since we call - # system with $head as the argument - # - if($^O =~ /mswin32/i) { - die '$head [' . $head . '] should match /^[A-Z]:/' - unless $head =~ /^[A-Z]:/; - } - else { - die '$head [' . $head . '] must start with a / or ./' - unless $head =~ /^(\/|\.\/)/; - } - - my $file = ${head} . '.f'; - my $prog = $head; - - open my $fh, ">", $file - or die "ERROR: Unable to write F77 code to $file: $!\n"; - print $fh $code; - close $fh; - - if (-e "$prog$Config{_exe}") { - my $success = unlink "$prog$Config{_exe}"; - if (!$success) { - # Antivirus might be holding on to the exec - # file so use a different name - # Poss path length issues? - # See $dalen in 2nd begin block$prog_iter++; - warn "Unable to delete $prog$Config{_exe}, generating new name" - if $Verbose; - $prog_iter++; - $prog .= $prog_iter; - }; - } - - system("$F77 $F77flags -o $prog$Config{_exe} $file". - (($Verbose || $DEBUG)?'': $null)); - - unlink $data if -f $data; - system( $prog ); - - die "ERROR: code did not create data file $data\n" - unless -e $data; - - unlink $prog.$Config{_exe}, $file; - -} # sub: createData() - -# Types to test the translation for, perl + f77 forms -my %types = ( 'float' => 'real*4', 'double' => 'real*8', 'long' => 'integer*4', - 'short' => 'integer*2', 'byte' => 'character' ); - -# Perl and f77 functions should be have the same net effect... -my $exprf = '100.*sin(0.01* i)'; -my $exprp = '100.*sin(0.01*$i)'; -#$exprf = 'i'; -#$exprp = '$i'; - -# Two dimensional functions -my $expr2f = '100.*sin(0.01* i)*cos(0.01* j)'; -# no output autocreation means have this mess... -my $expr2p = '(outer(sin(0.01*$i),cos(0.01*$j),$c=null),$c*100.)'; - -# need to define/declare variables used in the expressions above -my $j = sequence($ndata)+1; -my $i = $j; -my $c; -# 1 dimensional -- - -# -# f77, implied & explicit swapping for 4 byte types, with 2 separate -# writes; and header array as well as header file -# -foreach my $pdltype ('float', 'long') { - print STDERR "Type $pdltype swapped\n" if $Verbose; - my $f77type = $types{$pdltype}; - my $val = $exprf; - $val = "char(int($val))" if $pdltype eq 'byte'; - - my $code = <<"EOT"; - -c Program to test i/o of F77 unformatted files - program rawtest - implicit none - integer i - $f77type a($ndata) - do i = 1, $ndata - a(i) = $val - enddo - - open(8,file= - \$'$data' - \$,status='new',form='unformatted') - i = $ndata - write (8) i - write (8) a - close(8) - end - -EOT - - createData $head, codefold($code); - byte4swap($data); - open(FILE, "> $hdr"); - print FILE <<"EOT"; -# FlexRaw file header -f77 -long 1 1 -# Data -$pdltype 1 $ndata -EOT - close(FILE); - - my @a = readflex($data); - # print "@a\n"; - ok my $ok = ($a[0]->at(0) == $ndata); - my $res = eval "$pdltype $exprp"; - is_pdl $res, $a[1], "readflex $pdltype w hdr file"; - - open(FILE,">$hdr"); - print FILE <<"EOT"; -# FlexRaw file header -swap -f77 -# now for data specifiers -long 1 1 -# Data -$pdltype 1 $ndata -EOT - close(FILE); - - @a = readflex($data); - #print "@a\n"; - - unlink $hdr; - - ok $ok = ($a[0]->at(0) == $ndata); - $res = eval "$pdltype $exprp"; - is_pdl $res, $a[1], "readflex $pdltype w hdr file (explicit swap)"; - -# Now try header array - $ok = 1; - my $header = [ {Type => 'f77'}, - {Type => 'long', NDims => 1, Dims => [ 1 ] }, - {Type => $pdltype, NDims => 1, Dims => [ $ndata ] } ]; - @a = readflex($data,$header); - unlink $data; - ok $ok = ($a[0]->at(0) == $ndata); - $res = eval "$pdltype $exprp"; - is_pdl $res, $a[1], "readflex $pdltype w hdr array"; - -} # foreach: $pdltype == 'float', 'double' - -# 1d, all types, normal way round, f77 specifier -foreach my $pdltype (sort keys %types) { - print STDERR "Type $pdltype\n" if $Verbose; - my $f77type = $types{$pdltype}; - my $val = $exprf; - $val = "char(int($val))" if $pdltype eq 'byte'; - - my $code = <<"EOT"; - -c Program to test i/o of F77 unformatted files - program rawtest - implicit none - integer i - $f77type a($ndata) - do i = 1, $ndata - a(i) = $val - enddo - open(8,file= - \$'$data' - \$,status='new',form='unformatted') - i = $ndata - write (8) i,a - close(8) - end - -EOT - - createData $head, codefold($code); - - open(FILE, ">$hdr" ); - print FILE <<"EOT"; -# FlexRaw file header -f77 -long 1 1 -# Data -$pdltype 1 $ndata -EOT - close(FILE); - - my @a = readflex($data); - # print "@a\n"; - unlink $data, $hdr; - - ok my $ok = ($a[0]->at(0) == $ndata); - my $res = eval "$pdltype $exprp"; - is_pdl $res, $a[1], "f77 1D $pdltype data"; - # print $a[1]->getndims()," [",$a[1]->dims,"]\n"; - -} # foreach: $pdltype ( keys %types ) - -# 1 dimensional, no f77 specifier (format words explicitly ignored) -foreach my $pdltype (sort keys %types) { - print STDERR "Type $pdltype\n" if $Verbose; - my $f77type = $types{$pdltype}; - my $val = $exprf; - $val = "char(int($val))" if $pdltype eq 'byte'; - - my $code = <<"EOT"; - -c Program to test i/o of F77 unformatted files - program rawtest - implicit none - integer i - $f77type a($ndata) - do i = 1, $ndata - a(i) = $val - enddo - open(8,file= - \$'$data' - \$,status='new',form='unformatted') - i = $ndata - write (8) i,a - close(8) - end - -EOT - - createData $head, codefold($code); - - open(FILE,">$hdr"); - print FILE <<"EOT"; -# FlexRaw header file -byte 1 4 -long 1 # Test comments -1 Tricky comment -# Data -$pdltype 1 $ndata -byte 1 4 -# and hanging EOF - - -EOT - close(FILE); - - my @a = readflex($data); - # print "@a\n"; - unlink $data, $hdr; - - ok my $ok = ($a[1]->at(0) == $ndata); - my $res = eval "$pdltype $exprp"; - is_pdl $res,$a[2], "no f77, 1D $pdltype data"; - # print $a[2]->getndims()," [",$a[2]->dims,"]\n"; -} - -# 2 dimensional -foreach my $pdltype (sort keys %types) { - print STDERR "Type $pdltype\n" if $Verbose; - my $f77type = $types{$pdltype}; - my $val = $expr2f; - $val = "char(int($val))" if $pdltype eq 'byte'; - - my $code = <<"EOT"; - -c Program to test i/o of F77 unformatted files - program rawtest - implicit none - integer i, j - $f77type a($ndata, $ndata) - do i = 1, $ndata - do j = 1, $ndata - a(i,j) = $val - enddo - enddo - open(8,file= - \$'$data' - \$,status='new',form='unformatted') - i = $ndata - write (8) i,a - close(8) - end - -EOT - - createData $head, codefold($code); - - open(FILE,">$hdr"); - print FILE <<"EOT"; -# FlexRaw file header -f77 -long 1 1 -# Data -$pdltype 2 $ndata $ndata -EOT - close(FILE); - my @a = readflex($data); -# if ($pdltype eq 'byte') { -# print "$pdltype @a\n"; -# system("ls -l $data"); -# } - unlink $data, $hdr; - - ok my $ok = ($a[0]->at(0) == $ndata); - my $res = eval "$pdltype $expr2p"; - is_pdl $res, $a[1], "f77 format 2D $pdltype data"; - # print $a[1]->getndims()," [",$a[1]->dims,"]\n"; -} - -print STDERR "Combined types case\n" if $Verbose; - -my $code = <<"EOT"; - -c Program to test i/o of F77 unformatted files - program rawtest - implicit none - character a - integer*2 i - integer*4 l - real*4 f - real*8 d - d = 4*atan(1.) - f = d - l = 10**d - i = l - a = ' ' - open(8,file= - \$'$data' - \$,status='new',form='unformatted') -c Choose bad boundaries... - write (8) a,i,l,f,d - close(8) - end - -EOT - - createData $head, codefold($code); - -open(FILE,">$hdr"); -print FILE <<"EOT"; -# FlexRaw file header -byte 1 4 -byte 0 -short 0 -long 0 -float 0 -double 0 -byte 1 4 -EOT -close(FILE); - -my @a = readflex($data); -#print "@a\n"; -shift @a; - -my $d = double pdl (4*atan2(1,1)); -my $f = float ($d); -my $l = long (10**$f); -$i = short ($l); -my $x = byte (32); -my @req = ($x,$i,$l,$f,$d); -foreach (@req) { - my $h = shift @a; - is_pdl $h, $_, "readflex combined types"; -} - -SKIP: { - my $compress = File::Which::which('compress') ? 'compress' : 'gzip'; # some linuxes don't have compress - $compress = 'gzip' if $^O eq 'cygwin'; # fix bogus compress script prob - - if ( $^O eq 'MSWin32' ) { # fix for ASPerl + MinGW, needs to be more general - skip "No compress or gzip command on MSWin32", 1 unless File::Which::which($compress) and $^O; - } - -# Try compressed data - 0 == system "$compress -c $data > ${data}.Z" or diag "system $compress -c $data >${data}.Z failed: $?"; - unlink( $data ); - @a = readflex($data); - ok $#a==6; - @a = readflex("${data}.Z"); - ok $#a==6; - my $NULL = File::Spec->devnull(); - 0 == system "gunzip -q ${data}.Z >$NULL 2>&1" or diag "system gunzip -q ${data}.Z failed: $?"; - 0 == system "gzip -q $data >$NULL 2>&1" or diag "system gzip -q $data failed: $?"; - @a = readflex($data); - ok $#a==6; - @a = readflex("${data}.gz"); - ok $#a==6; - shift @a; - unlink "${data}.gz", $hdr; - $d = double pdl (4*atan2(1,1)); - $f = float ($d); - $l = long (10**$f); - $i = short ($l); - $x = byte (32); - @req = ($x,$i,$l,$f,$d); - foreach (@req) { - my $h = shift @a; - is_pdl $h, $_, "readflex compressed data"; - } -} - -# Try writing data -my $flexhdr = writeflex($data,@req); -writeflexhdr($data,$flexhdr) unless $PDL::IO::FlexRaw::writeflexhdr; -@a = readflex($data); -unlink $hdr; -foreach (@req) { - is_pdl shift @a, $_, "writeflex combined data types, hdr file"; -} -@a = readflex($data, $flexhdr); -foreach (@req) { - is_pdl shift @a, $_, "writeflex combined data types, readflex hdr array"; -} -unlink $data; - -$#a = -1; -foreach (@req) { - push @a,$_->dummy(0,10); -} -$flexhdr = writeflex($data,@a); -$flexhdr = [ {Type => 'byte', NDims => 1, Dims => 10}, - {Type => 'short', NDims => 1, Dims => 10}, - {Type => 'long', NDims => 1, Dims => 10}, - {Type => 'float', NDims => 1, Dims => 10}, - {Type => 'double', NDims => 1, Dims => 10} ]; -@a = readflex($data, $flexhdr); -unlink $data; -foreach (@req) { - is_pdl slice(shift @a,"(0)"), $_, "writeflex combined types[10], readflex explicit hdr array"; -} - -# Writing multidimensional data -map {$_ = $_->dummy(0,10)} @req; -$flexhdr = writeflex($data,@req); -writeflexhdr($data,$flexhdr) unless $PDL::IO::FlexRaw::writeflexhdr; -@a = readflex($data); -unlink $data; -unlink $hdr; -foreach (@req) { - is_pdl shift @a, $_, "multidimensional data"; -} - -# Use readflex with an open file handle -@req = (byte(1..3), - long(5..10), - float(10..15)->reshape(3,2)/100, - double(0..99)/1e8); -$flexhdr = writeflex($data, @req); - -open(IN, $data); -@a = readflex(\*IN, $flexhdr); -foreach (@req) { - is_pdl shift @a, $_, "readflex with file handle"; -} -close(IN); -unlink $data; - -# use writeflex with an open file handle -open(OUT, ">$data"); -$flexhdr = writeflex(\*OUT, @req); -close(OUT); -@a = readflex($data, $flexhdr); -foreach (@req) { - is_pdl shift @a, $_, "writeflex with file handle"; -} -unlink $data; - -done_testing;