Hi Monks.. I am trying to multiply matrices but i want someway to calculate the matrix to the power of x

If your matrices are square, you can decompose your matrix into a product of the form M = U * D * U^T where U is a unitary matrix, U^T is its transposition, and D is a diagonal matrix. Calculating M^x then becomes M^x = U * D^x * U^T, and since D is diagonal calculating the xth power is the same as taking each component to the xth power.

How does the magical decomposition work? My linear algebra fu is a bit rusty, but Matrix decomposition says the singular value decomposition (SVD) can help.

For many years, I too have suggested people look at PDL, but it's only now reflecting on zentara's suggestion that I see PDL could really use some example code up on perlmonks. Unfortunately, I'm not a big PDL user, so you might want to jump onto the PDL mailing list or peruse the Quantified Onion for more knowledgeable answers.

Taking a stab at an example, would the PDL code look something like this?

use PDL;

# amino acids in this order in the reaction rates pdl
# A R N D C Q E G H I L K M F P S T W Y V B Z X
$rates = pdl [
    [ 5, -2, -1 ],
    [ -2, 7, -1 ],
    [ -1, -1, 7 ],
];

$multiplied_rates = $rates x $rates;
print $multiplied_rates;
[download]

My calculation by hand seems to agree with the first two values. The rest is left as an exercise for the user.

Looks a lot easier than the OP's code, if you can overcome the fear of the new. And it's kinda cool :)

Edit:If I'd checked before posting, I would have linked to The Perl Data Language (PDL): A Quick Reference Guide or RFC: Getting Started with PDL (the Perl Data Language) both of which have lots of examples.

Sometimes I can think of 6 impossible LDAP attributes before breakfast.

I believe the perfect solution is in PDL::LinearAlgebra::Trans, specifically the mpow function: raise a matrix to a power.

@a = qw/2 2 2 5/;
#@a = qw/3 3 3 2 2 2 4 4 4/;
$n = 4;

powMatrix(\@a, $n);

sub powMatrix {
    $m = shift;
    @f = @$m;
    $i = shift;
    $l = sqrt(@$m);
    while(--$i){
        for($k = 0; $k < @$m; $k+=$l){        
            for($c = 0; $c < $l; $c++){
            $o = $c; $z = 0;
                for($b = $k; $b < $k+$l; $b++){
                $z += $f[$b] * $m->[$o];
                $o += $l;    
                }
            push (@g, $z);
            }
        }
    @f = @g; @g = ();
    }
    print "@f\n";
}
[download]