diff --git a/Manifest.toml b/Manifest.toml
index 043e57b..fc76e61 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -6,23 +6,35 @@ uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
 version = "0.0.1"
 
 [[AbstractFFTs]]
-deps = ["ChainRulesCore", "LinearAlgebra"]
-git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "1.1.0"
+version = "1.5.0"
+
+[[AbstractTrees]]
+git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+version = "0.4.4"
 
 [[Adapt]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.3.3"
+version = "3.7.2"
 
 [[ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
 
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
+[[Atomix]]
+deps = ["UnsafeAtomics"]
+git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be"
+uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+version = "0.1.0"
+
 [[BDIO]]
 deps = ["Documenter", "Nettle", "Test"]
 git-tree-sha1 = "45f43efe91dcda1939cfa5b0b028ea941afdc367"
@@ -33,45 +45,97 @@ version = "0.1.0"
 
 [[BFloat16s]]
 deps = ["LinearAlgebra", "Printf", "Random", "Test"]
-git-tree-sha1 = "a598ecb0d717092b5539dbbe890c98bac842b072"
+git-tree-sha1 = "dbf84058d0a8cbbadee18d25cf606934b22d7c66"
 uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
-version = "0.2.0"
+version = "0.4.2"
 
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[CEnum]]
-git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
+git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc"
 uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.4.2"
+version = "0.5.0"
 
 [[CUDA]]
-deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "TimerOutputs"]
-git-tree-sha1 = "429a1a05348ce948a96adbdd873fbe6d9e5e052f"
+deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "Statistics", "UnsafeAtomicsLLVM"]
+git-tree-sha1 = "76582ae19006b1186e87dadd781747f76cead72c"
 uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
