#!/usr/bin/perl
############################################################
# Module name: matrix;
#
# Usage/Description: See POD at the bottom
############################################################
package matrix;
use warnings;
use strict;

use overload(
	     '""'=>\&display,
	     '*'=>\&multiply,
	     '+'=>\&add,
	     '-'=>\&subtract) ;

############################################################
# method name: new
#
# Usage:
#
#   my $matrix=new matrix(rows=>$rows,columns=>$columns)
#
# Description: Creates a new matrix with the number
#   of rows and columns specified.
############################################################
sub new {
  my $class=shift;

  # place command line arguments into the class
  my $self=bless {@_},$class;

  # create the C data
  $self->{ _cdata } = matrix_allocate( $self->{rows},$self->{columns} );

  return $self;

}


############################################################
# method name: clone
#
# Usage:
# 
#   my $idential_matrix=$matrix->clone;
#
# Description: creates a clone of a matrix
#   and returns the new matrix. This is a deep copy.
#   no reference to the old one is kept in it.
############################################################
sub clone {
  my $matrix=shift;

  my $clone_matrix=new matrix(rows=>$matrix->rows,columns=>$matrix->columns);

  #loop through all rows and columns
  for(my $row_loop=0;$row_loop<$matrix->rows;$row_loop++) {
    for (my $column_loop=0;$column_loop<$matrix->columns;$column_loop++) {
      # set the clone value to the acutal value
      $clone_matrix->set($row_loop,$column_loop,
			 $matrix->get($row_loop,$column_loop) );
    }
  }
  return $clone_matrix;
}



############################################################
# method name: get
#
# Usage:
# 
#   my $value=$matrix->get($row,$column);
#
# Description: returns the value of (row,column).
############################################################
sub get {
  my ($self, $row,$column) = @_;
  return  matrix_get( $self->{_cdata}, $row,$column );
}


############################################################
# method name: set
#
# Usage:
#
#   $matrix->set($row,$column,$value);
#   $matrix->set($value); #sets all values
#   $matrix->set([[1,2],[2,2]]);
#
# Description: sets the value of (row,column) to value.
############################################################
sub set {
  my $self=shift;

  if (scalar @_ == 1) {
    my ($value)=@_;
    if (ref $value eq "array") { #If value passed is an array ref

      # then put data structure in array
      for (my $row_loop; #loop through the rows
	   $row_loop<@$value and $row_loop<$self->row;
	   $row_loop++) {	
	for (my $column_loop; # loop through the columns
	     $column_loop<(@{$value->[$row_loop]}) and $column_loop<$self->column;
	     $column_loop++) {
	  # put into the matrix coresponding data structure
	  matrix_set( $self->{_cdata},$row_loop,
		      $column_loop, $value->[$row_loop][$column_loop] );
	}
      }
      return;
    } else {
      # otherwise, set all values to the scalar
      matrix_set_all_values($self->{_cdata},$value);
      return;
    }
  }

  my ($row,$column,$value)=@_;
  return 0 if (not defined $row or
	       not defined $column or
	       not defined $value);

  matrix_set( $self->{_cdata}, $row, $column, $value );

}


############################################################
# method name: rows
#
# Usage:
#
#   my $rows=$matrix->rows;
#
# Description: returns the number of rows
############################################################
sub rows {
  my $self=shift;
  return  matrix_rows ( $self->{_cdata});
}

############################################################
# method name: columns
#
# Usage:
#
#   my $rows=$matrix->columns;
#
# Description: returns the number of columns.
############################################################
sub columns {
  my $self=shift;
  return  matrix_columns( $self->{_cdata});
}


############################################################
# method name: add
#
# Usage:
#
#   my $result=$matrix_1+$matrix_2;
#   my $result=$matrix_1->add($matrix2);
#
# Description: Adds 2 matricies together and returns the 
#   result.
############################################################
sub add {
  my ($matrix_1,$matrix_2)=@_;

  return undef if (($matrix_1->rows != $matrix_2->rows)||
		   ($matrix_1->columns != $matrix_2->columns));

  my $resultant_matrix=new matrix(rows=>$matrix_1->rows,
			       columns=>$matrix_1->columns);

  matrix_add($matrix_1->{_cdata},
	     $matrix_2->{_cdata},
	     $resultant_matrix->{_cdata});
  return $resultant_matrix;
}


