Merge pull request #114 from SciML/AP_Jacobian_recursion

fix ArrayPartition ZeroMatrix, i.e. stiff ODEs on recursive AP

Author: ChrisRackauckas

src/array_partition.jl

@@ -302,7 +302,7 @@ common_number(a, b) =
 
 ## Linear Algebra
 
-ArrayInterface.zeromatrix(A::ArrayPartition) = ArrayInterface.zeromatrix(reduce(vcat,vec.(A.x)))
+ArrayInterface.zeromatrix(A::ArrayPartition) = ArrayInterface.zeromatrix(Vector(A))
 
 LinearAlgebra.ldiv!(A::Factorization, b::ArrayPartition) = (x = ldiv!(A, Array(b)); copyto!(b, x))
 function LinearAlgebra.ldiv!(A::LU, b::ArrayPartition)

test/upstream.jl

@@ -28,3 +28,32 @@ F = ArrayPartition([0.0 0.0; 0.0 0.0],[0.0 0.0; 0.0 0.0])
 u0= ArrayPartition([0.1 1.2; 0.1 1.2], [0.1 1.2; 0.1 1.2])
 result = mymodel(F, u0)
 nlsolve(mymodel, u0)
+
+# Nested ArrayPartition solves
+
+dyn(u, p, t) = ArrayPartition(
+    ArrayPartition(zeros(1), [0.0]),
+    ArrayPartition(zeros(1), [0.0])
+)
+
+solve(
+    ODEProblem(
+        dyn,
+        ArrayPartition(
+            ArrayPartition(zeros(1), [-1.0]),
+            ArrayPartition(zeros(1), [0.75])
+        ),
+        (0.0, 1.0)
+    ),AutoTsit5(Rodas5())
+)
+
+@test_broken solve(
+    ODEProblem(
+        dyn,
+        ArrayPartition(
+            ArrayPartition(zeros(1), [-1.0]),
+            ArrayPartition(zeros(1), [0.75])
+        ),
+        (0.0, 1.0)
+    ),Rodas5()
+).retcode == :Success