MatrixOps::identity now preserves higher dims on multi-dim input

Basic/MatrixOps/matrixops.pd

@@ -130,21 +130,22 @@ pp_addpm(<<'EOD');
 
 Return an identity matrix of the specified size.  If you hand in a
 scalar, its value is the size of the identity matrix; if you hand in a
-dimensioned PDL, the 0th dimension is the size of the matrix.
+dimensioned PDL, the 0th dimension is the first two dimensions of the
+matrix, with higher dimensions preserved.
 
 =cut
 
 sub identity {
   my $n = shift;
-  my $out = ((UNIVERSAL::isa($n,'PDL')) ? 
-	  (  ($n->getndims > 0) ? 
-	     zeroes($n->dim(0),$n->dim(0)) : 
-	     zeroes($n->at(0),$n->at(0))
-	  ) :
-	  zeroes($n,$n)
-	  );
-  my $tmp; # work around perl -d "feature"
-  ($tmp = $out->diagonal(0,1))++;
+  my $out =
+    !UNIVERSAL::isa($n,'PDL') ? zeroes($n,$n) :
+    $n->getndims == 0 ? zeroes($n->at(0),$n->at(0)) :
+    undef;
+  if (!defined $out) {
+    my @dims = $n->dims;
+    $out = zeroes(@dims[0, 0, 2..$#dims]);
+  }
+  (my $tmp = $out->diagonal(0,1))++; # work around perl -d "feature"
   $out;
 }
 EOD
@@ -269,9 +270,7 @@ sub inv {
     barf("PDL::inv: got a singular matrix or LU decomposition\n");
   }
 
-  my $idenA = $x->zeros;
-  $idenA->diagonal(0,1) .= 1;
-  my $out = lu_backsub($lu,$perm,$par,$idenA)->transpose->sever;
+  my $out = lu_backsub($lu,$perm,$par,identity($x))->transpose->sever;
 
   return $out
     unless($x->is_inplace);

Changes

@@ -1,5 +1,6 @@
 - PDL::Complex->from_native / as_native
 - zap obsolete MatrixOps::eigen_c
+- MatrixOps::identity now preserves higher dims on multi-dim input
 
 2.048 2021-05-24
 - remove obsolete and inaccurate doc relating bad values to NaN

t/matrixops.t

@@ -145,6 +145,10 @@ ok(all($det == pdl([48,1],[-1,-216])), "threaded determinant");
 {
 ### Check identity and stretcher matrices...
 ok((identity(2)->flat == pdl(1,0,0,1))->all, "identity matrix");
+ok((identity(pdl 2)->flat == pdl(1,0,0,1))->all, "identity matrix with scalar ndarray");
+ok((identity(zeroes 2, 3)->flat == pdl(1,0,0,1))->all, "identity matrix with dimensioned ndarray");
+my @deep_identity_dims = identity(zeroes 2, 3, 4)->dims;
+is_deeply \@deep_identity_dims, [2, 2, 4], "identity matrix with multi-dimensioned ndarray" or diag 'got: ', explain \@deep_identity_dims;
 ok((stretcher(pdl(2,3))->flat == pdl(2,0,0,3))->all, "stretcher 2x2");
 ok((stretcher(pdl([2,3],[3,4]))->flat == pdl(2,0,0,3,3,0,0,4))->all, "stretcher 2x2x2");
 }