-version = "3.6.2"
+version = "5.1.1"
+
+[[CUDA_Driver_jll]]
+deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "1e42ef1bdb45487ff28de16182c0df4920181dc3"
+uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
+version = "0.7.0+0"
+
+[[CUDA_Runtime_Discovery]]
+deps = ["Libdl"]
+git-tree-sha1 = "bcc4a23cbbd99c8535a5318455dcf0f2546ec536"
+uuid = "1af6417a-86b4-443c-805f-a4643ffb695f"
+version = "0.2.2"
+
+[[CUDA_Runtime_jll]]
+deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
+git-tree-sha1 = "9704e50c9158cf8896c2776b8dbc5edd136caf80"
+uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
+version = "0.10.1+0"
 
 [[ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "9489214b993cd42d17f44c36e359bf6a7c919abf"
+git-tree-sha1 = "2118cb2765f8197b08e5958cdd17c165427425ee"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.15.0"
+version = "1.19.0"
 
-[[ChangesOfVariables]]
-deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
-git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
-uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.3"
+[[ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.11.4"
+
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.10"
 
 [[Compat]]
 deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "924cdca592bc16f14d2f7006754a621735280b74"
+git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.1.0"
+version = "4.10.1"
 
 [[CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.1+0"
+
+[[Crayons]]
+git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
+uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+version = "4.1.1"
+
+[[DataAPI]]
+git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.15.0"
+
+[[DataFrames]]
+deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
+git-tree-sha1 = "04c738083f29f86e62c8afc341f0967d8717bdb8"
+uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+version = "1.6.1"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.15"
+
+[[DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
 
 [[Dates]]
 deps = ["Printf"]
@@ -79,88 +143,156 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[DocStringExtensions]]
 deps = ["LibGit2"]
-git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.6"
+version = "0.9.3"
 
 [[Documenter]]
-deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "f7809f532671564e48cd81627ddcfb1ba670f87d"
+deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"]
+git-tree-sha1 = "2613dbec8f4748273bbe30ba71fd5cb369966bac"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.19"
+version = "1.2.1"
 
 [[Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[Expat_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c"
+uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
+version = "2.5.0+0"
 
 [[ExprTools]]
-git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
+git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
-version = "0.1.8"
+version = "0.1.10"
 
 [[FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
+[[FixedPointNumbers]]
+deps = ["Statistics"]
+git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.8.4"
+
+[[Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
 [[GMP_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
+version = "6.2.1+2"
 
 [[GPUArrays]]
-deps = ["Adapt", "LLVM", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"]
-git-tree-sha1 = "c783e8883028bf26fb05ed4022c450ef44edd875"
+deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
+git-tree-sha1 = "85d7fb51afb3def5dcb85ad31c3707795c8bccc1"
 uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
-version = "8.3.2"
+version = "9.1.0"
+
+[[GPUArraysCore]]
+deps = ["Adapt"]
+git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
+uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
+version = "0.1.5"
 
 [[GPUCompiler]]
-deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "TimerOutputs", "UUIDs"]
-git-tree-sha1 = "647a54f196b5ffb7c3bc2fec5c9a57fa273354cc"
+deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"]
+git-tree-sha1 = "a846f297ce9d09ccba02ead0cae70690e072a119"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
-version = "0.13.14"
+version = "0.25.0"
+
+[[Git]]
+deps = ["Git_jll"]
+git-tree-sha1 = "51764e6c2e84c37055e846c516e9015b4a291c7d"
+uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
+version = "1.3.0"
+
+[[Git_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
+git-tree-sha1 = "bb8f7cc77ec1152414b2af6db533d9471cfbb2d1"
+uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
+version = "2.42.0+0"
 
 [[IOCapture]]
 deps = ["Logging", "Random"]
-git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
+git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6"
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
-version = "0.2.2"
+version = "0.2.3"
+
+[[InlineStrings]]
+deps = ["Parsers"]
+git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461"
+uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
+version = "1.4.0"
 
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
-[[InverseFunctions]]
-deps = ["Test"]
-git-tree-sha1 = "b3364212fb5d870f724876ffcd34dd8ec6d98918"
-uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.7"
+[[InvertedIndices]]
+git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
+uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
+version = "1.3.0"
 
-[[IrrationalConstants]]
-git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
-uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
-version = "0.1.1"
+[[IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
 
 [[JLLWrappers]]
-deps = ["Preferences"]
-git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+deps = ["Artifacts", "Preferences"]
+git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.1"
+version = "1.5.0"
 
 [[JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
+git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.3"
+version = "0.21.4"
+
+[[JuliaNVTXCallbacks_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "af433a10f3942e882d3c671aacb203e006a5808f"
+uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e"
+version = "0.2.1+0"
+
+[[KernelAbstractions]]
+deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
+git-tree-sha1 = "653e0824fc9ab55b3beec67a6dbbe514a65fb954"
+uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+version = "0.9.15"
 
 [[LLVM]]
-deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"]
-git-tree-sha1 = "e7e9184b0bf0158ac4e4aa9daf00041b5909bf1a"
+deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
+git-tree-sha1 = "0678579657515e88b6632a3a482d39adcbb80445"
 uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "4.14.0"
+version = "6.4.1"
 
 [[LLVMExtra_jll]]
-deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg", "TOML"]
-git-tree-sha1 = "771bfe376249626d3ca12bcd58ba243d3f961576"
+deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
+git-tree-sha1 = "98eaee04d96d973e79c25d49167668c5c8fb50e2"
 uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
-version = "0.0.16+0"
+version = "0.0.27+1"
+
+[[LLVMLoopInfo]]
+git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea"
+uuid = "8b046642-f1f6-4319-8d3c-209ddc03c586"
+version = "1.0.0"
+
+[[LaTeXStrings]]
+git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
+uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+version = "1.3.1"
+
+[[LazilyInitializedFields]]
+git-tree-sha1 = "8f7f3cabab0fd1800699663533b6d5cb3fc0e612"
+uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf"
+version = "1.2.2"
 
 [[LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
@@ -169,10 +301,12 @@ uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.3"
 
 [[LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "7.84.0+0"
 
 [[LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
@@ -181,42 +315,75 @@ uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 [[LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.10.2+0"
 
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[[Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.17.0+0"
+
 [[LinearAlgebra]]
 deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[LogExpFunctions]]
-deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
-uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.15"
-
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.11"
+
 [[Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[MarkdownAST]]
+deps = ["AbstractTrees", "Markdown"]
+git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899"
+uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391"
+version = "0.1.2"
+
 [[MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.0+0"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "1.1.0"
 
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2022.2.1"
+
+[[NVTX]]
+deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"]
+git-tree-sha1 = "8bc9ce4233be3c63f8dcd78ccaf1b63a9c0baa34"
+uuid = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
+version = "0.3.3"
+
+[[NVTX_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "ce3269ed42816bf18d500c9f63418d4b0d9f5a3b"
+uuid = "e98f9f5b-d649-5603-91fd-7774390e6439"
+version = "3.1.0+2"
 
 [[Nettle]]
 deps = ["Libdl", "Nettle_jll"]
-git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
+git-tree-sha1 = "6fa48cbae828267848ee32c1bb31d1652e210d7d"
 uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
-version = "0.5.1"
+version = "1.0.0"
 
 [[Nettle_jll]]
 deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -226,36 +393,63 @@ version = "3.7.2+0"
 
 [[NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
 
 [[OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.20+0"
 
-[[OpenLibm_jll]]
+[[OpenSSL_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f"
+uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
+version = "3.0.12+0"
+
+[[OrderedCollections]]
+git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.6.3"
+
+[[PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
-uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-
-[[OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
-uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.5+0"
+uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
+version = "10.40.0+0"
 
 [[Parsers]]
-deps = ["Dates"]
-git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
+deps = ["Dates", "PrecompileTools", "UUIDs"]
+git-tree-sha1 = "a935806434c9d4c506ba941871b327b96d41f2bf"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.3.1"
+version = "2.8.0"
 
 [[Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.8.0"
+
+[[PooledArrays]]
+deps = ["DataAPI", "Future"]
+git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3"
+uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
+version = "1.4.3"
+
+[[PrecompileTools]]
+deps = ["Preferences"]
+git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+version = "1.2.0"
 
 [[Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.3.0"
+version = "1.4.1"
+
+[[PrettyTables]]
+deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"]
+git-tree-sha1 = "88b895d13d53b5577fd53379d913b9ab9ac82660"
+uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+version = "2.3.1"
 
 [[Printf]]
 deps = ["Unicode"]
@@ -271,9 +465,9 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[Random123]]
 deps = ["Random", "RandomNumbers"]
-git-tree-sha1 = "afeacaecf4ed1649555a19cb2cad3c141bbc9474"
+git-tree-sha1 = "552f30e847641591ba3f39fd1bed559b9deb0ef3"
 uuid = "74087812-796a-5b5d-8853-05524746bad3"
-version = "1.5.0"
+version = "1.6.1"
 
 [[RandomNumbers]]
 deps = ["Random", "Requires"]
@@ -286,6 +480,12 @@ git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
 version = "1.2.2"
 
+[[RegistryInstances]]
+deps = ["LazilyInitializedFields", "Pkg", "TOML", "Tar"]
+git-tree-sha1 = "ffd19052caf598b8653b99404058fce14828be51"
+uuid = "2792f1a3-b283-48e8-9a74-f99dce5104f3"
+version = "0.1.0"
+
 [[Requires]]
 deps = ["UUIDs"]
 git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
@@ -294,6 +494,19 @@ version = "1.3.0"
 
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
+
+[[Scratch]]
+deps = ["Dates"]
+git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
+uuid = "6c6a2e73-6563-6170-7368-637461726353"
+version = "1.2.1"
+
+[[SentinelArrays]]
+deps = ["Dates", "Random"]
+git-tree-sha1 = "0e7508ff27ba32f26cd459474ca2ede1bc10991f"
+uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
+version = "1.4.1"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -301,27 +514,58 @@ uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "5165dfb9fd131cf0c6957a3a7605dede376e7b63"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.2.0"
+
 [[SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[[SpecialFunctions]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "a9e798cae4867e3a41cae2dd9eb60c047f1212db"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.1.6"
+[[StaticArrays]]
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "fba11dbe2562eecdfcac49a05246af09ee64d055"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.8.1"
+
+[[StaticArraysCore]]
+git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.2"
 
 [[Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[[StringManipulation]]
+deps = ["PrecompileTools"]
+git-tree-sha1 = "a04cabe79c5f01f4d723cc6704070ada0b9d46d5"
+uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
+version = "0.3.4"
+
 [[TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
+
+[[TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits"]
+git-tree-sha1 = "cb76cf677714c095e535e3501ac7954732aeea2d"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.11.1"
 
 [[Tar]]
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.1"
 
 [[Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
@@ -340,18 +584,33 @@ uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
+[[UnsafeAtomics]]
+git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278"
+uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
+version = "0.2.1"
+
+[[UnsafeAtomicsLLVM]]
+deps = ["LLVM", "UnsafeAtomics"]
+git-tree-sha1 = "323e3d0acf5e78a56dfae7bd8928c989b4f3083e"
+uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249"
+version = "0.1.3"
+
 [[Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.12+3"
 
 [[libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.1.1+0"
 
 [[nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.48.0+0"
 
 [[p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+0"
diff --git a/docs/src/dirac.md b/docs/src/dirac.md
index 2f7fade..d0ecc3a 100644
--- a/docs/src/dirac.md
+++ b/docs/src/dirac.md
@@ -20,6 +20,16 @@ The workspace stores four fermion fields, namely `.sr`, `.sp`, `.sAp` and `.st`,
 for different purposes. If the representation is either `SU2fund` of `SU3fund`, an extra
 field with values in `U2alg`/`U3alg` is created to store the clover, used for the improvement.
 
+The functions using the fields allocated in [`DiracWorkspace`](@ref) are the following:
+
+- `dws.sr` : [`CG!`](@ref), [`flw_adapt`](@ref) (fermion case), [`bfl_error`](@ref)
+- `dws.st` : [`DwdagDw!`](@ref), [`bflw_step_vec!`](@ref)
+- `dws.sp` : [`CG!`](@ref), [`flw`](@ref) (fermion case), [`bflw_step!`](@ref), [`bflw_step_vec!`](@ref), [`propagator!`](@ref), [`bndpropagator!`](@ref), [`Tbndpropagator!`](@ref)
+- `dws.sAp` : [`CG!`](@ref), [`flw`](@ref) (fermion case), [`bflw_step!`](@ref), [`bflw_step_vec!`](@ref) 
+
+Note that other functions may call some of these functions, like [`flw_adapt`](@ref) depending on [`flw`](@ref), [`bflw!`](@ref) depending on [`bflw_step!`](@ref) or [`propagator!`](@ref) depending on [`CG!`](@ref). The fields used in the innermost function will also be modified by the outermost methods.
+
+
 ## Functions
 
 The functions [`Dw!`](@ref), [`g5Dw!`](@ref) and [`DwdagDw!`](@ref) are all related to the 
diff --git a/src/Dirac/Diracfields.jl b/src/Dirac/Diracfields.jl
index 5559ca3..9274d11 100644
--- a/src/Dirac/Diracfields.jl
+++ b/src/Dirac/Diracfields.jl
@@ -151,16 +151,17 @@ function krnl_assign_pf_su3!(f::AbstractArray, p , lp::SpaceParm, t::Int64)
         b = Int64(CUDA.threadIdx().x)
         r = Int64(CUDA.blockIdx().x)
 
-            if t == 0
+        if t == 0
             f[b,r] = Spinor(map(x->SU3fund(x[b,1,r,1] + im* x[b,1,r,2],
-                                        x[b,2,r,1] + im* x[b,2,r,2],
-                                        x[b,3,r,1] + im* x[b,3,r,2]),p))
-            elseif point_time((b,r),lp) == t
+                                           x[b,2,r,1] + im* x[b,2,r,2],
+                                           x[b,3,r,1] + im* x[b,3,r,2]),p))
+        elseif point_time((b,r),lp) == t
             f[b,r] = Spinor(map(x->SU3fund(x[b,1,r,1] + im* x[b,1,r,2],
-                                        x[b,2,r,1] + im* x[b,2,r,2],
-                                        x[b,3,r,1] + im* x[b,3,r,2]),p))
-            end
-
+                                           x[b,2,r,1] + im* x[b,2,r,2],
+                                           x[b,3,r,1] + im* x[b,3,r,2]),p))
+        else
+            f[b,r] = 0.0*f[b,r]
+        end
     end
 
     return nothing
@@ -197,14 +198,15 @@ function krnl_assign_pf_su2!(f::AbstractArray, p , lp::SpaceParm, t::Int64)
         b = Int64(CUDA.threadIdx().x)
         r = Int64(CUDA.blockIdx().x)
 
-            if t == 0
+        if t == 0
             f[b,r] = Spinor(map(x->SU2fund(x[b,1,r,1] + im* x[b,1,r,2],
-                                        x[b,2,r,1] + im* x[b,2,r,2]),p))
-            elseif point_time((b,r),lp) == t
+                                           x[b,2,r,1] + im* x[b,2,r,2]),p))
+        elseif point_time((b,r),lp) == t
             f[b,r] = Spinor(map(x->SU2fund(x[b,1,r,1] + im* x[b,1,r,2],
-                                        x[b,2,r,1] + im* x[b,2,r,2]),p))
-            end
-
+                                           x[b,2,r,1] + im* x[b,2,r,2]),p))
+        else
+            f[b,r] = 0.0*f[b,r]
+        end
     end
 
     return nothing
diff --git a/src/Dirac/Diracflow.jl b/src/Dirac/Diracflow.jl
index 8065840..ee4488e 100644
--- a/src/Dirac/Diracflow.jl
+++ b/src/Dirac/Diracflow.jl
@@ -41,13 +41,13 @@ flw(U, psi, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, dpar::DiracParam, lp:
 """
     function backflow(psi, U, Dt, nsave::Int64, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace)
 
-Performs one step back in flow time for the fermion field, according to 1302.5246. The fermion field must me that of the time-slice Dt and is flowed back to the first time-slice
+Performs the integration of the adjoint flow for the fermion field, according to 1302.5246. The fermion field must me that of the time-slice Dt and is flowed back to the first time-slice
 nsave is the total number of gauge fields saved in the process
 
 """
-function backflow(psi, U, Dt, maxnsave::Int64, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace)
+function backflow(psi, U, Dt, maxnsave::Int64, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm,int::FlowIntr, ymws::YMworkspace, dws::DiracWorkspace)
 
-    int = wfl_rk3(Float64,0.01,1.0) # Default integrator, it has to be order 3 rk but in can be zfl
+     # Default integrator is  wfl_rk3(Float64,0.01,1.0), it has to be order 3 rk but in can be zfl
 
     @timeit "Backflow integration" begin
         @timeit "GPU to CPU" U0 = Array(U)
@@ -98,6 +98,7 @@ function backflow(psi, U, Dt, maxnsave::Int64, gp::GaugeParm, dpar::DiracParam,
 
     return nothing
 end
+backflow(psi, U, Dt, maxnsave::Int64, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace) = backflow(psi, U, Dt, maxnsave, gp, dpar, lp, wfl_rk3(Float64,0.01,1.0), ymws, dws)
 
 """
 function bflw_step!(U, psi, eps, int::FlowIntr, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace)
@@ -108,8 +109,7 @@ function bflw_step!(psi, U,  eps, int::FlowIntr, gp::GaugeParm, dpar::DiracParam
 
     @timeit "Backflow step" begin
 
-        V = copy(U)
-        V .= U
+        @timeit "GPU to CPU" V = Array(U)
 
         force_gauge(ymws, U, int.c0, 1, gp, lp)
 
@@ -131,7 +131,7 @@ function bflw_step!(psi, U,  eps, int::FlowIntr, gp::GaugeParm, dpar::DiracParam
 
         Nablanabla!(dws.sp, U, 0.75*2*eps*psi, dpar, dws, lp)
 
-        U .= V
+        @timeit "CPU to GPU" copyto!(U,V)
 
         force_gauge(ymws, U, int.c0, 1, gp, lp)
 
@@ -144,7 +144,7 @@ function bflw_step!(psi, U,  eps, int::FlowIntr, gp::GaugeParm, dpar::DiracParam
         Nablanabla!(dws.sAp, U, 2*eps*dws.sp, dpar, dws, lp)
         dws.sAp .= psi + (8/9)*dws.sAp
 
-        U .= V
+        @timeit "CPU to GPU" copyto!(U,V)
 
         Nablanabla!(psi, U, 2*eps*(dws.sAp - (8/9)*dws.sp), dpar, dws, lp)
         psi .= (1/4)*psi + dws.sp + dws.sAp
@@ -166,8 +166,9 @@ function flw_adapt(U, psi, int::FlowIntr{NI,T}, tend::T, epsini::T, gp::GaugePar
         if ns > 10
             flw(U, psi, int, 9, eps, gp, dpar, lp, ymws, dws)
             ymws.U1 .= U
+            dws.sr .= psi
             flw(U, psi, int, 1, eps, gp, dpar, lp, ymws, dws)
-            flw(ymws.U1, int, 2, eps/2, gp, lp, ymws)
+            flw(ymws.U1,dws.sr, int, 2, eps/2, gp, dpar,lp, ymws,dws)
 
             dt = dt - 10*eps
             nstp = nstp + 10
@@ -175,8 +176,10 @@ function flw_adapt(U, psi, int::FlowIntr{NI,T}, tend::T, epsini::T, gp::GaugePar
 
             # adjust step size
             ymws.U1 .= ymws.U1 ./ U
+            dws.sr .= dws.sr .- psi
             maxd = CUDA.mapreduce(dev_one, max, ymws.U1, init=zero(tend))
-            eps  = min(int.max_eps, 2*eps, int.sft_fac*eps*(int.tol/maxd)^(one(tend)/3))
+            pfdist = sqrt(CUDA.mapreduce(norm2, max, dws.sr, init=zero(tend)))
+            eps  = min(int.max_eps, 2*eps, int.sft_fac*eps*(int.tol/maxd)^(one(tend)/3),int.sft_fac*eps*(int.tol/pfdist)^(one(tend)/3))
 
         else
             flw(U, psi, int, ns, eps, gp, dpar, lp, ymws, dws)
@@ -205,7 +208,7 @@ flw_adapt(U, psi, int::FlowIntr{NI,T}, tend::T, gp::GaugeParm, dpar::DiracParam,
 
     function Nablanabla!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D})
 
-Computes /`/` \\nabla^* \\nabla /`/` `si` and stores it in `si`.
+Computes /`/` \\nabla^* \\nabla /`/` `si` and stores it in `so`.
 
 """
 function Nablanabla!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,BC_PERIODIC,D}) where {D}
@@ -216,6 +219,7 @@ function Nablanabla!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Space
     end
     return nothing
 end
+
 function Nablanabla!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D},SpaceParm{4,6,BC_OPEN,D}}) where {D}
     SF_bndfix!(si,lp)
     @timeit "Laplacian" begin
@@ -238,7 +242,7 @@ function krnl_Nablanabla(so, U, si, th, lp::SpaceParm{4,6,BC_OPEN,D}) where {D}
 
             so[b,r] = -4*si[b,r]
 
-	        bu1, ru1 = up((b,r), 1, lp)
+            bu1, ru1 = up((b,r), 1, lp)
             bd1, rd1 = dw((b,r), 1, lp)
             bu2, ru2 = up((b,r), 2, lp)
             bd2, rd2 = dw((b,r), 2, lp)
@@ -313,9 +317,29 @@ function krnl_Nablanabla(so, U, si, th, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},Sp
 end
 
 
-
 export Nablanabla!, flw, backflow, flw_adapt, bflw_step!
 
+"""
+    function bfl_error(psi_t, psi_0, U, tend, int::FlowIntr, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace)
+
+Estimates the error of the backflow integration of `\\psi\\_t` into `\\psi\\_0` with a random noise source.
+"""
+function bfl_error(psi_t, psi_0, U, tend, int::FlowIntr, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace, dws::DiracWorkspace)
+
+    pfrandomize!(dws.sr,lp)
+    @timeit "GPU to CPU" V = Array(U)
+
+    R0 =  sum(dot.(psi_0,dws.sr))
+
+    flw_adapt(U, dws.sr, int, tend, int.eps_ini/2, gp, dpar, lp, ymws, dws)
+
+    R1 = sum(dot.(psi_t,dws.sr))
+    @timeit "CPU to GPU" copyto!(U,V)
+
+    return abs(R0-R1)
+end
+
+export bfl_error
 
 """
     function Dslash_sq!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D})
@@ -362,7 +386,6 @@ function Dslash_sq!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceP
     return nothing
 end
 
-
 function krnl_g5Dslsh!(so, U, si, th, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
@@ -393,7 +416,6 @@ function krnl_g5Dslsh!(so, U, si, th, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},Spac
     return nothing
 end
 
-
 function krnl_g5Dslsh!(so, U, si, th, lp::SpaceParm{4,6,B,D}) where {D,B}
 
     b = Int64(CUDA.threadIdx().x);  r = Int64(CUDA.blockIdx().x)
@@ -436,8 +458,6 @@ function krnl_g5Dslsh_impr!(so, Fcsw, csw, si, lp::SpaceParm{4,6,B,D}) where {B,
     return nothing
 end
 
-
-
 function krnl_g5Dslsh_impr!(so, Fcsw, csw, si, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D}
 
     @inbounds begin
diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index 46f5ac6..577eb11 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -40,7 +40,7 @@ include("YM/YM.jl")
 using .YM
 export ztwist
 export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2
-export force_gauge, MD!
+export force_gauge, force_gauge_flw, MD!
 export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
 export Eoft_clover, Eoft_plaq, Qtop
 export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
@@ -60,7 +60,7 @@ using .Dirac
 export DiracWorkspace, DiracParam
 export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!, mtwmdpar
 export read_prop, save_prop, read_dpar
-export Nablanabla!, flw, backflow
+export Nablanabla!, flw, backflow, bfl_error
 
 include("Solvers/Solvers.jl")
 using .Solvers
diff --git a/src/YM/YM.jl b/src/YM/YM.jl
index 32ddb71..5578730 100644
--- a/src/YM/YM.jl
+++ b/src/YM/YM.jl
@@ -165,7 +165,7 @@ include("YMfields.jl")
 export randomize!, zero!, norm2
 
 include("YMact.jl")
-export krnl_plaq!, force_gauge, force_wilson
+export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson
 
 include("YMhmc.jl")
 export gauge_action, hamiltonian, plaquette, HMC!, MD!
diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl
index 4c03097..272fe91 100644
--- a/src/YM/YMact.jl
+++ b/src/YM/YMact.jl
@@ -320,6 +320,22 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG,
     return nothing
 end
 
+function bnd_rescale_flw!(frc1, lp::SpaceParm{N,M,BC_OPEN,D}) where {N,M,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+        I = point_coord((b,r), lp)
+        it = I[N]
+
+        for id in 1:N-1
+            if (((it == 1) || (it == lp.iL[4])))
+                frc1[b,id,r] = 2*frc1[b,id,r]
+            end
+        end
+    end
+    return nothing
+end
 
 ##
 ## SF
@@ -874,7 +890,6 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG,
 end
 
 
-
 ##
 ## PERIODIC
 ##
@@ -1143,6 +1158,38 @@ end
 force_gauge(ymws::YMworkspace, U, c0, gp, lp) = force_gauge(ymws, U, c0, gp.cG[1], gp, lp)
 force_gauge(ymws::YMworkspace, U, gp, lp) = force_gauge(ymws, U, gp.c0, gp.cG[1], gp, lp)
 
+"""
+    function force_gauge_flw(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D})
+
+Computes the force for the gauge flow with Open Boundaries. An aditional factor two in the boundaries
+is included, see
+
+M. Luescher, S. Schaefer: "Lattice QCD with open boundary conditions and twisted-mass reweighting", Comput.Phys.Commun. 184 (2013) 519,
+
+for more details.
+
+"""
+function force_gauge_flw(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}) where {NI,N,M,D}
+
+    ztw = ztwist(gp, lp)
+    if abs(c0-1) < 1.0E-10
+        @timeit "Wilson gauge force" begin
+            force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm)
+        end
+    else
+        @timeit "Improved gauge force" begin
+            force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm, c0)
+        end
+    end
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz bnd_rescale_flw!(ymws.frc1,lp::SpaceParm)
+    end
+
+    return nothing
+end
+
+
 """ 
     function force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm)
 
diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl
index 42ed545..bcb6a35 100644
--- a/src/YM/YMflow.jl
+++ b/src/YM/YMflow.jl
@@ -93,7 +93,7 @@ function Base.show(io::IO, int::FlowIntr{N,T}) where {N,T}
     if N == 0
         println(io, " * Euler schem3")
     elseif N == 1
-        println(io, " * One stage scheme. Coefficients3")
+        println(io, " * One stage scheme. Coefficients")
         println(io, "    stg 1: ", int.e0[1], " ", int.e1[1])
     elseif N == 2
         println(io, " * Two stage scheme. Coefficients:")
@@ -201,6 +201,31 @@ function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpacePar
 end
 flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} = flw(U, int, ns, int.eps, gp, lp, ymws)
 
+function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D}
+    @timeit "Integrating flow equations" begin
+        for i in 1:ns
+            force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
+            if int.add_zth
+                add_zth_term(ymws::YMworkspace, U, lp)
+            end
+            ymws.mom .= ymws.frc1
+            U .= expm.(U, ymws.mom, 2*eps*int.r)
+
+            for k in 1:NI
+                force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
+                if int.add_zth
+                    add_zth_term(ymws::YMworkspace, U, lp)
+                end
+                ymws.mom .= int.e0[k].*ymws.mom .+ int.e1[k].*ymws.frc1
+                U .= expm.(U, ymws.mom, 2*eps)
+            end
+        end
+    end
+
+    return nothing
+end
+flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D} = flw(U, int, ns, int.eps, gp, lp, ymws)
+
 
 ##
 # Adaptive step size integrators
@@ -316,30 +341,30 @@ Eoft_plaq(U, gp::GaugeParm{T,G,NN}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) w
 
 
 function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
-    
+
     @inbounds begin
         b = Int64(CUDA.threadIdx().x)
         r = Int64(CUDA.blockIdx().x)
         I = point_coord((b,r), lp)
-        
+
         id1, id2 = lp.plidx[ipl]
         SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI)) && (id1 == N)
         TWP  = ((I[id1]==1)&&(I[id2]==1))
-        
+
         bu1, ru1 = up((b, r), id1, lp)
         bu2, ru2 = up((b, r), id2, lp)
-        
-        if SFBC && (ru1 != r)
+
+        if SFBC && (point_time((b,r),lp) == lp.iL[end])
             gt = Ubnd[id2]
         else
             gt = U[bu1,id2,ru1]
         end
-        
+
         if TWP
             plx[I] = ztw*tr(U[b,id1,r]*gt / (U[b,id2,r]*U[bu2,id1,ru2]))
         else
             plx[I] = tr(U[b,id1,r]*gt / (U[b,id2,r]*U[bu2,id1,ru2]))
-        end            
+        end
     end
     return nothing
 end
diff --git a/src/YM/YMsf.jl b/src/YM/YMsf.jl
index abd93d9..cd00ac7 100644
--- a/src/YM/YMsf.jl
+++ b/src/YM/YMsf.jl
@@ -92,7 +92,7 @@ end
 """
         function setbndfield(U, phi, lp::SpaceParm)
 
-Sets abelian boundary fields with phases `phi[1]` and `phi[2]` to the configuration `U` at time salice ``x_0=0``.
+Sets abelian boundary fields with phases `phi[1]` and `phi[2]` to the configuration `U` at time slice ``x_0=0``.
 """
 function setbndfield(U, phi, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}
 
diff --git a/test/dirac/test_adapt_ferm.jl b/test/dirac/test_adapt_ferm.jl
new file mode 100644
index 0000000..0368057
--- /dev/null
+++ b/test/dirac/test_adapt_ferm.jl
@@ -0,0 +1,48 @@
+using LatticeGPU, Test, CUDA
+
+T = Float64
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_PERIODIC, (0,0,0,0,0,0))
+gp = GaugeParm{T}(SU3{T}, 6.1, 1.0)
+dpar = DiracParam{T}(SU3fund,1.3,0.9,(1.0,1.0,1.0,1.0),0.0,0.0)
+ymws = YMworkspace(SU3, T, lp)
+dws = DiracWorkspace(SU3fund,T,lp);
+
+randomize!(ymws.mom, lp, ymws)
+U = exp.(ymws.mom)
+
+psi = scalar_field(Spinor{4,SU3fund{T}},lp);
+pfrandomize!(psi,lp)
+
+Ucp = deepcopy(U)
+psicp = deepcopy(psi)
+# First Integrate very precisely up to t=2 (Wilson)
+println(" # Very precise integration ")
+wflw = wfl_rk3(Float64, 0.0004, 1.0E-7)
+flw(U,psi, wflw, 5000, gp,dpar, lp, ymws, dws)
+pl_exact = Eoft_plaq(U, gp, lp, ymws)
+cl_exact = Eoft_clover(U, gp, lp, ymws)
+println("  - Plaq:   ", pl_exact)
+println("  - Clover: ", cl_exact)
+Ufin = deepcopy(U)
+psifin = deepcopy(psi)
+
+
+# Now use Adaptive step size integrator:
+for tol in (1.0E-4, 1.0E-5, 1.0E-6, 1.0E-7, 1.0E-8)
+    local wflw = wfl_rk3(Float64, 0.0001, tol)
+    U .= Ucp
+    psi .= psicp
+    ns, eps = flw_adapt(U,psi, wflw, 2.0, gp,dpar,lp, ymws,dws)
+    pl = Eoft_plaq(U, gp, lp, ymws)
+    cl = Eoft_clover(U, gp, lp, ymws)
+    psierr = sum(norm2.((psi.-psifin)))./prod(lp.iL)
+
+    println(" # Adaptive integrator (tol=$tol): ", ns, " steps")
+    U .= U ./ Ufin
+    maxd = CUDA.mapreduce(dev_one, max, U, init=0.0)
+    println("  - Plaq:   ", pl," [diff: ", abs(pl-pl_exact), "; ",
+            maxd, "]")
+    println("  - Clover: ", cl, " [diff: ", abs(cl-cl_exact), "; ",
+            maxd, "]")
+    println("  - Fermion diff: ", psierr)
+end
diff --git a/test/dirac/test_backflow.jl b/test/dirac/test_backflow.jl
index 601c276..0700b0e 100644
--- a/test/dirac/test_backflow.jl
+++ b/test/dirac/test_backflow.jl
@@ -1,13 +1,8 @@
-using CUDA
+using CUDA, LatticeGPU
 
-using Pkg
-
-Pkg.activate("/home/fperez/Git/LGPU_fork_ferflow")
-
-using LatticeGPU
-
-lp = SpaceParm{4}((4,4,4,4),(2,2,2,2),0,(0,0,0,0,0,0));
+println(" # Consistency condition for backflow")
 
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_PERIODIC, (0,0,0,0,0,0))
 pso = scalar_field(Spinor{4,SU3fund{Float64}},lp);
 psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
 psi2 = scalar_field(Spinor{4,SU3fund{Float64}},lp);
@@ -19,24 +14,27 @@ int = wfl_rk3(Float64, 0.01, 1.0)
 
 gp = GaugeParm{Float64}(SU3{Float64},6.0,1.0,(1.0,0.0),(0.0,0.0),lp.iL)
 
-dpar = DiracParam{Float64}(SU3fund,1.3,0.9,(1.0,1.0,1.0,1.0),0.0)
+dpar = DiracParam{Float64}(SU3fund,1.3,0.9,(1.0,1.0,1.0,1.0),0.0,0.0)
 
 randomize!(ymws.mom, lp, ymws)
 U = exp.(ymws.mom);
 
 pfrandomize!(psi,lp)
-for L in 4:19
+for L in 10:20:210
     pso .= psi
     V = Array(U)
-    a,b = flw_adapt(U, psi, int, L*int.eps, gp,dpar, lp, ymws,dws)
+    #a,b = flw_adapt(U, psi, int, L*int.eps, gp,dpar, lp, ymws,dws)
+    flw(U, psi, int, L,int.eps, gp,dpar, lp, ymws,dws)
      # for i in 1:a
      # flw(U, psi, int, 1 ,b[i], gp, dpar, lp, ymws, dws)
      # end
     pfrandomize!(psi2,lp)
 
-    foo = sum(dot.(psi,psi2))# field_dot(psi,psi2,sumf,lp)
+    foo = sum(dot.(psi,psi2))
     copyto!(U,V);
-    backflow(psi2,U,L*int.eps,7,gp,dpar,lp, ymws,dws)
-    println("Error:",(sum(dot.(pso,psi2))-foo)/foo)
+    backflow(psi2,U,L*int.eps,20,gp,dpar,lp, ymws,dws)
+    println("# Consistency backflow test for t=",L*int.eps)
+    println("Relative error:",abs((sum(dot.(pso,psi2))-foo)/foo))
     psi .= pso
 end
+
diff --git a/test/dirac/test_backflow_tl.jl b/test/dirac/test_backflow_tl.jl
index 2bf0d18..9330147 100644
--- a/test/dirac/test_backflow_tl.jl
+++ b/test/dirac/test_backflow_tl.jl
@@ -3,13 +3,14 @@ using LatticeGPU, CUDA, TimerOutputs
 #Test for the relation K(t,y;0,n)^+ Dw(n|m)^{-1} e^(ipm) = D(p)^{-1} exp(4t sin^2(p/2)) e^{ipn} with a given momenta (if p=0 its randomized), spin and color
 #Kernel en 1207.2096
 
+println(" # Free fermion propagator for backflow")
 
 @timeit "Plw backflow test" begin
 
     function Dwpw_test(;p=0,s=1,c=1)
         lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
         gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
-        dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0)
+        dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0)
         dws = DiracWorkspace(SU3fund,Float64,lp);
         ymws = YMworkspace(SU3,Float64,lp);
 
@@ -89,9 +90,7 @@ using LatticeGPU, CUDA, TimerOutputs
             g5Dw!(prop,U,pwave,dpar,dws,lp)
             CG!(prop,U,DwdagDw!,dpar,lp,dws,10000,1.0e-14)
 
-            for _ in 1:Nsteps
-                backflow(U,prop,1,int.eps,gp,dpar,lp, ymws,dws)
-            end
+            backflow(prop,U,Nsteps*int.eps,20,gp,dpar,lp, ymws,dws)
         end
 
 
@@ -103,15 +102,15 @@ using LatticeGPU, CUDA, TimerOutputs
 
 
     begin
-        dif = 0.0
+        global diff = 0.0
         for i in 1:3 for j in 1:4
-            dif += Dwpw_test(c=i,s=j)
+            global diff += Dwpw_test(c=i,s=j)
         end end
 
-        if dif < 1.0e-5
-            print("Backflow_tl test passed with average error ", dif/12,"!\n")
+        if diff < 1.0e-5
+            print("Backflow_tl test passed with average error ", diff/12,"\n")
         else
-            error("Backflow_tl test failed with difference: ",dif,"\n")
+            error("Backflow_tl test failed with difference: ",diff,"\n")
         end
 
 
diff --git a/test/dirac/test_flow_tl.jl b/test/dirac/test_flow_tl.jl
index 65ed678..b73be6b 100644
--- a/test/dirac/test_flow_tl.jl
+++ b/test/dirac/test_flow_tl.jl
@@ -1,15 +1,16 @@
 using LatticeGPU, CUDA, TimerOutputs
 
 #Test for the relation K(t,y;0,n) Dw(n|m)^{-1} e^(ipm) = D(p)^{-1} exp(-4t sin^2(p/2)) e^{ipn} with a given momenta (if p=0 its randomized), spin and color
-#Kernel en 1207.2096
+#Kernel from 1207.2096
 
+println(" # Free fermion propagator for frontflow")
 
 @timeit "Plw flow test" begin
 
     function Dwpw_test(;p=0,s=1,c=1)
         lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
         gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
-        dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0)
+        dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0)
         dws = DiracWorkspace(SU3fund,Float64,lp);
         ymws = YMworkspace(SU3,Float64,lp);
 
@@ -103,15 +104,15 @@ using LatticeGPU, CUDA, TimerOutputs
 
 
     begin
-        dif = 0.0
+        global diff = 0.0
         for i in 1:3 for j in 1:4
-            dif += Dwpw_test(c=i,s=j)
+            global diff += Dwpw_test(c=i,s=j)
         end end
 
-        if dif < 1.0e-4
-            print("Flow_tl test passed with average error ", dif/12,"!\n")
+        if diff < 1.0e-4
+            print("Flow_tl test passed with average error ", diff/12,"\n")
         else
-            error("Flow_tl test failed with difference: ",dif,"\n")
+            error("Flow_tl test failed with difference: ",diff,"\n")
         end
 
 
diff --git a/test/dirac/test_fp_fa.jl b/test/dirac/test_fp_fa.jl
index 5551da5..a7fa00c 100644
--- a/test/dirac/test_fp_fa.jl
+++ b/test/dirac/test_fp_fa.jl
@@ -2,6 +2,8 @@ using LatticeGPU
 using CUDA
 using TimerOutputs
 
+println(" # Free solution for SF correlation functions")
+
 @timeit "fA_fP test" begin
 
 
@@ -115,7 +117,7 @@ using TimerOutputs
         elseif difP > 1.0e-15
             error("fP test failed with error ", difP)
         else
-            print("fA & fP tests passed with errors: ", difA," and ",difP,"!\n")
+            print("fA & fP tests passed with errors: ", difA," and ",difP,"\n")
         end
 
 end
diff --git a/test/dirac/test_solver_plw.jl b/test/dirac/test_solver_plw.jl
index a4de0c7..7abc62a 100644
--- a/test/dirac/test_solver_plw.jl
+++ b/test/dirac/test_solver_plw.jl
@@ -2,6 +2,8 @@ using LatticeGPU, CUDA, TimerOutputs
 
 #Test for the relation Dw(n|m)^{-1} e^(ipm) = D(p)^{-1} e^{ipn} with a given momenta (if p=0 its randomized), spin and color
 
+println(" # Test for free fermion propagator")
+
 @timeit "Plw test" begin
 
 function Dwpw_test(;p=0,s=1,c=1)
@@ -84,12 +86,12 @@ end
 
 dif = sum(norm2.(prop - prop_th))
 
-if dif > 1.0e-15
+if dif > 1.0e-7
         error("Dwpl test for s=",s,", c=",c," failed with difference: ",dif,"\n")
 end
 
 
-return dif
+return sqrt(dif)
 end
 
 
@@ -101,8 +103,8 @@ for i in 1:3 for j in 1:4
     global diff += Dwpw_test(c=i,s=j)
 end end
 
-if diff < 1.0e-15
-    print("Dwpl test passed with average error ", diff/12,"!\n")
+if diff < 1.0e-7
+    print("Dwpl test passed with average error ", diff/12,"\n")
 else
     error("Dwpl test failed with difference: ",diff,"\n")
 end
diff --git a/test/dirac/test_solver_rand.jl b/test/dirac/test_solver_rand.jl
index 0714d8f..3370e67 100644
--- a/test/dirac/test_solver_rand.jl
+++ b/test/dirac/test_solver_rand.jl
@@ -2,6 +2,9 @@ using CUDA, LatticeGPU, TimerOutputs
 
 #Check that Dw ( (DwdagDw)^{-1} g5 Dw g5 ) psi = psi for random fields
 
+println(" # Test for the consistency of the solver")
+
+
 @timeit "Rand solver test" begin
 
 @timeit "Generate random fields" begin
@@ -46,7 +49,7 @@ res = sum(norm2.(rpsi-dws.sp))
 
 
 if res < 1.0e-6 
-    print("Drand test passed with ",res,"% error!\n")
+    print("Drand test passed with ",res,"% error\n")
     
 else
     error("Drand test failed with difference: ",res,"\n")
diff --git a/test/runtests.jl b/test/runtests.jl
index 2f68f70..6ffdddf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,10 @@
 
-include("SAD/test_sad.jl")
+#include("SAD/test_sad.jl")
 include("flow/test_adapt.jl")
 include("dirac/test_fp_fa.jl")
 include("dirac/test_solver_plw.jl")
 include("dirac/test_solver_rand.jl")
+include("dirac/test_flow_tl.jl")
+include("dirac/test_backflow_tl.jl")
+include("dirac/test_backflow.jl")
+include("dirac/test_adapt_ferm.jl")