#!/usr/bin/perl -w
# test_matrix.pl  uses  matrix.pm   in local directory
use strict;

print "test_matrix.pl running \n";

use matrix qw(&matadd &matsub &matmul &matinit &matrc &vecinit &matident
              &matrand &vecrand &matprt &vecprt &materr &vecerr
              &matcpy &veccpy &simeq &matinvert &matvecmul
              &vecadd &vecsub &vecdot &veccross &nuderiv);

my @mat1 = ([1, 2], [3, 4], [5, 6]);
my @mat2 = ([10, 11, 12], [13, 14, 15]);
my @mat3 = ([2, 3], [4, 5], [6, 7]);
my @sum =  ([0, 0], [0, 0], [0, 0]);
my @vec1 = (1, 2, 3);
my @vec2 = (2, 3, 4);
my @vec3;
vecinit(3, \@vec3);

# my @prod = ([0, 0, 0],[0, 0, 0],[0, 0, 0]);
my @prod;
matinit(3, 3, \@prod);

my @ident;
matident(3, 3, \@ident);

srand(time);
my @matr;
matrand(3, 3, \@matr);
my @vecr;
vecrand(3, \@vecr);

# my @mat3c = @mat3; # try copy
# $mat3[1][1] = 999; # no, just copy reference 
# $mat3c[0][0] = 111; # changes both names
# my $mat3c = $mat3; # no, need matcpy
my @mat3c;
matcpy(\@mat3, \@mat3c);
$mat3c[0][0] = 111;

my $max = 0;
my $rms = 0;

print "built mat1, mat2, mat3, sum, prod, ident \n";
matprt("mat1", \@mat1);
print "\n";
matprt("mat2", \@mat2);
print "\n";
matprt("mat3", \@mat3);
print "\n";
matprt("mat3c", \@mat3c);
print "\n";
matprt("matr", \@matr);
print "\n";
vecprt("vecr", \@vecr);
print "\n";
matprt("ident", \@ident);
print "\n";


my $nr1;
my $nc1;
($nr1, $nc1) = matrc(\@mat1);
print "nr1=$nr1, nc1=$nc1 \n";
my $nr2;
my $nc2;
($nr2, $nc2) = matrc(\@mat2);
print "nr2=$nr2, nc2=$nc2 \n";
print "\n";

matadd(\@mat1, \@mat3, \@sum);
print "matadd sum= \n";
matprt("sum", \@sum);
print "\n";

my @dif;
matsub(\@mat3, \@mat1, \@dif);
print "matsub \@ dif= \n";
matprt("dif", \@dif);
print "\n";

matmul(\@mat1, \@mat2, \@prod);
print "matmul \@ prod= \n";
matprt("prod", \@prod);
print "\n";

($max, $rms) = materr(\@prod, \@ident);
print "prod-ident max=$max, rms=$rms \n";
print "\n";

my @invr=([(0)x3],[(0)x3],[(0)x3]);
matinvert(\@matr, \@invr);
print "invert random matr \n";
matprt("invr", \@invr);
print "\n";

print "check matr * invr = prod == ident \n";
matmul(\@matr, \@invr, \@prod);
matprt("inv prod", \@prod);
($max, $rms) = materr(\@ident, \@prod);
print "ident-prod max=$max, rms=$rms \n";
print "\n";
# try Hilbert matrix A[i][j] = 1.0/(i+j+2)

my @vecx=(0,0,0);
my @vecyc=(0,0,0);
simeq(\@matr, \@vecr, \@vecx);
print "simeq random matr vecr vecx \n";
vecprt("vecx", \@vecx);
print "\n";
matvecmul(\@matr, \@vecx, \@vecyc);
vecprt("vecyc", \@vecyc);
($max, $rms) = vecerr(\@vecr, \@vecyc);
print "vecr-vecyc max=$max, rms=$rms \n";
print "\n";

vecadd(\@vec1, \@vec2, \@vec3);
print "vec1 + vec2 = vec3 \n";
vecprt("vec1+vec2", \@vec3);
print "\n";

vecsub(\@vec1, \@vec2, \@vec3);
print "vec1 + vec2 = vec3 \n";
vecprt("vec1-vec2", \@vec3);
print "\n";

my $dotp = vecdot(\@vec1, \@vec2);
print "vec1 dot vec2 = $dotp \n";
print "\n";

my @crossp;
veccross(\@vec1, \@vec2, \@crossp);
print "vec1 cross vec2 = crossp \n";
vecprt("crossp", \@crossp);
print "\n";
print "\n";

my @x;
my @c;
my $order=2;
my $npoint=5;
my $point =1;
for(my $i=0; $i<$npoint; $i++) { $x[$i] = $i; }
nuderiv($order, $npoint, $point, \@x, \@c);
print "test nuderiv x spacing, c derivative corfficients \n";
print "order= $order, npoint= $npoint, point= $point \n";
vecprt("x", \@x);
print "\n";
vecprt("c", \@c);
print "\n";

print "test_matrix.pl finished \n";