############################################################
# method name: subtract
#
# Usage:
#
#   my $result=$matrix_1-$matrix_2;
#   my $result=$matrix_1->subtract($matrix2);
#
# Description: subtracts 2 matricies and returns the 
#   result.
############################################################
sub subtract {
  my ($matrix_1,$matrix_2)=@_;
  return undef if (($matrix_1->rows != $matrix_2->rows)||
		   ($matrix_1->columns != $matrix_2->columns));

  my $resultant_matrix=new matrix(rows=>$matrix_1->rows,
			       columns=>$matrix_1->columns);
  matrix_subtract($matrix_1->{_cdata},
		  $matrix_2->{_cdata},
		  $resultant_matrix->{_cdata});
  return $resultant_matrix;
}


############################################################
# method name: display
#
# Usage:
#
#   print $matrix;
#   print $matrix->display;
#
# Description: returns a string representation of a matrix.
############################################################
sub display {
  my $self=shift;
  my $return_string='';

  # calculate the largest string
  my $max_length=0;
  foreach my $row (0..$self->rows-1) {
    foreach my $column (0..$self->columns-1) {
      $max_length=length $self->get($row,$column)
	if length $self->get($row,$column)> $max_length;
    }
  }

  foreach my $row (0..$self->rows-1) {
    $return_string.="[ ";
    foreach my $column (0..$self->columns-1) {
      # format string so all items align.
      $return_string.=sprintf "%".$max_length."d ",$self->get($row,$column);
    }
    $return_string.="]\n";
  }
  return $return_string;
}


############################################################
# method name: multiply
#
# Usage:
#
#  my $result=$matrix_1 * $matrix_2.
#  my $result=$scalar*$matrix;
#  my $result=$matrix*$scalar;
#
#  # or alternativly
#
#  my $result=$matrix_2->multiply($matrix_2);
#  my $result=$matrix_1->multiplY($scalar);
#
# Description: Performs either matrix or scalar 
#   multiplication and returns the result.
############################################################
sub multiply {
  my ($LHS,$RHS)=@_;
  my $resultant_matrix;

  if (not ref $LHS and not ref $RHS) { #if neither LHS or RHS are matricies
    $resultant_matrix=$LHS*$RHS;

  } elsif (not ref $LHS) { #if the LHS is a scalar and RHS is a matrix
    $resultant_matrix=new matrix(rows=>$RHS->rows,
			  columns=>$RHS->columns);
    scalar_multiply($RHS->{_cdata},$LHS,$resultant_matrix->{_cdata});

  } elsif (not ref $RHS) { # if the RHS is a scalar and RHS is a matrix
    $resultant_matrix=new matrix(rows=>$LHS->rows,
			  columns=>$LHS->columns);
    scalar_multiply($LHS->{_cdata},$RHS,$resultant_matrix->{_cdata});

  } else { # if both LHS and RHS are matricies

    # cols on left must equal rows on right
    return undef if ($LHS->columns != $RHS->rows);

    $resultant_matrix=new matrix(rows=>$LHS->rows,columns=>$RHS->columns);
    #matrix_multiply($LHS->{_cdata},$RHS->{_cdata},$resultant_matrix->{_cdata});
    matrix_multiply($LHS->{_cdata},$RHS->{_cdata},$resultant_matrix->{_cdata});
  }
  return $resultant_matrix;
}




