MatrixSolver/form.jl

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