Merge pull request #106 from SciML/ap_linalg_arrays

Handle and test ArrayPartition linear algebra over Arrays

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "RecursiveArrayTools"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "2.4.2"
+version = "2.4.3"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

src/array_partition.jl

@@ -331,11 +331,11 @@ function LinearAlgebra.ldiv!(A::T, bb::ArrayPartition) where T<:Union{UnitUpperT
     lens = map(length, b)
     @inbounds for j in n:-1:1
         Ajj = T(getblock(A, lens, j, j))
-        xj = ldiv!(Ajj, b[j])
+        xj = ldiv!(Ajj, vec(b[j]))
         for i in j-1:-1:1
             Aij = getblock(A, lens, i, j)
             # bi = -Aij * xj + bi
-            mul!(b[i], Aij, xj, -1, true)
+            mul!(vec(b[i]), Aij, xj, -1, true)
         end
     end
     return bb
@@ -348,11 +348,11 @@ function LinearAlgebra.ldiv!(A::T, bb::ArrayPartition) where T<:Union{UnitLowerT
     lens = map(length, b)
     @inbounds for j in 1:n
         Ajj = T(getblock(A, lens, j, j))
-        xj = ldiv!(Ajj, b[j])
+        xj = ldiv!(Ajj, vec(b[j]))
         for i in j+1:n
             Aij = getblock(A, lens, i, j)
             # bi = -Aij * xj + b[i]
-            mul!(b[i], Aij, xj, -1, true)
+            mul!(vec(b[i]), Aij, xj, -1, true)
         end
     end
     return bb

test/upstream.jl

@@ -14,19 +14,17 @@ sol = solve(prob,AutoTsit5(Rosenbrock23()))
 
 @test all(Array(sol) .== sol)
 
-function f!(F, vars)
-    x = vars.x[1]
-    F.x[1][1] = (x[1]+3)*(x[2]^3-7)+18
-    F.x[1][2] = sin(x[2]*exp(x[1])-1)
-    y=vars.x[2]
-    F.x[2][1] = (y[1]+3)*(y[2]^3-7)+18
-    F.x[2][2] = sin(y[2]*exp(y[1])-1)
+function mymodel(F, vars)
+    for i in 1:2
+        x = vars.x[i]
+        F.x[i][1,1] = (x[1,1]+3)*(x[1,2]^3-7)+18.0
+        F.x[i][1,2] = sin(x[1,2]*exp(x[1,1])-1)
+        F.x[i][2,1] = (x[2,1]+3)*(x[2,2]^3-7)+19.0
+        F.x[i][2,2] = sin(x[2,2]*exp(x[2,1])-3)
+    end
 end
-
 # To show that the function works
-F = ArrayPartition([0.0 0.0],[0.0, 0.0])
-u0= ArrayPartition([0.1; 1.2], [0.1; 1.2])
-result = f!(F, u0)
-
-# To show the NLsolve error that results with ArrayPartitions:
-nlsolve(f!, ArrayPartition([0.1; 1.2], [0.1; 1.2]))
+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)