############################################################
# method name: lu_decompose
#
# Usage:
#
#   lu_decompose has 2 usages. 
#
#   my $lu=$matrix->lu_decompose;
#   # or
#   $matrix->lu_decompose($lower,$upper)
#
#   # in the first case, the flattend LU decomposition
#   # is returned. In the second, the 2 lu decomposed
#   # matricies are returned through $lower and $upper.
#   # in that case, $lower*$upper=$matrix
#
# Description: Performs the LU decomposition using 
#   crout's algorithm. This algorithm works by assigning
#   the diagonal of the lower matrix to 1. The rest
#   of the values are determined so that when the lower
#   and upper matrix are multiplied the input matrix
#   is recieved.
#
#   If this function is called with no parameters,
#   the lu decomposition, flatted into one NxN matrix,
#   is returned (what is missing is the overlap of the 1s
#   in the lower matrix).
#
#   If this funciton is called with 2 matricies,
#   the LU decomposition is reutrned.
############################################################
sub lu_decompose {
  my $self=shift;

  return undef if  ($self->rows != $self->columns);
  my ($lower,$upper)=(shift,shift);

  # if either of the passed argumetns are not matricies,
  # then, return a flattened decomposition matrix
  if (ref $lower ne __PACKAGE__ or 
      ref $upper ne __PACKAGE__) {
    my $lu_decompose=new matrix(rows=>$self->rows,
				columns=>$self->columns); 
    # make a copy of the matrix to do the decomposition on.
    $lu_decompose=$self->clone;
    matrix_lu_decompose($lu_decompose->{_cdata});
    return $lu_decompose;
  }
  # otherwsie, we are passed 2 matricies.
  # the decomposition should placed inside 
  # the 2 matricies.

  # if the lower matrix is not the right dimension
  if ($lower->rows != $self->rows or
      $lower->columns != $self->columns) {
    my $temp_matrix=new matrix(rows=>$self->rows,
			       columns=>$self->columns);
    $lower->{_cdata}=$temp_matrix->{_cdata};
    # we do this so the caller still pointers at lower
  }

  $upper->{_cdata}=$self->clone->{_cdata};
  $lower->set(0);
  matrix_lu_decompose($upper->{_cdata});

  # now, we must move the lower part of $upper into $lower
  for (my $column_loop=1;$column_loop<$upper->columns;$column_loop++) {
    for (my $row_loop=0;$row_loop<$column_loop;$row_loop++) {
      # move lower triangular part of upper triangular matrix
      # onto the $lower. Then erase it
      $lower->set($column_loop,$row_loop,
		  $upper->get($column_loop,$row_loop));
      $upper->set($column_loop,$row_loop,0);
    }
  }

  for (my $diagonal_loop=0;$diagonal_loop<$upper->columns;$diagonal_loop++) {
    $lower->set($diagonal_loop,$diagonal_loop,1)
  }

  # $lower and $upper are now LU decomposed $self
  return;
}



# sub print_multiply_tables
# simply used in debuggin to print out tables
# creatind doing the chain multiplication
sub print_multiply_tables {
  my ($costs)=@_;

  my $return_string;
  # calculate the largest string
  my $max_length=0;
  foreach my $row (0..scalar @{$costs}-1) {
    foreach my $column (0..scalar @{$costs->[0]}-1) {
      # so we can print all values
      $costs->[$column][$row]=0 if not defined $costs->[$column][$row];
      $max_length=length $costs->[$column][$row]
	if length $costs->[$column][$row] > $max_length;
    }
  }
  foreach my $row (0..scalar @{$costs}-1) {
    $return_string.="[ ";
    foreach my $column (0..scalar @{$costs->[0]}-1) {
      # format string so all items align.

      $return_string.=sprintf "%".$max_length."d ",$costs->[$column][$row];
    }
    $return_string.="]\n";
  }
  return $return_string;
}



sub multiply_chain_solve {
  my ($matricies,$best,$lower_index,$upper_index)=@_;
  
  # base case, one matrix, return it
  return "m_$lower_index" 
    if (scalar@{$matricies}==1 or $lower_index==$upper_index);
  
  # find split point in best array;
  my $split=$best->[$lower_index][$upper_index];
  
  # multiply left and right parts
  my $left_product=multiply_chain_solve([@{$matricies}[$lower_index..$split]],
					$best,
					$lower_index,
					$split);
  my $right_product=multiply_chain_solve([@{$matricies}[$split+1..$upper_index]],
					 $best,
					 $split+1,
					 $upper_index);
  # return their product multiplied
  return "( $left_product * $right_product )";
}


