1#! /bin/perl 2# 3# From: Reuben Settergren <rjs@jhu.edu> 4# 5# I wrote a perl script to randomly generate matrices with smallish integers, 6# and test for a suitable determinant (typically +/-1), and ouput the matrix 7# and its inverse. 8# 9# There are switches to make it -s[ymmetric] and/or -d[iagonally strong], 10# which are necessary and good properties for covariance matrices. You may 11# need to grab Math::MatrixReal and Math::Factor::XS from cpan. 12 13use Math::MatrixReal; 14use Math::Factor::XS qw(prime_factors); 15use Getopt::Std; 16$opt{n} = 1; 17getopts('hsdn:', \%opt); 18 19$usage .= "easyinverse.pl [-hsd] [-n N]\n"; 20$usage .= " -h this message\n"; 21$usage .= " -s symmetric\n"; 22$usage .= " -d diagonally dominant\n"; 23$usage .= " -n N number to generate (DEF 1)\n"; 24if ($opt{h}) { 25 print $usage; 26 exit; 27} 28 29sub randint{ # but not 0 30 my $min = shift or -5; 31 my $max = shift or +5; 32 33 my $ansa = 0; 34 while ($ansa == 0) { 35 my $r = rand; 36 $r *= ($max - $min + 1); 37 my $i = int($r); 38 $ansa = $min + $i; 39 } 40 return $ansa; 41} 42 43$size = 4; 44 45sub printmateq { 46 my $lbl = shift; 47 my $mat = shift; 48 my $round = shift; 49 50 if ($round) { 51 for my $r (1..$size) { 52 for my $c (1..$size) { 53 my $elt = $mat->element($r,$c); 54 my $pelt = abs($elt); 55 my $prnd = int($pelt * $round + 0.5) / $round; 56 my $rnd = $prnd * ($elt < 0 ? -1 : 1); 57 $mat->assign($r, $c, $rnd); 58 }} 59 } 60 61 my @rows; 62 for my $r (1..$size) { 63 my @row = map {$mat->element($r, $_)} (1..$size); 64 my $rowstr = '{' . (join ',', @row) . '}'; 65 push @rows, $rowstr; 66 } 67 my $matstr = '{' . (join ',', @rows) . '}'; 68 $matstr =~ s!\.?000+\d+!!g; 69 print "$lbl = $matstr;\n"; 70} 71 72 73 74for (1..$opt{n}) { # requested number of matrices 75 76 $M = new Math::MatrixReal($size,$size); 77 78 while (1) { 79 for $r (1..$size) { 80 if ($opt{d}) { $M->assign($r,$r, randint(5,10)) } 81 else { $M->assign($r,$r, randint(-3, 3)) } 82 83 for $c ($r+1..$size) { 84 $i = randint(-5,5); 85 $M->assign($r, $c, $i); 86 $M->assign($c, $r, ($opt{s} ? $i : randint(-5,5))); 87 } 88 } 89 90 $det = int($M->det + 0.5); 91 next if $det == 0; 92 @factors = prime_factors(abs($det)); 93 $all_good_factors = 1; # cross your fingers! 94 for $f (@factors) { 95 if ($f != 2 || $f != 5) { 96 $all_good_factors = 0; 97 last; 98 } 99 } 100 next unless $all_good_factors; 101 102 # test for positive semidefinite ~ all positive eigenvalues 103 $evals = $M->sym_eigenvalues; 104 $all_positive_evals = 1; # keep hope alive! 105 for $r (1..$size) { 106 if ($evals->element($r,1) <= 0) { 107 $all_positive_evals = 0; 108 last; 109 } 110 } 111 last if $all_positive_evals; 112 } 113 114 printmateq(".mat", $M); 115 printmateq(".inv", $M->inverse, 1000); 116 117} 118 119 120