home *** CD-ROM | disk | FTP | other *** search
/ PC Professionell 2004 December / PCpro_2004_12.ISO / files / webserver / tsw / TSW_3.4.0.exe / Apache2 / perl / Matrix.pm < prev    next >
Encoding:
Perl POD Document  |  2002-09-09  |  9.2 KB  |  445 lines
  1. package Math::Cephes::Matrix;
  2. use strict;
  3. use vars qw(@EXPORT_OK $VERSION);
  4.  
  5. require Exporter;
  6. *import = \&Exporter::import;
  7. @EXPORT_OK = qw(mat);
  8.  
  9. $VERSION = '0.36';
  10.  
  11. require Math::Cephes;
  12.  
  13. sub new {
  14.   my ($caller, $arr) = @_;
  15.   my $refer = ref($caller);
  16.   my $class = $refer || $caller;
  17.   die "Must supply data for the matrix" 
  18.     unless ($refer or $arr);
  19.   unless ($refer) {
  20.     die "Please supply an array of arrays for the matrix data"
  21.       unless (ref($arr) eq 'ARRAY' and ref($arr->[0]) eq 'ARRAY');
  22.     my $n = scalar @$arr;
  23.     my $m = scalar @{$arr->[0]};
  24.     die "Matrices must be square" unless $m == $n;
  25.   }
  26.   my ($coef, $n);
  27.   if ($refer) {
  28.     $n = $caller->{n};
  29.     my $cdata = $caller->{coef};
  30.     foreach (@$cdata) {
  31.       push @$coef, [ @$_];
  32.     }
  33.   }
  34.   else {
  35.     ($coef, $n) = ($arr, scalar @$arr);
  36.   }
  37.   bless { coef => $coef,
  38.       n => $n,
  39.     }, $class;
  40. }
  41.  
  42. sub mat {
  43.   return Math::Cephes::Matrix->new(shift);
  44. }
  45.  
  46. sub mat_to_vec {
  47.   my $self = shift;
  48.   my ($M, $n) = ($self->{coef}, $self->{n});
  49.   my $A = [];
  50.   for (my $i=0; $i<$n; $i++) {
  51.     for (my $j=0; $j<$n; $j++) {
  52.       my $index = $i*$n+$j;
  53.       $A->[$index] = $M->[$i]->[$j];
  54.     }
  55.   }
  56.   return $A;
  57. }
  58.  
  59. sub vec_to_mat {
  60.   my ($self, $X) = @_;
  61.   my $n = $self->{n};
  62.   my $I = [];
  63.   for (my $i=0; $i<$n; $i++) {
  64.     for (my $j=0; $j<$n; $j++) {
  65.       my $index = $i*$n+$j;
  66.       $I->[$i]->[$j] = $X->[$index];
  67.     }
  68.   }
  69.   return $I;
  70. }
  71.  
  72. sub check {
  73.   my ($self, $B) = @_;
  74.   my $na = $self->{n};
  75.   my $ref = ref($B);
  76.   if ($ref eq 'Math::Cephes::Matrix') {
  77.     die "Matrices must be of the same size"
  78.       unless $B->{n} == $na;
  79.     return $B->coef;
  80.   }
  81.   elsif ($ref eq 'ARRAY') {
  82.     my $nb = scalar @$B;
  83.     my $ref0 = ref($B->[0]);
  84.     if ($ref0 eq 'ARRAY') {
  85.       my $m = scalar @{$B->[0]};
  86.       die "Can only use square matrices" unless $m == $nb;
  87.       die "Can only use matrices of the same size"
  88.     unless $na == $nb;
  89.       return $B;
  90.     }
  91.     elsif (not $ref0) {
  92.       die "Can only use vectors of the same size" unless $nb == $na;
  93.       return $B;
  94.     }
  95.     else {
  96.       die "Unknown reference '$ref0' for data";
  97.     }
  98.   }
  99.   else {
  100.     die "Unknown reference '$ref' for data";    
  101.   }
  102. }
  103.  
  104. sub coef {
  105.   return $_[0]->{coef};
  106. }
  107.  
  108. sub clr {
  109.   my $self = shift;
  110.   my $n = $self->{n};
  111.   my $B = [];
  112.   for (my $i=0; $i<$n; $i++) {
  113.     for (my $j=0; $j<$n; $j++) {
  114.       $B->[$i]->[$j] = 0;
  115.     }
  116.   }
  117.   $self->{coef} = $B;
  118. }
  119.  
  120. sub simq {
  121.   my ($self, $B) = @_;
  122.   $B = $self->check($B);
  123.   my ($M, $n) = ($self->{coef}, $self->{n});
  124.   die "Must supply an array reference for B" unless ref($B) eq 'ARRAY';
  125.   my $A = $self->mat_to_vec();
  126.   my $X = [split //, 0 x $n];
  127.   my $IPS = [split //, 0 x $n];
  128.   my $flag = 0;
  129.   my $ret = Math::Cephes::simq($A, $B, $X, $n, $flag, $IPS);
  130.   return $ret ? undef : $X;
  131. }
  132.  
  133.  
  134. sub inv {
  135.   my $self = shift;
  136.   my ($M, $n) = ($self->{coef}, $self->{n});
  137.   my $A = $self->mat_to_vec();
  138.   my $X = [split //, 0 x ($n*$n)];
  139.   my $B = [split //, 0 x $n];
  140.   my $IPS = [split //, 0 x $n];
  141.   my $flag = 0;
  142.   my $ret = Math::Cephes::minv($A, $X, $n, $B, $IPS);
  143.   return undef if $ret;
  144.   my $I = $self->vec_to_mat($X);
  145.   return Math::Cephes::Matrix->new($I);
  146. }
  147.  
  148. sub transp {
  149.   my $self = shift;
  150.   my ($M, $n) = ($self->{coef}, $self->{n});
  151.   my $A = $self->mat_to_vec();
  152.   my $T = [split //, 0 x ($n*$n)];
  153.   Math::Cephes::mtransp($n, $A, $T);
  154.   my $R = $self->vec_to_mat($T);
  155.   return Math::Cephes::Matrix->new($R);
  156. }
  157.  
  158. sub add {
  159.   my ($self, $B) = @_;
  160.   $B = $self->check($B);
  161.   my ($A, $n) = ($self->{coef}, $self->{n});
  162.   my $C = [];
  163.   for (my $i=0; $i<$n; $i++) {
  164.     for (my $j=0; $j<$n; $j++) {
  165.       $C->[$i]->[$j] = $A->[$i]->[$j] + $B->[$i]->[$j];
  166.     }
  167.   }
  168.   return Math::Cephes::Matrix->new($C);
  169. }
  170.  
  171. sub sub {
  172.   my ($self, $B) = @_;
  173.   $B = $self->check($B);
  174.   my ($A, $n) = ($self->{coef}, $self->{n});
  175.   my $C = [];
  176.   for (my $i=0; $i<$n; $i++) {
  177.     for (my $j=0; $j<$n; $j++) {
  178.       $C->[$i]->[$j] = $A->[$i]->[$j] - $B->[$i]->[$j];
  179.     }
  180.   }
  181.   return Math::Cephes::Matrix->new($C);
  182. }
  183.  
  184. sub mul {
  185.   my ($self, $B) = @_;
  186.   $B = $self->check($B);
  187.   my ($A, $n) = ($self->{coef}, $self->{n});
  188.   my $C = [];
  189.   if (ref($B->[0]) eq 'ARRAY') {
  190.     for (my $i=0; $i<$n; $i++) {
  191.       for (my $j=0; $j<$n; $j++) {
  192.     for (my $m=0; $m<$n; $m++) {
  193.       $C->[$i]->[$j] += $A->[$i]->[$m] * $B->[$m]->[$j];
  194.     }
  195.       }
  196.     }
  197.     return Math::Cephes::Matrix->new($C);
  198.   }
  199.   else {
  200.     for (my $i=0; $i<$n; $i++) {
  201.       for (my $m=0; $m<$n; $m++) {
  202.     $C->[$i] += $A->[$i]->[$m] * $B->[$m];
  203.       }
  204.     }
  205.     return $C;
  206.   }
  207. }
  208.  
  209. sub div {
  210.   my ($self, $B) = @_;
  211.   $B = $self->check($B);
  212.   my $C = Math::Cephes::Matrix->new($B)->inv();
  213.   my $D = $self->mul($C);
  214.   return $D;
  215. }
  216.  
  217. sub eigens {
  218.   my $self = shift;
  219.   my ($M, $n) = ($self->{coef}, $self->{n});
  220.   my $A = [];
  221.   for (my $i=0; $i<$n; $i++) {
  222.     for (my $j=0; $j<$n; $j++) {
  223.       my $index = ($i*$i+$i)/2 + $j;
  224.       $A->[$index] = $M->[$i]->[$j];
  225.     }
  226.   }
  227.   my $EV1 = [split //, 0 x ($n*$n)];
  228.   my $E = [split //, 0 x $n];
  229.   my $IPS = [split //, 0 x $n];
  230.   Math::Cephes::eigens($A, $EV1, $E, $n);
  231.   my $EV = $self->vec_to_mat($EV1);
  232.   return ($E, Math::Cephes::Matrix->new($EV));
  233. }
  234.  
  235. 1;
  236.  
  237. __END__
  238.  
  239. =head1 NAME
  240.  
  241. Math::Cephes::Matrix - Perl interface to the cephes matrix routines
  242.  
  243. =head1 SYNOPSIS
  244.  
  245.   use Math::Cephes::Matrix qw(mat);
  246.   # 'mat' is a shortcut for Math::Cephes::Matrix->new
  247.   my $M = mat([ [1, 2, -1], [2, -3, 1], [1, 0, 3]]);
  248.   my $C = mat([ [1, 2, 4], [2, 9, 2], [6, 2, 7]]);
  249.   my $D = $M->add($C);          # D = M + C
  250.   my $Dc = $D->coef;
  251.   for (my $i=0; $i<3; $i++) {
  252.     print "row $i:\n";
  253.     for (my $j=0; $j<3; $j++) {
  254.         print "\tcolumn $j: $Dc->[$i]->[$j]\n";
  255.     }
  256.   }
  257.  
  258. =head1 DESCRIPTION
  259.  
  260. This module is a layer on top of the basic routines in the
  261. cephes math library for operations on square matrices. In 
  262. the following, a Math::Cephes::Matrix object is created as
  263.  
  264.   my $M = Math::Cephes::Matrix->new($arr_ref);
  265.  
  266. where C<$arr_ref> is a reference to an array of arrays, as 
  267. in the following example:
  268.  
  269.   $arr_ref = [ [1, 2, -1], [2, -3, 1], [1, 0, 3] ]
  270.  
  271. which represents
  272.  
  273.      / 1   2  -1  \
  274.      | 2  -3   1  |
  275.      \ 1   0   3  /
  276.  
  277. A copy of a I<Math::Cephes::Matrix> object may be done as
  278.  
  279.   my $M_copy = $M->new();
  280.  
  281. =head2 Methods
  282.  
  283. =over 4
  284.  
  285. =item I<coef>: get coefficients of the matrix
  286.  
  287.  SYNOPSIS:
  288.  
  289.  my $c = $M->coef;
  290.  
  291.  DESCRIPTION:
  292.  
  293. This returns an reference to an array of arrays
  294. containing the coefficients of the matrix. 
  295.  
  296. =item I<clr>: set all coefficients equal to zero
  297.  
  298.  SYNOPSIS:
  299.  
  300.  $M->clr();
  301.  
  302.  DESCRIPTION:
  303.  
  304. This sets all the coefficients of the matrixidentically to 0.
  305.  
  306. =item I<add>: add two matrices
  307.  
  308.  SYNOPSIS:
  309.  
  310.  $P = $M->add($N);
  311.  
  312.  DESCRIPTION:
  313.  
  314. This sets $P equal to $M + $N.
  315.  
  316. =item I<sub>: subtract two matrices
  317.  
  318.  SYNOPSIS:
  319.  
  320.  $P = $M->sub($N);
  321.  
  322.  DESCRIPTION:
  323.  
  324. This sets $P equal to $M - $N.
  325.  
  326. =item I<mul>: multiply two matrices or a matrix and a vector
  327.  
  328.  SYNOPSIS:
  329.  
  330.  $P = $M->mul($N);
  331.  
  332.  DESCRIPTION:
  333.  
  334. This sets $P equal to $M * $N. This method can handle
  335. matrix multiplication, when $N is a matrix, as well
  336. as matrix-vector multiplication, where $N is an
  337. array reference representing a column vector.
  338.  
  339. =item I<div>: divide two matrices
  340.  
  341.  SYNOPSIS:
  342.  
  343.  $P = $M->div($N);
  344.  
  345.  DESCRIPTION:
  346.  
  347. This sets $P equal to $M * ($N)^(-1).
  348.  
  349. =item I<inv>: invert a matrix
  350.  
  351.  SYNOPSIS:
  352.  
  353.  $I = $M->inv();
  354.  
  355.  DESCRIPTION:
  356.  
  357. This sets $I equal to ($M)^(-1).
  358.  
  359. =item I<transp>: transpose a matrix
  360.  
  361.  SYNOPSIS:
  362.  
  363.  $T = $M->transp();
  364.  
  365.  DESCRIPTION:
  366.  
  367. This sets $T equal to the transpose of $M.
  368.  
  369. =item I<simq>: solve simultaneous equations
  370.  
  371.  SYNOPSIS:
  372.  
  373.  my $M = Math::Cephes::Matrix->new([ [1, 2, -1], [2, -3, 1], [1, 0, 3]]);
  374.  my $B = [2, -1, 10];
  375.  my $X = $M->simq($B);
  376.  for (my $i=0; $i<3; $i++) {
  377.     print "X[$i] is $X->[$i]\n";
  378.   }
  379.  
  380. where $M is a I<Math::Cephes::Matrix> object, $B
  381. is an input array reference, and $X is an output
  382. array reference.
  383.  
  384.  DESCRIPTION:
  385.  
  386. A set of N simultaneous equations may be represented
  387. in matrix form as
  388.  
  389.   M X = B
  390.  
  391. where M is an N x N square matrix and X and B are column
  392. vectors of length N. 
  393.  
  394. =item I<eigens>: eigenvalues and eigenvectors of a real symmetric matrix
  395.  
  396.  SYNOPSIS:
  397.  
  398.  my $S = Math::Cephes::Matrix->new([ [1, 2, 3], [2, 2, 3], [3, 3, 4]]);
  399.  my ($E, $EV1) = $S->eigens();
  400.  my $EV = $EV1->coef;
  401.  for (my $i=0; $i<3; $i++) {
  402.    print "For i=$i, with eigenvalue $E->[$i]\n";
  403.    my $v = [];
  404.    for (my $j=0; $j<3; $j++) {
  405.      $v->[$j] = $EV->[$i]->[$j];
  406.    }
  407.    print "The eigenvector is @$v\n";
  408.  }
  409.  
  410. where $M is a I<Math::Cephes::Matrix> object representing
  411. a real symmetric matrix. $E is an array reference containing
  412. the eigenvalues of $M, and $EV is a I<Math::Cephes::Matrix> object 
  413. representing the eigenvalues, the I<ith> row corresponding
  414. to the I<ith> eigenvalue.
  415.  
  416.  DESCRIPTION:
  417.  
  418. If M is an N x N real symmetric matrix, and X is an N component
  419. column vector, the eigenvalue problem
  420.  
  421.   M X = lambda X
  422.  
  423. will in general have N solutions, with X the eigenvectors
  424. and lambda the eigenvalues.
  425.  
  426. =back
  427.  
  428. =head1 BUGS
  429.  
  430. Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
  431.  
  432. =head1 COPYRIGHT
  433.  
  434. The C code for the Cephes Math Library is
  435. Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier, 
  436. and is available at http://www.netlib.org/cephes/.
  437. Direct inquiries to 30 Frost Street, Cambridge, MA 02140.
  438.  
  439. The perl interface is copyright 2000, 2002 by Randy Kobes.
  440. This library is free software; you can redistribute it and/or
  441. modify it under the same terms as Perl itself.
  442.  
  443. =cut
  444.  
  445.