package matrix;

use strict;
# use messages; # strange
use Exporter qw(import);

our $VERSION     = 1.00;
our @ISA         = qw(Exporter);
our @EXPORT_OK   = qw(matadd matsub matmul matrc matinit vecinit matident
                      matrand vecrand matprt vecprt materr vecerr
                      matcpy veccpy simeq matinvert matvecmul
                      vecadd vecsub vecdot veccross nuderiv);

sub matadd {
    my ($mat1, $mat2, $sum) = @_; # parameters
    my $r1 = @$mat1;
    my $c1 = @{$mat1->[0]};
    my $r2 = @$mat2;
    my $c2 = @{$mat2->[0]};
    print "in matadd r1=$r1, c1=$c1, r2=$r2, c2=$c2 \n";

    for (my $i = 0; $i < $r1; $i++) {
        for (my $j = 0; $j < $c1; $j++) {
           $sum->[$i][$j] = $mat1->[$i][$j] + $mat2->[$i][$j];
        }
    }
} # end matadd

sub matsub {
    my ($mat1, $mat2, $sum) = @_; # parameters
    my $r1 = @$mat1;
    my $c1 = @{$mat1->[0]};
    my $r2 = @$mat2;
    my $c2 = @{$mat2->[0]};
    print "in matsub r1=$r1, c1=$c1, r2=$r2, c2=$c2 \n";

    for (my $i = 0; $i < $r1; $i++) {
        for (my $j = 0; $j < $c1; $j++) {
           $sum->[$i][$j] = $mat1->[$i][$j] - $mat2->[$i][$j];
        }
    }
} # end matsub

sub matmul {
    my ($mat1, $mat2, $product) = @_; # parameters
    my $r1 = @$mat1;
    my $c1 = @{$mat1->[0]};
    my $r2 = @$mat2;
    my $c2 = @{$mat2->[0]};
    print "in matmul r1=$r1, c1=$c1, r2=$r2, c2=$c2 \n";

    for (my $i = 0; $i < $r1; $i++) {
        for (my $j = 0; $j < $c2; $j++) {
            my $sum = 0;
            for (my $k = 0; $k < $c1; $k++) {
                $sum += $mat1->[$i][$k]* $mat2->[$k][$j];
            }
            $product->[$i][$j] = $sum;
        }
    }
} # end matmul

sub matrc {  # return number of rows and columns
    my ($mat1)  = @_;
    my $num_rows = @$mat1;
    my $num_cols =  @{$mat1->[0]};
    # Assume all rows have an equal number of columns.
    return ($num_rows, $num_cols);
} # end matrc

sub matinit {
    my ($nrow, $ncol, $mat) = @_; # parameters; # my @mat; then matinit
    print "in matinit nrow=$nrow, ncol=$ncol \n";

    for (my $i = 0; $i < $nrow; $i++) {
        for (my $j = 0; $j < $ncol; $j++) {
           $mat->[$i][$j] = 0;
        }
    }
} # end matinit

sub matident {
    my ($nrow, $ncol, $mat) = @_; # parameters; # my @mat; then matident
    print "in matident nrow=$nrow, ncol=$ncol \n";

    for (my $i = 0; $i < $nrow; $i++) {
        for (my $j = 0; $j < $ncol; $j++) {
           $mat->[$i][$j] = 0;
           if ($i == $j) {$mat->[$i][$j] = 1;}
        }
    }
} # end matident

sub vecinit {
    my ($n, $vec) = @_; # parameters; # my @vec; then vecinit
    print "in vecinit n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
       $vec->[$i] = 0;
    }
} # end vecinit

sub matrand {
    my ($nrow, $ncol, $mat) = @_; # parameters; # my @mat; then matinit
    print "in matrand nrow=$nrow, ncol=$ncol \n";
    for (my $i = 0; $i < $nrow; $i++) {
        for (my $j = 0; $j < $ncol; $j++) {
           $mat->[$i][$j] = rand;
        }
    }
} # end matrand

sub vecrand {
    my ($n, $vec) = @_; # parameters; # my @vec; then vecinit
    print "in vecrand n=$n \n";
    for (my $i = 0; $i < $n; $i++) {
       $vec->[$i] = rand;
    }
} # end vecand