############################################################
# method name: multiply_chain
#
# Usage:
#
#  my $product=multiply_chain($matrix_1,$matrix_2,$matrix_3
#                             $matrix_4,$matrix_5,$matrix6);
#
# Description: does the product of all matricies inputted.#
#    return the product. Multiplies the matricies in 
#    a specail order to optomize the multiplication.
############################################################
sub multiply_chain {
  my @matricies=@_;
  # all input must be matricies
  foreach (@matricies){
    return (undef) if (ref $_ ne __PACKAGE__);
  }

  # ensure all matrix multiplications are valid
  for(my $loop=0;$loop<scalar @matricies-1;$loop++) {
    return undef if ($matricies[$loop]->columns !=
		     $matricies[$loop+1]->rows);
  }
  my ($costs,$dims,$best)=([[]],[[]],[[]]);

  # initalize the three auxiliary matriceis that we'll use to
  # store the costs, the dimensions and the best
  for (my $i=0;$i<scalar @matricies;$i++) {
    $costs->[$i][$i]=0;
    $dims->[$i][$i]=[$matricies[$i]->rows,$matricies[$i]->columns];
  }

  # i loop for pairs, triples, quadruples
  for (my $i=1;$i<scalar @matricies;$i++) {
    # j loops for each pair, triple, quad
    for (my $j=$i;$j<scalar @matricies;$j++) {
      for (my $k=$j-$i;$k<$j;$k++) {
	my ($row1,$col1) = @{$dims->[$j-$i][$k]};
	my ($row2,$col2) = @{$dims->[$k+1][$j]};
	# compute the cost of parenthesization
	my $try=($costs->[$j-$i][$k]+$costs->[$k+1][$j]+
		 $row1*$col1*$col2);


	# if it's the lowest we've seen (or the first),
	# store the cost, dimension, and best
	if (not defined $costs->[$j-$i][$j] or 
	    $try<$costs->[$j-$i][$j]) {
	  $costs->[$j-$i][$j]=$try;
	  $dims->[$j-$i][$j]=[$row1,$col2];
	  $best->[$j-$i][$j]=$k;
	}
      } # end for k
    } # end for j
  } # end for i

  #print "costs=\n",print_multiply_tables($costs);

  #print "best=\n",print_multiply_tables($best);

  print multiply_chain_solve([@matricies],
			      $best,0,
			      scalar @matricies-1);
  print "\n";
			      

}


############################################################




if (not caller) {
  my @matricies=(new matrix(rows=>10,columns=>100),
		 new matrix(rows=>100,columns=>5),
		 new matrix(rows=>5,columns=>50),
		 new matrix(rows=>50,columns=>70),
		 new matrix(rows=>70,columns=>7),
		 new matrix(rows=>7,columns=>10),);
  $_->set(0) foreach @matricies;
  multiply_chain(@matricies);
  
 
}



sub extra_crap {
  my $matrix_1= new matrix(rows=>5,columns=>10);
  my $matrix_2= new matrix(rows=>10,columns=>5);

  for my $row (0..4) {
    for my $column (0..9) {
      $matrix_1->set($row,$column,$row);
      $matrix_2->set($column,$row,$column);
    }
  }

  print "$matrix_1*\n$matrix_2=\n";
  my $result=$matrix_1*$matrix_2;
  print "$result\n";
  print "4*\$result=\n";
  $result*=4;
  print "$result\n\n";

  $matrix_2=new matrix(rows=>5,columns=>10);
  $matrix_2->set(5);
  print "5x10 matrix of 5s:\n$matrix_2\n\n";

  print "$matrix_1-5*\n$matrix_2=\n",$matrix_1-5*$matrix_2,"\n";

  $matrix_1=new matrix(rows=>1,columns=>10);
  $matrix_2=new matrix(columns=>1,rows=>10);
  $matrix_1->set(1);
  $matrix_2->set(1);
  print"$matrix_2*\n$matrix_1=\n",$matrix_2*$matrix_1,"\n\n";
  print"$matrix_1*\n$matrix_2=\n",$matrix_1*$matrix_2,"\n\n";


  $matrix_1=new matrix(rows=>4,columns=>4);
  $matrix_1->set(0,0,2);
  print "matrix_1(0,0)=",$matrix_1->get(0,0),"\n\n";
  $matrix_1->set(0,1,0);
  $matrix_1->set(0,2,0);
  $matrix_1->set(0,3,0);

  $matrix_1->set(1,0,-2);
  $matrix_1->set(1,1,1);
  $matrix_1->set(1,2,-1);
  $matrix_1->set(1,3,0);

  $matrix_1->set(2,0,6);
  $matrix_1->set(2,1,2);
  $matrix_1->set(2,2,1);
  $matrix_1->set(2,3,0);

  $matrix_1->set(3,0,0);
  $matrix_1->set(3,1,0);
  $matrix_1->set(3,2,0);
  $matrix_1->set(3,3,-1);
  print "\n\nmatrix_1:\n$matrix_1\nlu decomposition:\n";
  my $decomposed=$matrix_1->lu_decompose;
  print $decomposed;

  my ($lower,$upper);
  $lower=$matrix_1->clone;
  $upper=$matrix_1->clone;
  $matrix_1->lu_decompose($lower,$upper);
  print "\n\n$lower*\n$upper=\n".$lower*$upper;

}

