Typo one

src/Dirac/Dirac.jl

@@ -53,7 +53,7 @@ struct DiracWorkspace{T}
 end
 export DiracWorkspace, DiracParam
 
-function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMWorkspace) where {B,D}
+function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
 
     if B == BC_SF_AFWB || B == BC_SF_ORBI
 
@@ -94,7 +94,7 @@ function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}, ymws::YMWorksp
     return nothing
 end
 
-function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws::YMWorkspace) where {B,D}
+function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
 
     if B == BC_SF_AFWB || B == BC_SF_ORBI
         @timeit "DwdagDw" begin
@@ -183,7 +183,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}, ymws:
 end
 
 
-function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}, ymws::YMWorkspace) where {B,D}
+function g5Dw!(so, U, si, dpar, lp::SpaceParm{4,6,B,D}, ymws::YMworkspace) where {B,D}
 
     if B == BC_SF_AFWB || B == BC_SF_ORBI
         @timeit "SF fix" begin

src/Solvers/CG.jl

@@ -34,7 +34,7 @@ return sum(sumf)
 end
 
 
-function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, ymws::YMWorkspace, maxiter::Int64 = 10, tol=1.0) where {T}
+function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, ymws::YMworkspace, maxiter::Int64 = 10, tol=1.0) where {T}
 
     dws.sr .= si
     dws.sp .= si
@@ -78,3 +78,68 @@ function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T},
     
     return niter
 end
+
+
+function CG_Dw!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, ymws = "ymws", gp = "gp", maxiter::Int64 = 10, tol=1.0) where {T}
+
+    # r -> sr , s -> sp , p -> si , B^+ r -> st
+
+    iterations = 0
+    sumf = scalar_field(Complex{T}, lp) 
+
+    tol = tol * CUDA.mapreduce(x -> norm2(x), +, si)
+
+    nini = CUDA.mapreduce(x -> norm2(x), +, si)
+
+    dws.sr .= si
+    fill!(si,zero(eltype(si)))
+    err = 0.0
+
+    A(dws.sp, U, dws.sr, dpar, lp, ymws, gp)
+    
+
+    dws.st .= dws.sp
+
+    Brn = CUDA.mapreduce(x -> norm2(x), +, dws.st)
+
+    niter = 0
+
+    print("Normalized tol (tol*|v|^2) = ", tol ,"\n")
+
+    for i in 1:maxiter
+        A(dws.sAp, U, dws.sp, dpar, lp, ymws, gp)
+
+        prod = field_dot(dws.sAp,dws.sAp,sumf,lp)
+        
+        alpha = Brn/prod
+
+        si     .= si     .+ alpha .* dws.sp
+        dws.sr .= dws.sr .- alpha .* dws.sAp
+
+        A(dws.st, U, dws.sr, dpar, lp, ymws, gp)
+
+        err = CUDA.mapreduce(x -> norm2(x), +, dws.sr)
+        
+        print(i," : |r_n|/|v| = ", sqrt(err/nini) , "\n")
+
+        if err < tol
+	    niter = i
+            break
+        end
+
+        Brn = CUDA.mapreduce(x -> norm2(x), +, dws.st)
+
+        beta = Brn/(alpha*prod)
+        dws.sp .= dws.st .+ beta .* dws.sp
+        
+    end
+
+    if err > tol
+	error("CG! not converged after $maxiter iterations (Residuals: $err)")
+    end
+    
+    return niter
+end
+
+
+