add MatrixOps::tritosquare

Basic/MatrixOps/matrixops.pd

@@ -1352,6 +1352,28 @@ Ignores upper half of input.
 ',
 );
 
+pp_def("tritosquare",
+	Pars => 'a(m); [o]b(n,n)',
+	GenericTypes => [ppdefs_all],
+	RedoDimsCode => '
+          float n = sqrtf(0.25 + 2*$SIZE(m)) - 0.5;
+          $SIZE(n) = roundf(n);
+          if (fabsf($SIZE(n) - n) > 0.1) $CROAK("Non-triangular vector size=%"IND_FLAG, $SIZE(m));
+	',
+	Code => '
+            register PDL_Indx mna=0, nb=0;
+            loop(m) %{
+                  $b(n0 => mna, n1 => nb) = $a();
+                  mna++; if(mna > nb) {mna = 0; nb ++;}
+            %}
+	',
+	Doc => '=for ref
+
+Convert a triangular vector to lower-triangular square matrix storage.
+Does not touch upper half of output.
+',
+);
+
 pp_addpm({At=>'Bot'},<<'EOD');
 
 =head1 AUTHOR

Changes

@@ -1,6 +1,7 @@
 - test, document PDL::string, make more consistent (#459) - thanks @vadim-160102 for report
 - fix pre-2002 rangeb bug (#457) - thanks @jo-37 for report
 - add datachgd method
+- add MatrixOps::tritosquare
 
 2.085 2024-01-30
 - switch FFT code to use heap, not VLA (#436) - thanks @HaraldJoerg for report

t/matrixops.t

@@ -332,6 +332,13 @@ eval {squaretotri($x, zeroes(7))};
 like $@, qr/dim has size 7/;
 $y = squaretotri($x);
 is $y.'', "[0 3 4 6 7 8]", 'squaretotri with no output arg given';
+is tritosquare($y).'', '
+[
+ [0 0 0]
+ [3 4 0]
+ [6 7 8]
+]
+', 'tritosquare';
 $y = squaretotri(sequence(3,3,2));
 is $y.'', "
 [