1;


############################################################




use Inline C => <<'END_C';

#include <math.h>


typedef struct {
  int       rows;
  int       columns;
  // pointer to a ponter of doubles, for a 2 dimensional array
  double**   thedoubles;
} cMatrixDouble;


/***********************************************************
 * Function Name: matrix_allocate
 *
 * Usage: 
 *
 *   $c_data=matrix_allocate(10,10);
 *
 * Description: This function allocates enough space to hold
 *   an NxN matrix of doubles in C. This information is 
 *   returned in a chunk of data which can only be accessed 
 *   by other C functions.
 *
 *  Note that no actual values are set during the allocation.
 *  This could change in the future.
 *
 *   As Jim mahoney describes it in his notes: (some of 
 *   the code was adapted from a similar program of his)
 *
 * Return a perl SV who's IV (integer value) is a pointer
 * to a cMatrixDouble struct.  This technique (from 
 *  Inline::C-Cookbook) embeds an arbitrary C structure with a
 *  single perl scalar that we can then pass back and forth.  
 *  Memory management of these internal structures must be 
 * coded explicity.  
 **********************************************************/
SV* matrix_allocate(int rows,int columns) {
  
  cMatrixDouble*   vector  = malloc(sizeof(cMatrixDouble));
  SV*              data    = newSViv(0);
  double**         nums;

  vector->rows = rows;
  vector->columns = columns;
  /* Allocate in the pointer an array of pointers
     This is useful because you can have each of
     Those pointers represent a row, who then points
     at a list of columns. This is a slighly kluge way
     Of dynamicall allocating an NxN matrix. But it is
     nice because afterwords, we can say nums[row][column] 
  */
  nums = malloc( rows * sizeof(double*) );
  /* If memory did not allocate, error out */
  if( nums == NULL ) {  
    printf("malloc error \n" );
    exit; 
  }
  
  /* Loop over each pointer in nums, pointing it at a row matrix */
  int row_loop;
  for (row_loop=0; row_loop < rows; row_loop++) {
    nums[row_loop] = malloc(columns *sizeof(double) );    
    /* If memory did not allocate, error out */
    if( nums[row_loop]  == NULL ) { 
      fprintf(stdout,"malloc error \n" );
      exit; 
    }    
  }
  
  /* point the pointer thedoubles at the 2 dimensional array */
  vector->thedoubles = nums;
  sv_setiv( data, (IV)vector );
  SvREADONLY_on(data);
  return data;
}


/***********************************************************
 * Function Name: matrix_multiply
 *
 * Usage: 
 *
 *   $product_matrix=new matrix(rows=>$matrix_1->rows,
 *                              columns=>matrix_2->$columns);
 *   matrix_multiply($matrix_1->{_cdata}, $matrix_2->{_cdata},
 *                   $product_matrix->{_cdata} );
 *
 * Description: This function performs matrix multiplication
 *   using the classical definition of matrix multiplication.
 *   It takes in _cdata from 3 matricies. If the dimensions
 *   line up, the third piece of data passed back contains
 *   the product matrix.  
 **********************************************************/
