diff --git a/Project.toml b/Project.toml
index 1d5f1077..9f83fc39 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FinEtools"
 uuid = "91bb5406-6c9a-523d-811d-0644c4229550"
 authors = ["Petr Krysl <pkrysl@ucsd.edu>"]
-version = "7.1.3"
+version = "7.1.4"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
diff --git a/src/AlgoBaseModule.jl b/src/AlgoBaseModule.jl
index 088bb97a..bc4cdc45 100644
--- a/src/AlgoBaseModule.jl
+++ b/src/AlgoBaseModule.jl
@@ -629,16 +629,30 @@ function vector_blocked(V, nfreedofs)
 end
 
 """
-    solve!(u::F, K::M, F::V) where {F<:AbstractField, M<:AbstractMatrix, V<:AbstractVector}
+    solve_blocked(A::M, b::VB, x::VX, nfreedofs::IT) where {M<:AbstractMatrix, VB<:AbstractVector, VX<:AbstractVector, IT<:Integer}
+
+Solve a blocked system of linear algebraic equations.
+
+b_f and x_d are known, x_f and b_d need to be computed.
+"""
+function solve_blocked(A::M, b::VB, x::VX, nfreedofs::IT) where {M<:AbstractMatrix, VB<:AbstractVector, VX<:AbstractVector, IT<:Integer}
+    A_b = matrix_blocked(A, nfreedofs)
+    b_b = vector_blocked(b, nfreedofs)
+    x_b = vector_blocked(x, nfreedofs)
+    x_f = A_b.ff \ (b_b.f - A_b.fd * x_b.d)
+    b_d = A_b.df * x_f + A_b.dd * x_b.d
+    return x_f, b_d
+end
+
+"""
+    solve_blocked!(u::AF, K::M, F::V) where {AF<:AbstractField, M<:AbstractMatrix, V<:AbstractVector}
 
 Solve a system of linear algebraic equations.
 """
-function solve!(u::U, K::M, F::V) where {U<:AbstractField, M<:AbstractMatrix, V<:AbstractVector}
-    K_ff, K_fd = matrix_blocked(K, nfreedofs(u), nfreedofs(u), (:ff, :fd))[(:ff, :fd)]
-    F_f = vector_blocked(F, nfreedofs(u))[:f]
-    U_d = gathersysvec(u, :d)
-    U_f =  K_ff\(F_f - K_fd * U_d)
-    scattersysvec!(u, U_f)
+function solve_blocked!(u::AF, K::M, F::V) where {AF<:AbstractField, M<:AbstractMatrix, V<:AbstractVector}
+    U = gathersysvec(u, :a)
+    x_f, b_d = solve_blocked(K, F, U, nfreedofs(u))
+    scattersysvec!(u, x_f)
     return u
 end
 
diff --git a/test/test_basics.jl b/test/test_basics.jl
index fad2f55a..8c704fdf 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -2527,7 +2527,7 @@ end
 
 module mblocked1
 using Test
-using FinEtools.AlgoBaseModule: solve!, matrix_blocked, vector_blocked
+using FinEtools.AlgoBaseModule: matrix_blocked, vector_blocked
 using LinearAlgebra
 
 function test(n, nfreedofs)
@@ -2575,7 +2575,7 @@ end
 
 module mblocked2
 using Test
-using FinEtools.AlgoBaseModule: solve!, vector_blocked
+using FinEtools.AlgoBaseModule: vector_blocked
 using LinearAlgebra
 
 function test(n, nfreedofs)
@@ -2617,3 +2617,74 @@ nothing
 end
 
 
+module mblocked3
+using Test
+using FinEtools.AlgoBaseModule: solve_blocked, solve_blocked!, vector_blocked, matrix_blocked
+using LinearAlgebra
+
+function test(n, nfreedofs)
+    A = rand(n, n) + I(n)
+    b = rand(n)
+    x = A \ b
+
+    b_b = vector_blocked(b, nfreedofs)
+    x_b = vector_blocked(x, nfreedofs)
+    A_b = matrix_blocked(A, nfreedofs)
+
+    x_f = A_b.ff \ (b_b.f - A_b.fd * x_b.d)
+    b_d = A_b.df * x_f + A_b.dd * x_b.d
+
+    @test norm(x - vcat(x_f, x_b.d)) / norm(x)  < 1.0e-9
+    @test norm(x_f - x_b.f) / norm(x_f)  < 1.0e-9
+    @test norm(b_b.d - b_d) / norm(b_d)  < 1.0e-9
+
+    x_f, b_d  = solve_blocked(A, b, x, nfreedofs)
+
+    @test norm(x - vcat(x_f, x_b.d)) / norm(x)  < 1.0e-9
+    @test norm(x_f - x_b.f) / norm(x_f)  < 1.0e-9
+    @test norm(b_b.d - b_d) / norm(b_d)  < 1.0e-9
+
+    true
+end
+test(10, 9)
+test(100, 9)
+test(1000, 999)
+test(200, 99)
+nothing
+end
+
+module mblocked4
+using Test
+using FinEtools
+using FinEtools.AlgoBaseModule: solve_blocked, solve_blocked!, vector_blocked, matrix_blocked
+using LinearAlgebra
+
+function test(n, nfreed)
+    A = rand(n, n) + 2*I(n)
+    b = rand(n)
+    x = A \ b
+
+    u = NodalField(x)
+    u.is_fixed[nfreed+1:end] .= true
+    numberdofs!(u)
+
+    solve_blocked!(u, A, b)
+
+    x_f = gathersysvec(u, :f)
+    x_b = vector_blocked(x, nfreed)
+
+    @test norm(x_f - x_b.f) / norm(x_f)  < 1.0e-9
+
+    true
+end
+test(10, 9)
+test(100, 9)
+test(1000, 999)
+test(200, 99)
+test(130, 99)
+test(130, 1)
+nothing
+end
+
+
+