sub matprt {
    my ($name, $mat) = @_; # parameters
    my $nrow = @$mat;
    my $ncol = @{$mat->[0]};
    print "in matprt nrow=$nrow, ncol=$ncol \n";

    for (my $i = 0; $i < $nrow; $i++) {
        for (my $j = 0; $j < $ncol; $j++) {
           print "$name \[ $i \] \[ $j \] = $mat->[$i][$j] \n";
        }
    }
} # end matprt

sub vecprt {
    my ($name, $vec) = @_; # parameters
    my $n = @$vec;
    print "in vecprt n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
       print "$name \[ $i \] = $vec->[$i] \n";
    }
} # end vecprt

sub materr {
    my ($mat1, $mat2) = @_; # parameters
    my $nrow = @$mat1;
    my $ncol = @{$mat1->[0]};
    my $n = $nrow * $ncol;
    my $sum = 0;
    my $sumsq = 0;
    my $err = 0;
    my $max = 0;
    my $rms = 0;
    print "in materr nrow=$nrow, ncol=$ncol \n";

    for (my $i = 0; $i < $nrow; $i++) {
      for (my $j = 0; $j < $ncol; $j++) {
        $err = abs($mat1->[$i][$j]-$mat2->[$i][$j]);
        $sum = $sum + $err;
        $sumsq = $sumsq + $err*$err;
        if ($err > $max) { $max = $err; }
      }
    }
    $rms = sqrt(($sumsq - $sum*$sum/$n)/$n);
    return ($max, $rms);
} # end materr

sub vecerr {
    my ($vec1, $vec2) = @_; # parameters
    my $n = @$vec1;
    my $sum = 0;
    my $sumsq = 0;
    my $err = 0;
    my $max = 0;
    my $rms = 0;
    print "in vecerr n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
       $err = abs($vec1->[$i]-$vec2->[$i]);
       $sum = $sum + $err;
       $sumsq = $sumsq + $err*$err;
       if ($err > $max) { $max = $err; }
    }
    $rms = sqrt(($sumsq - $sum*$sum/$n)/$n);
    return ($max, $rms);
} # end vecerr

sub matcpy {
    my ($mat1, $mat2) = @_; # parameters
    my $nrow = @$mat1;
    my $ncol = @{$mat1->[0]};
    print "in matcpy nrow=$nrow, ncol=$ncol \n";

    for (my $i = 0; $i < $nrow; $i++) {
        for (my $j = 0; $j < $ncol; $j++) {
           $mat2->[$i][$j] = $mat1->[$i][$j];
        }
    }
} # end matcpy

sub veccpy {
    my ($vec1, $vec2) = @_; # parameters
    my $n = @$vec1;
    print "in veccpy n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
      $vec2->[$i] = $vec1->[$i];
    }
} # end veccpy
 