void matrix_multiply(SV* data_1, SV* data_2, SV* resultant_data) {
  cMatrixDouble* matrix_1 = (cMatrixDouble*)SvIV(data_1);
  cMatrixDouble* matrix_2 = (cMatrixDouble*)SvIV(data_2);
  cMatrixDouble* resultant_matrix = (cMatrixDouble*)SvIV(resultant_data);

  /* rows of matrix_1 must equal columns of matrix_2 */
  if (matrix_1->columns != matrix_2->rows) 
    return;
  /* resultant matrix must already be of proper size */
  if ((matrix_1->rows != resultant_matrix->rows) ||
      (matrix_2->columns != resultant_matrix->columns))
    return;
      
  int row_loop;
  int column_loop;
  int current_item;

  /* loop through all rows */
  for (row_loop=0;row_loop<resultant_matrix->rows;row_loop++) {
    /* loop through all columns */
    for (column_loop=0;column_loop<resultant_matrix->columns;column_loop++) {
      
      int value=0;
      int num_of_items_to_multiply=matrix_1->columns;
      /* loop through all items in each row & column */
      for (current_item=0;current_item<num_of_items_to_multiply;current_item++) {
	/* sum the values */
	value+=matrix_1->thedoubles[row_loop][current_item]*
                        matrix_2->thedoubles[current_item][column_loop];
      }
      resultant_matrix->thedoubles[row_loop][column_loop]=value;
    }
  }
  /* resultant_matrix is now modified and passed back */
}


/***********************************************************
 * Function Name: scalar_multiply
 *
 * Usage: 
 *
 *
 *   $scalar_multiple=new matrix(rows=>$matrix->rows,
 *                               columns=>$matrix->columns); 
 *   scalar_multiply($matrix->{_cdata},$scalar,
 *                          $scalar_multiple->{_cdata});
 *   // scalar*$matrix is stored inside $scalar_multiple 
 *
 * Description: This function performs scalar multiplication
 *   on a matrix. It does so by multiplying every item
 *   in the matrix by a scalar. The results are passed 
 *   back as the third parameter.   
 **********************************************************/
void scalar_multiply(SV* data, double scalar, SV* resultant_data) {
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);
  cMatrixDouble* resultant_matrix = (cMatrixDouble*)SvIV(resultant_data);

  /* both matricies must be of same dimension */
  if ((matrix->rows != resultant_matrix->rows) ||
      (matrix->columns != resultant_matrix->columns))
    return;

  int row_loop;
  int column_loop;
  /* loop through all rows */
  for (row_loop=0;row_loop<resultant_matrix->rows;row_loop++) {
    /* loop through all columns */
    for (column_loop=0;column_loop<resultant_matrix->columns;column_loop++) {
      resultant_matrix->thedoubles[row_loop][column_loop]=
	scalar*matrix->thedoubles[row_loop][column_loop];
    }
  }
  /* resultant_matrix is now modified and passed back */
}


/***********************************************************
 * Function Name: matrix_add
 *
 * Usage: 
 *
 *   $sum=new matrix(rows=>$matrix->rows,
 *                   columns=>$matrix->columns); 
 *   matrix_add($matrix_1->{_cdata},$matrix_2->{_cdata},
 *              $sum->{_cdata});
 *   // $matrix_1+$matrix_2 is stored inside $sum 
 *
 * Description: Adds up every term from the first 2 matricies
 *   and returns the sum in the 3rd matrix. All 3 matricies
 *   must be of the same size.
 **********************************************************/
void matrix_add(SV* data_1, SV* data_2, SV* resultant_data) {
  cMatrixDouble* matrix_1 = (cMatrixDouble*)SvIV(data_1);
  cMatrixDouble* matrix_2 = (cMatrixDouble*)SvIV(data_2);
  cMatrixDouble* resultant_matrix = (cMatrixDouble*)SvIV(resultant_data);

  /* all matricies must be of same dimension */
  if ((matrix_1->rows != matrix_2->rows) ||
      (matrix_1->columns != matrix_2->columns) ||
      (matrix_1->rows != resultant_matrix->rows) ||
      (matrix_1->columns != resultant_matrix->columns))
    return;

  int row_loop;
  int column_loop;
  /* loop through all rows */
  for (row_loop=0;row_loop<resultant_matrix->rows;row_loop++) {
    /* loop through all columns */
    for (column_loop=0;column_loop<resultant_matrix->columns;column_loop++) {
      resultant_matrix->thedoubles[row_loop][column_loop]=
	matrix_1->thedoubles[row_loop][column_loop]+
	matrix_2->thedoubles[row_loop][column_loop];;
    }
  }
  /* resultant_matrix is now modified and passed back */
}


