Fix sf.net bug #256 again

The identity matrix being passed to lu_backsub() by inv() was not maintaining the extra thread dimensions. Added a test to t/matrixops.t to verify the fix.
PDLPorters · Mar 5, 2015 · 538f34a · 538f34a
1 parent aedd1ea
commit 538f34a
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 3 deletions.
diff --git a/Basic/MatrixOps/matrixops.pd b/Basic/MatrixOps/matrixops.pd
@@ -299,7 +299,9 @@ sub inv {
     barf("PDL::inv: got a singular matrix or LU decomposition\n");
   }
 
-  my $out = lu_backsub($lu,$perm,$par,identity($a))->xchg(0,1)->sever;
+  my $idenA = $a->zeros;
+  $idenA->diagonal(0,1) .= 1;
+  my $out = lu_backsub($lu,$perm,$par,$idenA)->xchg(0,1)->sever;
 
   return $out
     unless($a->is_inplace);

diff --git a/t/matrixops.t b/t/matrixops.t
@@ -16,7 +16,7 @@ sub near {
 	return ($dist <= $tol)->all;
 }
 
-BEGIN { plan tests => 30,
+BEGIN { plan tests => 34,
 }
 
 my $tol = 1e-14;
@@ -64,6 +64,15 @@ ok(defined $a1);
 ok(ref ($opt->{lu}->[0]) eq 'PDL');
 ok(near(matmult($a1,$a),$identity,$tol));
 
+### Check inv() with added thread dims (simple check)
+my $C22 = pdl([5,5],[5,7.5]);
+my $C22inv = eval { $C22->inv };
+ok(!$@);                                                    # ran OK
+ok(near($C22inv,pdl([0.6, -0.4], [-0.4, 0.4])));            # right answer
+my $C222 = $C22->dummy(2,2);
+my $C222inv = eval { $C222->inv };
+ok(!$@);                                                         # ran OK
+ok(near($C222inv,pdl([0.6, -0.4], [-0.4, 0.4])->dummy(2,2)));    # right answer
 
 ### Check inv() for matrices with added thread dims (bug #3172882 on sf.net)
 $a94 = pdl( [  1,  0,  4, -1, -1, -3,  0,  1,  0 ],
@@ -134,7 +143,7 @@ ok(!$@);
 
 ### Check that it really returns eigenvectors
 $c = float(($a x $vec) / $vec); 
-print "c is $c\n";
+#print "c is $c\n";
 ok(all($c->slice(":,0") == $c->slice(":,1")));
 
 ### Check that the eigenvalues are correct for this matrix