package Math::Cephes::Matrix;
use strict;
use warnings;
use vars qw(@EXPORT_OK $VERSION);

require Exporter;
*import = \&Exporter::import;
@EXPORT_OK = qw(mat);

$VERSION = '0.5308';

require Math::Cephes;

sub new {
  my ($caller, $arr) = @_;
  my $refer = ref($caller);
  my $class = $refer || $caller;
  die "Must supply data for the matrix"
    unless ($refer or $arr);
  unless ($refer) {
    die "Please supply an array of arrays for the matrix data"
      unless (ref($arr) eq 'ARRAY' and ref($arr->[0]) eq 'ARRAY');
    my $n = scalar @$arr;
    my $m = scalar @{$arr->[0]};
    die "Matrices must be square" unless $m == $n;
  }
  my ($coef, $n);
  if ($refer) {
    $n = $caller->{n};
    my $cdata = $caller->{coef};
    foreach (@$cdata) {
      push @$coef, [ @$_];
    }
  }
  else {
    ($coef, $n) = ($arr, scalar @$arr);
  }
  bless { coef => $coef,
	  n => $n,
	}, $class;
}

sub mat {
  return Math::Cephes::Matrix->new(shift);
}

sub mat_to_vec {
  my $self = shift;
  my ($M, $n) = ($self->{coef}, $self->{n});
  my $A = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      my $index = $i*$n+$j;
      $A->[$index] = $M->[$i]->[$j];
    }
  }
  return $A;
}

sub vec_to_mat {
  my ($self, $X) = @_;
  my $n = $self->{n};
  my $I = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      my $index = $i*$n+$j;
      $I->[$i]->[$j] = $X->[$index];
    }
  }
  return $I;
}

sub check {
  my ($self, $B) = @_;
  my $na = $self->{n};
  my $ref = ref($B);
  if ($ref eq 'Math::Cephes::Matrix') {
    die "Matrices must be of the same size"
      unless $B->{n} == $na;
    return $B->coef;
  }
  elsif ($ref eq 'ARRAY') {
    my $nb = scalar @$B;
    my $ref0 = ref($B->[0]);
    if ($ref0 eq 'ARRAY') {
      my $m = scalar @{$B->[0]};
      die "Can only use square matrices" unless $m == $nb;
      die "Can only use matrices of the same size"
	unless $na == $nb;
      return $B;
    }
    elsif (not $ref0) {
      die "Can only use vectors of the same size" unless $nb == $na;
      return $B;
    }
    else {
      die "Unknown reference '$ref0' for data";
    }
  }
  else {
    die "Unknown reference '$ref' for data";
  }
}

sub coef {
  return $_[0]->{coef};
}

sub clr {
  my $self = shift;
  my $what = shift || 0;
  my $n = $self->{n};
  my $B = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      $B->[$i]->[$j] = $what;
    }
  }
  $self->{coef} = $B;
}

