Math-MatrixReal

 view release on metacpan or  search on metacpan

lib/Math/MatrixReal.pm  view on Meta::CPAN

        my $m;
        OUTER:
        do {
            for ($m = $l; $m < ($rows - 1); $m++)
            {
                my $dd = abs($d->[$m]) + abs($d->[$m+1]);
                last if ((abs($e->[$m]) + $dd) == $dd);
            }
            if ($m != $l)
            {
                ## why only allow 30 iterations?
                croak("Too many iterations!") if ($iter++ >= 30);
                my $g = ($d->[$l+1] - $d->[$l])
                    / (2.0 * $e->[$l]);
                my $r = _pythag($g, 1.0);
                $g = $d->[$m] - $d->[$l]
                    + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
                my ($p,$s,$c) = (0.0, 1.0,1.0);
            for (my $i = ($m-1); $i >= $l; $i--)
            {
                my $ii = $i + 1;
                my $f = $s * $e->[$i];

lib/Math/MatrixReal.pm  view on Meta::CPAN

        my $m;
        OUTER:
        do {
            for ($m = $l; $m < ($rows - 1); $m++)
            {
                my $dd = abs($d->[$m]) + abs($d->[$m+1]);
                last if ((abs($e->[$m]) + $dd) == $dd);
            }
            if ($m != $l)
            {
                croak("Too many iterations!") if ($iter++ >= 30);
                my $g = ($d->[$l+1] - $d->[$l])
                    / (2.0 * $e->[$l]);
                my $r = _pythag($g, 1.0);
                $g = $d->[$m] - $d->[$l]
                    + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
                my ($p,$s,$c) = (0.0, 1.0,1.0);
                for (my $i = ($m-1); $i >= $l; $i--)
                {
                    my $ii = $i + 1;
                    my $f = $s * $e->[$i];

lib/Math/MatrixReal.pm  view on Meta::CPAN

the native C<inverse()> function. Here is a small benchmark:

    # $matrix1 is 15x15
    $det = $matrix1->det;
    timethese( 10,
        {'inverse' => sub { $matrix1->inverse(); },
          'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
        } );


    Benchmark: timing 10 iterations of LR, cofactor, inverse...
        inverse:  1 wallclock secs ( 0.56 usr +  0.00 sys =  0.56 CPU) @ 17.86/s (n=10)
    cofactor: 36 wallclock secs (36.62 usr +  0.01 sys = 36.63 CPU) @  0.27/s (n=10)

=item *

C<$adjoint = $matrix-E<gt>adjoint();>

The adjoint is just the transpose of the cofactor matrix. This method is 
just an alias for C< ~($matrix-E<gt>cofactor)>.

lib/Math/MatrixReal.pm  view on Meta::CPAN

equation system "C<A * x = b>" using an analytical algorithm like
the "decompose_LR()" and "solve_LR()" method pair.

In fact in some cases, due to the numerical properties (the "condition")
of the matrix "A", the numerical error of the obtained result can be
greater than by using an approximative (iterative) algorithm like one
of the three implemented here.

All three methods, GSM ("Global Step Method" or "Gesamtschrittverfahren"),
SSM ("Single Step Method" or "Einzelschrittverfahren") and RM ("Relaxation
Method" or "Relaxationsverfahren"), are fix-point iterations, that is, can
be described by an iteration function "C<x(t+1) = Phi( x(t) )>" which has
the property:

  Phi(x)  =  x    <==>    A * x  =  b

We can define "C<Phi(x)>" as follows:

  Phi(x)  :=  ( En - A ) * x  +  b

where "En" is a matrix of the same size as "A" ("n" rows and columns)

t/eigen_3x3.t  view on Meta::CPAN

  my $b = Math::MatrixReal->new($rows, 1);
  for (my $i = 1; $i <= $rows; $i++)
    {
      $b->assign($i, 1, rand());
    }
  # Normalize it
  my $l = $b->length();
  $b->multiply_scalar($b, (1.0 / $l));
  # Now do LR decomposition for linear system
  my $inv_it_LR = $inv_it->decompose_LR();
  # Check iterations
  my $iter = 0;
  my $delta;
  do {
    my ($dim, $b_base, $base) = $inv_it_LR->solve_LR($b);
    # Normalize
    my $l = $b_base->length();
    $b_base->multiply_scalar($b_base, (-1.0 / $l));
#    print "b_base=\n$b_base";
    $b->subtract($b_base,$b);
    $delta = $b->norm_one();