sub matinvert {
    my ($mat1, $mat2) = @_; # parameters
    my $n = @$mat1; # must be square
    my @krow;
    my @kcol;
    my @temp;
    my $pivot = 0;
    my $i_pivot = 0;
    my $j_pivot = 0;
    my $abs_pivot = 0;
    my $hold = 0;
    print "in matinvert n by n = $n \n";
    for (my $i = 0; $i < $n; $i++) {
      for (my $j = 0; $j < $n; $j++) {
	$mat2->[$i][$j] = $mat1->[$i][$j];
      }
    }
    for (my $k = 0; $k < $n; $k++) {
       $krow[$k] = $k;
       $kcol[$k] = $k;
       $temp[$k] = 0.0;
    }
    # begin main reduction loop
    for (my $k = 0; $k < $n; $k++) {
      $pivot = $mat2->[$krow[$k]][$kcol[$k]];
      $i_pivot = $k;
      $j_pivot = $k;
      for (my $i = $k; $i < $n; $i++) {
        for (my $j = $k; $j < $n; $j++) {
	  $abs_pivot = abs($pivot);
          if (abs($mat2->[$krow[$i]][$kcol[$j]]) > $abs_pivot) {
	    $i_pivot = $i;
	    $j_pivot = $j;
            $pivot = $mat2->[$krow[$i]][$kcol[$j]];
	  } # end if
	} # end i
      } # end j
      $hold = $krow[$k];
      $krow[$k] = $krow[$i_pivot];
      $krow[$i_pivot] = $hold;
      $hold = $kcol[$k];
      $kcol[$k] = $kcol[$j_pivot];
      $kcol[$j_pivot] = $hold;
      # reduce about pivot
      $mat2->[$krow[$k]][$kcol[$k]] = 1.0/$pivot;
      for (my $j = 0; $j < $n; $j++) {
	if ($j != $k) {
          $mat2->[$krow[$k]][$kcol[$j]] = $mat2->[$krow[$k]][$kcol[$j]]*
                                          $mat2->[$krow[$k]][$kcol[$k]];
	}
      }
      # inner reduction loop
      for (my $i = 0; $i < $n; $i++) {
	if ($k != $i) {
          for (my $j = 0; $j < $n; $j++) {
	    if($k != $j) {
              $mat2->[$krow[$i]][$kcol[$j]] = $mat2->[$krow[$i]][$kcol[$j]]-
                                              $mat2->[$krow[$i]][$kcol[$k]]* 
                                              $mat2->[$krow[$k]][$kcol[$j]];

	    } # end if j
	  } # end j
          $mat2->[$krow[$i]][$kcol[$k]] = - $mat2->[$krow[$i]][$kcol[$k]]*
                                          $mat2->[$krow[$k]][$kcol[$k]];
	} # end if i
      } # end i
    } # end k main
    # unscramble rows
    for (my $j = 0; $j < $n; $j++) {
      for (my $i = 0; $i < $n; $i++) {
        $temp[$kcol[$i]] = $mat2->[$krow[$i]][$j];
      } # end i
      for (my $i = 0; $i < $n; $i++) {
         $mat2->[$i][$j] = $temp[$i];
      } # end i
    } # end j
    # unscramble columns
    for (my $i = 0; $i < $n; $i++) {
      for (my $j = 0; $j < $n; $j++) {
        $temp[$krow[$j]] = $mat2->[$i][$kcol[$j]];
      } # end j
      for (my $j = 0; $j < $n; $j++) {
         $mat2->[$i][$j] = $temp[$j];
      } # end j
    } # end i
} # end matinvert

sub simeq {
    # solve linear equations for X where Y = A * X
    # method: Gauss-Jordan elimination using maximum pivot
    # usage:  simeq(\@A, \@Y, \@X);
    #    Translated to perl by : Jon Squire , 20 August 2016
    #    First written by Jon Squire December 1959 for IBM 650, translated to
    #    other languages  e.g. Fortran converted to Ada converted to C
    #    converted to java then converted to perl

    my ($mat, $yvec, $xvec) = @_; # parameters solve A x = y
    my $n = @$mat; # must be square
    my $m = $n + 1;
    my $B->[$n-1][$n] = 0;
    my @krow;
    my $pivot = 0;
    my $i_pivot = 0;
    my $abs_pivot = 0;
    my $hold = 0;
    print "in simeq n by n = $n \n";
    for ( my $i = 0; $i < $n; $i++) {
      for ( my $j = 0; $j < $n; $j++) {
        $B->[$i][$j] = $mat->[$i][$j];
      } 
      $B->[$i][$n] = $yvec->[$i];
    }
    # set up row interchange
    for (my $k = 0; $k < $n; $k++) {
       $krow[$k] = $k;
    }
    # begin main reduction loop
    for (my $k = 0; $k < $n; $k++) {
      $pivot = $B->[$krow[$k]][$k];
      $abs_pivot = abs($pivot);
      $i_pivot = $k;
      for (my $i = $k; $i < $n; $i++) {
	$abs_pivot = abs($pivot);
        if (abs($B->[$krow[$i]][$k]) > $abs_pivot) {
	  $i_pivot = $i;
          $pivot = $B->[$krow[$i]][$k];
	  $abs_pivot = abs($pivot);
	} # end if
      } # end i
      $hold = $krow[$k];
      $krow[$k] = $krow[$i_pivot];
      $krow[$i_pivot] = $hold;
      # reduce about pivot
      for (my $j = $k+1; $j < $m; $j++) {
	$B->[$krow[$k]][$j] = $B->[$krow[$k]][$j]/$B->[$krow[$k]][$k];
      } # end j
      # inner reduction loop
      for (my $i = 0; $i < $n; $i++) {
	if ($i != $k) {
          for (my $j = $k+1; $j < $m; $j++) {
	    $B->[$krow[$i]][$j] = $B->[$krow[$i]][$j] -
                                  $B->[$krow[$i]][$k]*$B->[$krow[$k]][$j];
	  } # end j
	} # end if
      } # end i
    } # end k
    # build xvec
    for (my $i = 0; $i < $n; $i++) {
      $xvec->[$i] = $B->[$krow[$i]][$n];
    }
} # end simeq