/***********************************************************
 * Function Name: matrix_add
 *
 * Usage: 
 *
 *
 *   $difference=new matrix(rows=>$matrix->rows,
 *                   columns=>$matrix->columns); 
 *   matrix_subtract($matrix_1->{_cdata},$matrix_2->{_cdata},
 *              $difference->{_cdata});
 *   // $matrix_1-$matrix_2 is stored inside $difference 
 *
 * Description: subtracts every term from the first 2 matricies
 *   and returns the difference in the 3rd matrix. All 3 
 *   matricies must be of the same size.
 **********************************************************/
void matrix_subtract(SV* data_1, SV* data_2, SV* resultant_data) {
  cMatrixDouble* matrix_1 = (cMatrixDouble*)SvIV(data_1);
  cMatrixDouble* matrix_2 = (cMatrixDouble*)SvIV(data_2);
  cMatrixDouble* resultant_matrix = (cMatrixDouble*)SvIV(resultant_data);

  /* all matricies must be of same dimension */
  if ((matrix_1->rows != matrix_2->rows) ||
      (matrix_1->columns != matrix_2->columns) ||
      (matrix_1->rows != resultant_matrix->rows) ||
      (matrix_1->columns != resultant_matrix->columns))
    return;

  int row_loop;
  int column_loop;
  /* loop through all rows */
  for (row_loop=0;row_loop<resultant_matrix->rows;row_loop++) {
    /* loop through all columns */
    for (column_loop=0;column_loop<resultant_matrix->columns;column_loop++) {
      resultant_matrix->thedoubles[row_loop][column_loop]=
	matrix_1->thedoubles[row_loop][column_loop]-
	matrix_2->thedoubles[row_loop][column_loop];;
    }
  }
  /* resultant_matrix is now modified and passed back */
}


/***********************************************************
 * Function Name: matrix_get
 *
 * Usage: 
 *
 *   $item=matrix_gete($matrix->{_cdata}, $row, $column);
 *
 * Description: 
 *
 *   If you have a matrix held in a variable like cdata (which
 *   is used only by being passed to C functions) 
 *   you can call this function with a row and column
 *   to get back that value in the array.
 **********************************************************/
double  matrix_get( SV* data, int row,int column ){
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);
  row = row % matrix->rows;
  column=column % matrix->columns;
  return matrix->thedoubles[row][column];
}


/***********************************************************
 * Function Name: matrix_set
 *
 * Usage: 
 *
 *   matrix_set($matrix->{_cdata}, $row, $column,$value);
 *
 * Description: 
 *
 *   If you have a matrix held in a variable like cdata (which
 *   is used only by being passed to C functions), 
 *   you can call this function with a row and column and
 *   value to set that spot to the given value.
 **********************************************************/
double  matrix_set( SV* data, int row,int column, double value){
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);
  row = row % matrix->rows;
  column=column % matrix->columns;
  matrix->thedoubles[row][column] = value;
  return value;
}


/***********************************************************
 * Function Name: matrix_rows
 *
 * Usage: 
 * 
 *   $rows=matrix_rows($matrix->{_cdata});
 *
 * Description: Returns the number of rows in the matrix
 **********************************************************/
int  matrix_rows( SV* data) {
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);
  return matrix->rows;
}


/***********************************************************
 * Function Name: matrix_columns
 *
 * Usage: 
 * 
 *   $rows=matrix_columns($matrix->{_cdata});
 *
 * Description: Returns the number of columns in the matrix
 **********************************************************/
int  matrix_columns( SV* data) {
  cMatrixDouble* matrix = (cMatrixDouble *)SvIV(data);
  return matrix->columns;
}


/***********************************************************
 * Function Name: matrix_set_all_values
 *
 * Usage: 
 * 
 *   matrix_set_all_values($matrix->{_cdata},$value);
 *
 * Description: sets every value in the matrix to
 *   what is inputted
 **********************************************************/
int matrix_set_all_values(SV* data, double value) {
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);

  int row_loop;
  int column_loop;

  /* loop through all rows */
  for (row_loop=0;row_loop<matrix->rows;row_loop++) {
    /* loop through all columns */
    for (column_loop=0;column_loop<matrix->columns;column_loop++) {
      matrix->thedoubles[row_loop][column_loop]=value;
    }
  }
}


