#!/usr/bin/perl

$M1=3;
$M2=3;
if ($#ARGV>=0) {$M1 = $M2 = $ARGV[0];}
if ($#ARGV>=1) {$M2 = $ARGV[1];}


# $N1m = $M;
# $N2m = $M;
# $N3m = $M;

print "\#include \"sfunction.h\"\n\n";

for ($N1=1; $N1<=$M1; $N1++){
    for ($N2=1; $N2<=$M2; $N2++){
	for ($N3=1; $N3<=$M1; $N3++){

	    $ii = $N1 + 16*($N2 + 16*$N3);
	    $alredy{$ii}++;
	    @a = ($N1,$N2,$N3);
	    $case{$ii} =  [($N1,$N2,$N3)];
	    
	    print "template<>\n";
	    print "void Multiply_<$N1,$N2,$N3>(function2D<double>& C, const function2D<double>& A, const function2D<double>& B)\n";
	    print "{\n";
	    for ($i=0; $i<$N1; $i++){
		for ($j=0; $j<$N3; $j++){
		    $c0 = $i*$N3+$j;
		    print "  C($i,$j) = ";
		    for ($k=0; $k<$N2-1; $k++) {
			$a0 = $i*$N2+$k;
			$b0 = $k*$N3+$j;
			print "A($i,$k)*B($k,$j)+";
		    }
		    $k = $N2-1;
		    $a0 = $i*$N2+$k;
		    $b0 = $k*$N3+$j;
		    print "A($i,$k)*B($k,$j);";
		    print "\n";
		}
	    }
	    print "}\n";
	}
    }
}
for ($N1=1; $N1<=$M2; $N1++){
    for ($N2=1; $N2<=$M1; $N2++){
	for ($N3=1; $N3<=$M2; $N3++){
    
	    $ii = $N1 + 16*($N2 + 16*$N3);
	    if ($alredy{$ii}) {next;}
	    @a = ($N1,$N2,$N3);
	    $case{$ii} =  [($N1,$N2,$N3)];
	    
	    print "template<>\n";
	    print "void Multiply_<$N1,$N2,$N3>(function2D<double>& C, const function2D<double>& A, const function2D<double>& B)\n";
	    print "{\n";
	    for ($i=0; $i<$N1; $i++){
		for ($j=0; $j<$N3; $j++){
		    $c0 = $i*$N3+$j;
		    print "  C($i,$j) = ";
		    for ($k=0; $k<$N2-1; $k++) {
			$a0 = $i*$N2+$k;
			$b0 = $k*$N3+$j;
			print "A($i,$k)*B($k,$j)+";
		    }
		    $k = $N2-1;
		    $a0 = $i*$N2+$k;
		    $b0 = $k*$N3+$j;
		    print "A($i,$k)*B($k,$j);";
		    print "\n";
		}
	    }
	    print "}\n";
	}
    }
}


print "\n\n";
print "void Multiply(function2D<double>& C, const function2D<double>& A, const function2D<double>& B)\n";
print "{\n";
print "  int N = A.size_N();\n";
print "  int K = A.size_Nd();\n";
print "  int M = B.size_Nd();\n";
print "  int Kp = B.size_N();\n";
print "  int ii = N + 16*(K + M*16);\n";
print "  C.N = N; C.Nd = M;\n";
print "  if (K!=Kp) {std::cerr<<\"Matrix sizes not correct in Multiply matrices!\"<<std::endl;}\n";
print "  switch(ii){\n";

foreach $f (sort { $a <=> $b } keys %case){
    @Ns = @{$case{$f}};
    print "  case $f : Multiply_<$Ns[0],$Ns[1],$Ns[2]>(C,A,B); break;\n";
}
print "  default :   C.MProduct(A,B); break;\n";
print "  }\n";
print "}\n";