sub matvecmul  {
    my ($mat, $vec, $product) = @_; # parameters
    my $r1 = @$mat;
    my $c1 = @{$mat->[0]};
    print "in matvecmul r1=$r1, c1=$c1 \n";

    for (my $i = 0; $i < $r1; $i++) {
      my $sum = 0;
      for (my $j = 0; $j < $c1; $j++) {
        $sum += $mat->[$i][$j] * $vec->[$j];
      } # end j
      $product->[$i] = $sum;
    } # end i
} # end matvecmul

sub vecadd {
    my ($vec1, $vec2, $sum) = @_; # parameters
    my $n = @$vec1;
    print "in vecadd n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
      $sum->[$i] = $vec1->[$i] + $vec2->[$i];
    }
} # end vecadd

sub vecsub {
    my ($vec1, $vec2, $dif) = @_; # parameters
    my $n = @$vec1;
    print "in vecsub n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
      $dif->[$i] = $vec1->[$i] - $vec2->[$i];
    }
} # end vecsub

sub vecdot {
    my ($vec1, $vec2) = @_; # parameters
    my $n = @$vec1;
    my $sum = 0;
    print "in vecdot n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
      $sum += $vec1->[$i] * $vec2->[$i];
    }
    return $sum;
} # end vecdot

sub veccross {
    my ($vec1, $vec2, $vec3) = @_; # parameters
    my $n = @$vec1;
    my $sign = 1.0;
    my $k = 0;
    print "in veccross n=$n \n";

    for (my $i = 0; $i < $n; $i++) {
      $vec3->[$i] = 0.0;
      for (my $j = 0; $j < $n; $j++) {
	$k = $j + $i;
        if ( $k >= $n ) { $k = $k - $n; }
          $vec3->[$i] += $vec1->[$j] * $vec2->[$k];
      } # end j
      $vec3->[$i] = $sign * $vec3->[$i];
      $sign = - $sign;
    } # end i
} # end veccross

sub nuderiv { # see nuderiv.c for comments
  my ($order, $npoint, $point, $x, $c) = @_; # paramenters 

  my $n = $npoint;
  my $A->[$n-1][$n-1] = 0;
  my $AI->[$n-1][$n-1] = 0;
  my $fct->[$n-1] = 0;
  my $pwr;
  # A = (double *)calloc(n*n, sizeof(double));
  # AI = (double *)calloc(n*n, sizeof(double));
  # fct = (double *)calloc(n, sizeof(double));

  for(my $i=0; $i<$n; $i++) { # build matrix that will be inverted 
    $pwr=1.0;
    for(my $j=0; $j<$n; $j++) {
      $A->[$i][$j] = $pwr;
      $pwr = $pwr*$x->[$i];
    } # end j
  } # end i
  matinvert($A, $AI); # invert the matrix
  # form derivative for order and point
  for(my $i=0; $i<$n; $i++) { # compute factors for order
    $fct->[$i] = 1.0;
  }
  for(my $j=0; $j<$order; $j++) {
    for(my $i=0; $i<$n; $i++) {
      $fct->[$i] = $fct->[$i]*($i-$j);
    } # end i
  } # end j
  for(my $i=0; $i<$n; $i++) {
    $c->[$i] = 0.0;
    $pwr = 1.0;
    for(my $j=$order; $j<$n; $j++) {
      $c->[$i] = $c->[$i] + $fct->[$j]*$AI->[$j][$i]*$pwr;
      $pwr = $pwr * $x->[$point];
    } # end j
  } # end i
} # end nuderiv

1;