/***********************************************************
 * Function Name: matrix_lu_decompose
 *
 * Usage: 
 * 
 *   matrix_lu_decompose($matrix->{_cdata},$value);
 *
 * Description: Performs the LU decomposition using 
 *   crout's algorithm.
 **********************************************************/
void matrix_lu_decompose(SV* data) {
  cMatrixDouble* matrix = (cMatrixDouble*)SvIV(data);


  /* both matricies must be of same dimension */
  if (matrix->rows!=matrix->columns) return;

  // i=rows
  // j=columns

  int column_loop;
  int row_loop;
  int sum;
  int k;

  // loop through all columns
  for (column_loop=0;column_loop<matrix->columns;column_loop++) {
    // loop finding betas.
    for (row_loop=0;row_loop<=column_loop;row_loop++) {
      sum = matrix->thedoubles[row_loop][column_loop];
      
      for (k=1;k<row_loop;k++) 
	sum -= matrix->thedoubles[row_loop][k]*matrix->thedoubles[k][column_loop];
      matrix->thedoubles[row_loop][column_loop]=sum;

    }
    // loop finding alphas
    for (row_loop=column_loop+1;row_loop<matrix->columns;row_loop++) {
      sum = matrix->thedoubles[row_loop][column_loop];
      

      for (k=1;k<column_loop;k++) 
	sum -= matrix->thedoubles[row_loop][k]*matrix->thedoubles[k][column_loop];
      sum = sum / matrix->thedoubles[column_loop][column_loop];
      matrix->thedoubles[row_loop][column_loop]=sum; 
    }
  }
}



END_C


=head1 NAME

Matrix - An object orientend framework for creating and using matricies.

=head1 SYNOPSIS

  use matrix;
  my $A = new matrix(rows=>2,columns=>2);
  $A->set(0); #sets all values to 0
  my $B = new matrix(rows=>2,columns=>2);

  print "A=$A\n";
  print "B=$B\n";
  print "5*A+A*B=",5*$A+$A*$B;

=head1 DESCRIPTION

The following methods are available:

=head2 new

Constructor method which creates a new matrix. Arguments: Number of rows and columns in the matrix. 

  my $matrix=new matrix(rows=>3,columns=>3);


=head2 display

To display a matrix:

  print $matrix,"\n; 
  print "$matrix\n";


=head2 columns

columns returns the number of columns in the matrix

  my $number_of_columns=$matrix->columns;


=head2 rows

rows returns the number of rows in the matrix.

  my $number_of_rows=$matrix->rows;


=head2 set

Set is used to set values of elements. It can be called in 3 distinct contexts to perform 3 distinct functions. 

If you want to set one individual item in the matrix, pass the row, column, and value.

  $matrix->set(0,0,5) # sets the element (0,0) to 5

If you want to set all values at once

  $matrix->set(5) # sets all elements to 5

If you want to set several distinct values at once, pass in data structure:

  $matrix->set([[1,0],[0,1]]); # sets the matrix to the identity matrix


=head2 get

Get the value of an element in the matrix.

  my $value=$matrix->get(0,0) #get value in element (0,0)


=head2 add

This is a method for adding 2 equal dimension matricies. The eaisyest way to do this is as follows

  $sum_matrix=$matrix_1 + $matrix_2;


=head2 subtract

Basically the same as add, but with a minus sign.

  my $difference=$matrix_1 - $matrix2;

  # is equivalent to
  my $difference=$matrix_1 + (-1)*$matrix_2;


=head2 multiply

multiply can multiply either 2 matricies or a matrix and a scalar. The syntax for matrix multiplication is as follows:

  $result=$matrix_1 * $matrix_2.

The syntax for scalar multiplication is
 
  my $result=$scalar*$matrix;
  # or
  my $result=$matrix*$scalar;


=head1 AUTHOR

Joshua Lande, (Jolande@marlboro.edu)

=head1 BUGS

None that I can find. Please contact me if you find any

=head1 SEE ALSO

Math::MatrixReal

PDL

  http://search.cpan.org/~csoe/PDL-2.4.2/Basic/PDL.pm

Math::MatrixReal

  http://search.cpan.org/~leto/Math-MatrixReal-1.9/MatrixReal.pm


=head1 COPYRIGHT

  This program is free software. You may copy or redistribute it under the same terms as perl itself.
