form.jl

@@ -0,0 +1,146 @@
+function form{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,1}; kwargs...)
+	Ax = zeros(promote_type(T1,T2),size(A,1))
+	return form(x -> A_mul_B!(1.0,A,x,0.0,Ax),b)
+end
+
+"""
+x,flag,err,iter,resvec = form(A,b,tol=1e-2,maxIter=100,M=1,x=[],out=0)
+
+Generalized Minimal residual (GMRESm) method with  restarts applied to A*x = b.
+
+Input:
+
+A                - function computing A*xor
+
+M                -preconditioner,function computing M\\x
+"""
+function form(A::Function,b::Vector,tol::Real=1e-2,maxIter::Int=100,M::Function=identity,x::Vector=[],out::Int=2,storeInterm::Bool=false)
+    n = length(b)
+    if norm(b)==0;return zeros(eltype(b),n),-9,0.0,0,[0.0];end
+    if isempty(x)
+        x = zeros(n)
+        r = M(b)
+    else
+        r = M(b-A(x))
+    end
+    if storeInterm
+        x = zeros(n,maxIter)
+    end
+
+    bnrm2 = norm(b)
+    if bnrm2 == 0.0;bnrm2 = 1.0;end
+
+    err = norm( r ) / bnrm2
+
+    if err < tol; return x,err;end
+
+    #Arnoldi向量 m+1
+    V = zeros(n,maxIter+1)
+    #m+1,m
+    H = zeros(maxIter+1,maxIter)
+
+    #e1
+    e1 = zeros(n)
+    e1[1] = 1.0
+
+    #暂时不考虑complex
+
+    #残差数组
+    resvec = zeros(maxIter)
+    if out==2
+        println(@sprintf("=== gmres ===\n%4s\t%7s\n","iter","relres"))
+    end
+
+    #初始化
+    iter = 0
+    flag = -1
+    cnt = 1
+    y = 0
+    # v1 = r0/belta
+    V[:,1] = r / norm(r)
+    #theta = belta*e1
+    s = norm(r)*e1;
+    #显示迭代次数
+    i = 1
+    for i = 1:maxIter
+        if out==2;;print(@sprintf("%3d\t",i));end
+        #wi = Avi
+        w = A(V[:,i])
+        w = M(w)
+
+        #Arnoldi process
+        for k = 1:i
+            # hij = （wj,vi)
+            H[k,i] = dot(w,V[:,k])
+            # wj = wj - hijVi
+            w -= H[k,i] * V[:,k]
+        end
+        #hj+1,j = ||wj||2
+        H[i+1,i] = norm(w)
+
+        #求解方程Hjy = belta* e1
+        y  = H[1:i,1:i]\s[1:i]
+        err = H[i+1,i]*abs(y[i])/bnrm2
+        if out == 2;print(@sprintf("%1.1e\n",err));end
+        resvec[cnt] = err
+        cnt = cnt + 1
+        if  err <= tol
+            flag = 0
+            break
+        end
+        #vj+1 = wj/h(j+1,j)
+        V[:,i+1] = w/H[i+1,i]
+    end
+    x += V[:,1:i]*y
+
+    r = b - A(x)
+    r = M(r)
+
+
+    if out==2; print(@sprintf("\t %1.1e\n", err)); end
+
+    if out>=0
+        if flag==-1
+            println(@sprintf("gmres iterated maxIter (=%d) times without achieving the desired tolerance.",maxIter))
+        elseif flag==0 && out>=1
+            println(@sprintf("gmres achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,resvec[cnt]))
+        end
+    end
+    return x,flag,resvec[cnt-1],iter,resvec[1:cnt-1]
+end
+
+
+"""
+c,s,r = SymOrtho(a,b)
+
+Computes a Givens rotation
+
+Implementation is based on Table 2.9 in
+
+Choi, S.-C. T. (2006).
+Iterative Methods for Singular Linear Equations and Least-squares Problems.
+Phd thesis, Stanford University.
+"""
+function symOrtho(a,b)
+    c = 0.0; s = 0.0; r = 0.0
+	if b==0
+		s = 0.0
+		r = abs(a)
+		c = (a==0) ? c=1.0 : c = sign(a)
+	elseif a == 0
+		c = 0.0
+		s = sign(b)
+		r = abs(b)
+	elseif abs(b) > abs(a)
+		tau = a/b
+		s   = sign(b)/sqrt(1+tau^2)
+		c   = s*tau
+		r   = b/s
+	elseif abs(a) > abs(b)
+		tau = b/a
+		c   = sign(a)/sqrt(1+tau^2)
+		s   = c*tau
+		r   = a/c
+	end
+	return c,s,r
+end

gmres.jl

@@ -0,0 +1,194 @@
+function gmres{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,1},restrt::Int; kwargs...)
+	Ax = zeros(promote_type(T1,T2),size(A,1))
+	return gmres(x -> A_mul_B!(1.0,A,x,0.0,Ax),b,restrt;kwargs...)
+end
+
+"""
+x,flag,err,iter,resvec = gmres(A,b,restrt,tol=1e-2,maxIter=100,M=1,x=[],out=0)
+
+Generalized Minimal residual (GMRESm) method with  restarts applied to A*x = b.
+
+Input:
+
+A                - function computing A*xor
+
+M                -preconditioner,function computing M\\x
+"""
+function gmres(A::Function,b::Vector,restrt::Int,tol::Real=1e-2,maxIter::Int=100,M::Function=identity,x::Vector=[],out::Int=0,storeInterm::Bool=false)
+    n = length(b)
+    if norm(b)==0;return zeros(eltype(b),n),-9,0.0,0,[0.0];end
+    if isempty(x)
+        x = zeros(n)
+        r = M(b)
+    else
+        r = M(b-A(x))
+    end
+    if storeInterm
+        x = zeros(n,maxIter)
+    end
+
+    bnrm2 = norm(b)
+    if bnrm2 == 0.0;bnrm2 = 1.0;end
+
+    err = norm( r ) / bnrm2
+
+    if err < tol; return x,err;end
+
+    #重启次数 m
+    restrt = min(restrt,n-1)
+    #Arnoldi向量 m+1
+    V = zeros(n,restrt+1)
+    #m+1,m
+    H = zeros(restrt+1,restrt)
+    #Gj的cs
+    cs = zeros(restrt)
+    #Gj的sn
+    sn = zeros(restrt)
+    #e1
+    e1 = zeros(n)
+    e1[1] = 1.0
+
+    #暂时不考虑complex
+
+    #残差数组
+    resvec = zeros((1+restrt)*maxIter)
+    if out==2
+        println(@sprintf("=== gmres ===\n%4s\t%7s\n","iter","relres"))
+    end
+
+    #初始化
+    iter = 0
+    flag = -1
+    cnt = 1
+
+    for iter = 1:maxIter
+        # v1 = r0/belta
+        V[:,1] = r / norm(r)
+        #theta = belta*e1
+        s = norm(r)*e1;
+        #显示迭代次数
+        if out==2;;print(@sprintf("%3d\t",iter));end
+        #开始迭代 带重启
+        for i = 1:restrt
+            #wi = Avi
+            w = A(V[:,i])
+            w = M(w)
+
+            #Arnoldi process
+            for k = 1:i
+                # hij = （wj,vi)
+                H[k,i] = dot(w,V[:,k])
+                # wj = wj - hijVi
+                w -= H[k,i] * V[:,k]
+            end
+            #hj+1,j = ||wj||2
+            H[i+1,i] = norm(w)
+
+
+            #Apply Givens rotation
+            for k = 1:i-1
+                temp = cs[k]*H[k,i] + sn[k]*H[k+1,i]
+                H[k+1,i] = -sn[k] *H[k,i] + cs[k]*H[k+1,i]
+                H[k,i] = temp
+            end
+
+            #跳出
+            if H[i+1,i] == 0
+
+            end
+
+            #vj+1 = wj/h(j+1,j)
+            V[:,i+1] = w/H[i+1,i]
+
+            #From the Givens rotation Gj
+            cs[i],sn[i] = symOrtho(H[i,i],H[i+1,i])
+
+            #Apply Gj to last column of Hj+1,j
+            H[i,i] = cs[i]*H[i,i]+sn[i,i]*H[i+1,i]
+            H[i+1,i] = 0.0
+
+            #Apply Gj to right-hand side
+            s[i+1] = -sn[i]*s[i]
+            s[i] = cs[i]*s[i]
+
+            #check convergence
+            err = abs(s[i+1])/bnrm2
+
+            if out == 2;print(@sprintf("%1.1e",err));end
+
+            resvec[cnt] = err
+
+            if err <= tol
+                y = H[1:i,1:i] \ s[1:i]
+                x += V[:,1:i]*y
+                if out == 2;print("\n"); end
+                flag = 0;break;
+            end
+            cnt = cnt + 1
+        end
+        if  err <= tol
+            flag = 0
+            break
+        end
+        y  = H[1:restrt,1:restrt]\s[1:restrt]
+        x += V[:,1:restrt]*y
+        if storeInterm; X[:,iter] = x; end
+
+        r = b - A(x)
+        r = M(r)
+
+        s[restrt+1] = norm(r)
+        resvec[cnt] = abs(s[restrt+1]) / bnrm2
+
+        if out==2; print(@sprintf("\t %1.1e\n", err)); end
+
+    end
+    if out>=0
+        if flag==-1
+            println(@sprintf("gmres iterated maxIter (=%d) times without achieving the desired tolerance.",maxIter))
+        elseif flag==0 && out>=1
+            println(@sprintf("gmres achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,resvec[cnt]))
+        end
+    end
+    if storeInterm
+	    return X[:,1:iter],flag,resvec[cnt],iter,resvec[1:cnt]
+	else
+	    return x,flag,resvec[cnt],iter,resvec[1:cnt]
+	end
+end
+
+
+"""
+c,s,r = SymOrtho(a,b)
+
+Computes a Givens rotation
+
+Implementation is based on Table 2.9 in
+
+Choi, S.-C. T. (2006).
+Iterative Methods for Singular Linear Equations and Least-squares Problems.
+Phd thesis, Stanford University.
+"""
+function symOrtho(a,b)
+    c = 0.0; s = 0.0; r = 0.0
+	if b==0
+		s = 0.0
+		r = abs(a)
+		c = (a==0) ? c=1.0 : c = sign(a)
+	elseif a == 0
+		c = 0.0
+		s = sign(b)
+		r = abs(b)
+	elseif abs(b) > abs(a)
+		tau = a/b
+		s   = sign(b)/sqrt(1+tau^2)
+		c   = s*tau
+		r   = b/s
+	elseif abs(a) > abs(b)
+		tau = b/a
+		c   = sign(a)/sqrt(1+tau^2)
+		s   = c*tau
+		r   = a/c
+	end
+	return c,s,r
+end

jacobi.jl

@@ -0,0 +1,35 @@
+function jacobi(A::SparseMatrixCSC,b::Vector;tol::Real=1e-2,maxIter::Int=1,out::Int=0,storeInterm::Bool=false;)
+	n = length(b)
+	OmInvD = 1./diag(A) #D矩阵 并乘以omega
+	x = zeros(eltype(A),n)
+	r = copy(b)
+	rnorm0 = norm(b)
+	res = 0
+	value = 0.0
+	flag   = -1
+	iter   = 0
+	dx = 0
+	for iter = 1:maxIter
+		println(iter)
+        dr = zeros(eltype(A),n)
+		for j = 1:n
+			dx =  r[j] * OmInvD[j]
+            x[j] += dx
+            for gidx = A.colptr[j] : A.colptr[j+1]-1
+                dr[A.rowval[gidx]] -= A.nzval[gidx] * dx
+            end
+		end
+        r += dr
+
+		res = norm(r) / rnorm0
+
+		println(res)
+
+		if res < tol
+			flag = 0; break
+		end
+	end
+
+	return x,flag,iter
+
+end

sor.jl

@@ -0,0 +1,59 @@
+
+"""
+x,flag,err,iter,resvec = sor(A::SparseMatrixCSC, b::Vector; omega::Real=2/3, tol::Real=1e-2, maxIter::Int=1,out::Int=0)
+
+Symmetric successive over-relaxation applied to the linear system A*x = b.
+
+!! Note that upper/lower triangular matrix are never build to save memory !!
+
+Input:
+
+	A       - sparse matrix CSC
+	b       - right hand side vector
+	omega   - relaxation parameter (omega=1.0: Gauss Seidel)
+	tol     - error tolerance (default=1e-2)
+	maxIter - maximum number of iterations (default=1)
+	out     - flag for output (-1 : no output, 0 : only errors, 1 : final status, 2: error at each iteration)
+
+Output:
+
+	x       - solution
+	flag    - exit flag (  0 : desired tolerance achieved,
+	                      -1 : maxIter reached without converging)
+	err     - error, i.e., norm(A*x-b)/norm(b)
+	iter    - number of iterations
+	resvec  - error at each iteration
+"""
+function sor(A::SparseMatrixCSC,b::Vector;tol::Real=1e-2,maxIter::Int=1,omega::Real=1.0,out::Int=0,storeInterm::Bool=false;)
+	n = length(b)
+	OmInvD = omega./diag(A) #D矩阵 并乘以omega
+	x = zeros(eltype(A),n)
+	r = copy(b)
+	rnorm0 = norm(b)
+	res = 0
+	value = 0.0
+	flag   = -1
+	iter   = 0
+	dx = 0
+	for iter = 1:maxIter
+		println(iter)
+		#第一次迭代 按照列来
+		for j = 1:n
+			dx =  r[j] * OmInvD[j]
+			for gidx = A.colptr[j] : A.colptr[j+1]-1
+				r[A.rowval[gidx]] -= A.nzval[gidx] * dx
+			end
+			x[j] += dx
+		end
+		res = norm(A*x-b) / rnorm0
+
+		println(res)
+
+		if res < tol
+			flag = 0; break
+		end
+	end
+
+	return x,flag,iter
+
+end

ssor.jl

@@ -29,48 +29,39 @@ function ssor(A::SparseMatrixCSC,b::Vector;tol::Real=1e-2,maxIter::Int=1,omega::
 	n = length(b)
 	OmInvD = omega./diag(A) #D矩阵 并乘以omega
 	x = zeros(eltype(A),n)
-	r = b - A * x
-	rnorm0 = norm(r)
+	r = copy(b)
+	rnorm0 = norm(b)
+	res = 0
 	value = 0.0
 	flag   = -1
 	iter   = 0
+	dx = 0
 	for iter = 1:maxIter
 		println(iter)
-		#第一次迭代
-		for i = 1:n
-			value = b[i]
-			for j = 1:n
-				for gidx = A.colptr[j] : A.colptr[j+1]-1
-					if A.rowval[gidx] == i
-						value -= A.nzval[gidx] * x[j]
-					end
-				end
+		#第一次迭代 按照列来
+		for j = 1:n
+			dx =  r[j] * OmInvD[j]
+			for gidx = A.colptr[j] : A.colptr[j+1]-1
+				r[A.rowval[gidx]] -= A.nzval[gidx] * dx
 			end
-			x[i]  = x[i] + value*OmInvD[i]
+			x[j] += dx
 		end
-
-		#第二次迭代
-		for i = n:-1:1
-			value = b[i]
-			for j = 1:n
-				for gidx = A.colptr[j] : A.colptr[j+1]-1
-					if A.rowval[gidx] == i
-						value -= A.nzval[gidx] * x[j]
-					end
-				end
+		for j = n:-1:1
+			dx =  r[j] * OmInvD[j]
+			for gidx = A.colptr[j] : A.colptr[j+1]-1
+				r[A.rowval[gidx]] -= A.nzval[gidx] * dx
 			end
-			x[i]  = x[i] + value*OmInvD[i]
+			x[j] += dx
 		end
-		r = b - A*x
-		println(norm(r) / rnorm0 )
+		res = norm(A*x-b) / rnorm0
+
+		println(res)
 
-		if norm(r) / rnorm0 < tol
+		if res < tol
 			flag = 0; break
 		end
 	end
 
-	println(r)
-	println(x)
-	return x,flag, norm(r) / rnorm0,iter
+	return x,flag,iter
 
 end