#!/usr/bin/perl -w

use strict;
use Getopt::Std;
use Math::MatrixReal;
my (%opts, $i, $j, $A, $x, $c, $v);

%opts = (
    f => "%+5.1f",	# format string
    m =>  3,	# height
    n =>  3,	# width
    r =>  undef,# rank
    l => -9,	# lower bound
    u =>  9,	# upper bound
    v => undef,	# variable names
);

print "# options: @ARGV\n";
getopts("f:m:n:r:l:u:v:", \%opts);
unless (defined $opts{v} and length($opts{v})==$opts{n}) {
    if ($opts{n} <= 3) {
	$opts{v} = substr("xyz", 0, $opts{n});
    } elsif ($opts{n} <= 6) {
	$opts{v} = substr("uvw", 0, $opts{n}-3) . "xyz";
    } else {
	$opts{v} = map { chr(ord('a')+$_) } 0..$opts{n}-1
    }
}
@$v = (undef, split //, $opts{v});

$A = rand_mat($opts{m}, $opts{m}) *  lower_rank($opts{m}, $opts{n}, $opts{r});
$x = rand_mat($opts{n}, 1);
$c = $A * $x;
print "#";
for ($j=1; $j<=$opts{n}; ++$j) {
    printf "  $opts{f} ", $x->element($j, 1);
}
print "\n";
for ($i=1; $i<=$opts{m}; ++$i) {
    for ($j=1; $j<=$opts{n}; ++$j) {
	printf " $opts{f} $v->[$j]", $A->element($i, $j);
    }
    printf " = $opts{f}\n", $c->element($i, 1);
}

sub rndint {
    return int(rand($opts{u}-$opts{l}+1))+$opts{l};
}

sub lower_rank {
    my ($m, $n, $r) = @_;
    my ($i, $j, $LC);
    $r = $m > $n ? $n : $m unless defined $r;
    die if $r > $m or $r > $n;
    $r = [ sort @{ rand_permute($m) }[0..$r-1] ];
    $i = 0;
    for ($j=0; $j<$n; ++$j) {
	my ($k);
	for ($k=0; $k<$i; ++$k) {
	    $LC->[$k][$j] = int(rand(5))-2;
	}
	for ($k=$i; $k<$m; ++$k) {
	    $LC->[$k][$j] = 0;
	}
	$LC->[$i++][$j] = 1 if ($i<=$#$r and $j == $r->[$i]);
    }
    $LC = Math::MatrixReal->new_from_rows($LC);
    return $LC;
}

sub rand_mat {
    my ($m, $n) = @_;
    my ($i, $j, $mat);
    for ($j=0; $j<$m; ++$j) {
	for ($i=0; $i<$m; ++$i) {
	    $mat->[$i][$j] = rndint();
	}
    }
    return Math::MatrixReal->new_from_rows($mat);
}

sub rand_permute {
    my ($n) = @_;
    my ($i, $r);
    [ map { $_->[0] } sort { $a->[1]<=>$b->[1] } map { [$_, rand()] } 0..$n-1 ];
}