sub simq {
  my ($self, $B) = @_;
  $B = $self->check($B);
  my ($M, $n) = ($self->{coef}, $self->{n});
  die "Must supply an array reference for B" unless ref($B) eq 'ARRAY';
  my $A = $self->mat_to_vec();
  my $X = [split //, 0 x $n];
  my $IPS = [split //, 0 x $n];
  my $flag = 0;
  my $ret = Math::Cephes::simq($A, $B, $X, $n, $flag, $IPS);
  return $ret ? undef : $X;
}


sub inv {
  my $self = shift;
  my ($M, $n) = ($self->{coef}, $self->{n});
  my $A = $self->mat_to_vec();
  my $X = [split //, 0 x ($n*$n)];
  my $B = [split //, 0 x $n];
  my $IPS = [split //, 0 x $n];
  my $flag = 0;
  my $ret = Math::Cephes::minv($A, $X, $n, $B, $IPS);
  return undef if $ret;
  my $I = $self->vec_to_mat($X);
  return Math::Cephes::Matrix->new($I);
}

sub transp {
  my $self = shift;
  my ($M, $n) = ($self->{coef}, $self->{n});
  my $A = $self->mat_to_vec();
  my $T = [split //, 0 x ($n*$n)];
  Math::Cephes::mtransp($n, $A, $T);
  my $R = $self->vec_to_mat($T);
  return Math::Cephes::Matrix->new($R);
}

sub add {
  my ($self, $B) = @_;
  $B = $self->check($B);
  my ($A, $n) = ($self->{coef}, $self->{n});
  my $C = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      $C->[$i]->[$j] = $A->[$i]->[$j] + $B->[$i]->[$j];
    }
  }
  return Math::Cephes::Matrix->new($C);
}

sub sub {
  my ($self, $B) = @_;
  $B = $self->check($B);
  my ($A, $n) = ($self->{coef}, $self->{n});
  my $C = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      $C->[$i]->[$j] = $A->[$i]->[$j] - $B->[$i]->[$j];
    }
  }
  return Math::Cephes::Matrix->new($C);
}

sub mul {
  my ($self, $B) = @_;
  $B = $self->check($B);
  my ($A, $n) = ($self->{coef}, $self->{n});
  my $C = [];
  if (ref($B->[0]) eq 'ARRAY') {
    for (my $i=0; $i<$n; $i++) {
      for (my $j=0; $j<$n; $j++) {
	for (my $m=0; $m<$n; $m++) {
	  $C->[$i]->[$j] += $A->[$i]->[$m] * $B->[$m]->[$j];
	}
      }
    }
    return Math::Cephes::Matrix->new($C);
  }
  else {
    for (my $i=0; $i<$n; $i++) {
      for (my $m=0; $m<$n; $m++) {
	$C->[$i] += $A->[$i]->[$m] * $B->[$m];
      }
    }
    return $C;
  }
}

sub div {
  my ($self, $B) = @_;
  $B = $self->check($B);
  my $C = Math::Cephes::Matrix->new($B)->inv();
  my $D = $self->mul($C);
  return $D;
}

sub eigens {
  my $self = shift;
  my ($M, $n) = ($self->{coef}, $self->{n});
  my $A = [];
  for (my $i=0; $i<$n; $i++) {
    for (my $j=0; $j<$n; $j++) {
      my $index = ($i*$i+$i)/2 + $j;
      $A->[$index] = $M->[$i]->[$j];
    }
  }
  my $EV1 = [split //, 0 x ($n*$n)];
  my $E = [split //, 0 x $n];
  my $IPS = [split //, 0 x $n];
  Math::Cephes::eigens($A, $EV1, $E, $n);
  my $EV = $self->vec_to_mat($EV1);
  return ($E, Math::Cephes::Matrix->new($EV));
}

1;

__END__

=head1 NAME

Math::Cephes::Matrix - Perl interface to the cephes matrix routines

=head1 SYNOPSIS

  use Math::Cephes::Matrix qw(mat);
  # 'mat' is a shortcut for Math::Cephes::Matrix->new
  my $M = mat([ [1, 2, -1], [2, -3, 1], [1, 0, 3]]);
  my $C = mat([ [1, 2, 4], [2, 9, 2], [6, 2, 7]]);
  my $D = $M->add($C);          # D = M + C
  my $Dc = $D->coef;
  for (my $i=0; $i<3; $i++) {
    print "row $i:\n";
    for (my $j=0; $j<3; $j++) {
        print "\tcolumn $j: $Dc->[$i]->[$j]\n";
    }
  }

=head1 DESCRIPTION

This module is a layer on top of the basic routines in the
cephes math library for operations on square matrices. In
the following, a Math::Cephes::Matrix object is created as

  my $M = Math::Cephes::Matrix->new($arr_ref);

where C<$arr_ref> is a reference to an array of arrays, as
in the following example:

  $arr_ref = [ [1, 2, -1], [2, -3, 1], [1, 0, 3] ]

which represents

     / 1   2  -1  \
     | 2  -3   1  |
     \ 1   0   3  /

A copy of a I<Math::Cephes::Matrix> object may be done as

  my $M_copy = $M->new();

=head2 Methods

=over 4

=item I<coef>: get coefficients of the matrix

 SYNOPSIS:

 my $c = $M->coef;

 DESCRIPTION:

This returns an reference to an array of arrays
containing the coefficients of the matrix.

=item I<clr>: set all coefficients equal to a value.

 SYNOPSIS:

 $M->clr($n);

 DESCRIPTION:

This sets all the coefficients of the matrix identically to I<$n>.
If I<$n> is not given, a default of 0 is used.

=item I<add>: add two matrices

 SYNOPSIS:

 $P = $M->add($N);

 DESCRIPTION:

This sets $P equal to $M + $N.

=item I<sub>: subtract two matrices

 SYNOPSIS:

 $P = $M->sub($N);

 DESCRIPTION:

This sets $P equal to $M - $N.

=item I<mul>: multiply two matrices or a matrix and a vector

 SYNOPSIS:

 $P = $M->mul($N);

 DESCRIPTION:

This sets $P equal to $M * $N. This method can handle
matrix multiplication, when $N is a matrix, as well
as matrix-vector multiplication, where $N is an
array reference representing a column vector.

=item I<div>: divide two matrices

 SYNOPSIS:

 $P = $M->div($N);

 DESCRIPTION:

This sets $P equal to $M * ($N)^(-1).

=item I<inv>: invert a matrix

 SYNOPSIS:

 $I = $M->inv();

 DESCRIPTION:

This sets $I equal to ($M)^(-1).

=item I<transp>: transpose a matrix

 SYNOPSIS:

 $T = $M->transp();

 DESCRIPTION:

This sets $T equal to the transpose of $M.

=item I<simq>: solve simultaneous equations

 SYNOPSIS:

 my $M = Math::Cephes::Matrix->new([ [1, 2, -1], [2, -3, 1], [1, 0, 3]]);
 my $B = [2, -1, 10];
 my $X = $M->simq($B);
 for (my $i=0; $i<3; $i++) {
    print "X[$i] is $X->[$i]\n";
  }

where $M is a I<Math::Cephes::Matrix> object, $B
is an input array reference, and $X is an output
array reference.

 DESCRIPTION:

A set of N simultaneous equations may be represented
in matrix form as

  M X = B

where M is an N x N square matrix and X and B are column
vectors of length N.

=item I<eigens>: eigenvalues and eigenvectors of a real symmetric matrix

 SYNOPSIS:

 my $S = Math::Cephes::Matrix->new([ [1, 2, 3], [2, 2, 3], [3, 3, 4]]);
 my ($E, $EV1) = $S->eigens();
 my $EV = $EV1->coef;
 for (my $i=0; $i<3; $i++) {
   print "For i=$i, with eigenvalue $E->[$i]\n";
   my $v = [];
   for (my $j=0; $j<3; $j++) {
     $v->[$j] = $EV->[$i]->[$j];
   }
   print "The eigenvector is @$v\n";
 }

where $M is a I<Math::Cephes::Matrix> object representing
a real symmetric matrix. $E is an array reference containing
the eigenvalues of $M, and $EV is a I<Math::Cephes::Matrix> object
representing the eigenvalues, the I<ith> row corresponding
to the I<ith> eigenvalue.

 DESCRIPTION:

If M is an N x N real symmetric matrix, and X is an N component
column vector, the eigenvalue problem

  M X = lambda X

will in general have N solutions, with X the eigenvectors
and lambda the eigenvalues.

=back

=head1 BUGS

Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>

=head1 COPYRIGHT

The C code for the Cephes Math Library is
Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier,
and is available at http://www.netlib.org/cephes/.
Direct inquiries to 30 Frost Street, Cambridge, MA 02140.

The perl interface is copyright 2000, 2002 by Randy Kobes.
This library is free software; you can redistribute it and/or
modify it under the same terms as Perl itself.

=cut