Testing static condensation for simple matrices

In [4]:

    

function make_matrix(n, s=5.0)
    A = randn(n,n)
    for i = 1:n
        A[i,i] += s
    end
    A
end


    

    
Out[4]:





make_matrix (generic function with 2 methods)

In [5]:

    

nonsymm2symm(A) = (A + A')/2.0


    

    
Out[5]:





nonsymm2symm (generic function with 1 method)

In [6]:

    

make_symmat(n) = nonsymm2symm(make_matrix(n))


    

    
Out[6]:





make_symmat (generic function with 1 method)

In [7]:

    

n = 5
A = make_symmat(n)


    

    
Out[7]:





5x5 Array{Float64,2}:
  3.22619   -0.305419   0.346868  -1.51704   -1.04451 
 -0.305419   6.01351    0.50301   -1.14465   -0.614683
  0.346868   0.50301    5.52941    0.50905   -0.629272
 -1.51704   -1.14465    0.50905    5.96433    0.952459
 -1.04451   -0.614683  -0.629272   0.952459   4.07297 

In [9]:

    

type MatrixDecomp
    nb::Int
    ni::Int
    Abb::Array{Float64,2}
    Aii::Array{Float64,2}
    Abi::Array{Float64,2}
    Aib::Array{Float64,2}
    
    function MatrixDecomp(A, nb, ni)
        Abb = A[1:nb, 1:nb]
        Abi = A[1:nb, (nb+1):end]
        Aib = A[(nb+1):end, 1:nb]
        Aii = A[(nb+1):end, (nb+1):end]
        
        new(nb, ni, Abb, Aii, Abi, Aib)
    end
end


    

In [10]:

    

X = MatrixDecomp(A, 3, 2)


    

    
Out[10]:





MatrixDecomp(3,2,3x3 Array{Float64,2}:
  3.22619   -0.305419  0.346868
 -0.305419   6.01351   0.50301 
  0.346868   0.50301   5.52941 ,2x2 Array{Float64,2}:
 5.96433   0.952459
 0.952459  4.07297 ,3x2 Array{Float64,2}:
 -1.51704  -1.04451 
 -1.14465  -0.614683
  0.50905  -0.629272,2x3 Array{Float64,2}:
 -1.51704  -1.14465    0.50905 
 -1.04451  -0.614683  -0.629272)

In [12]:

    

using Base.LinAlg.LAPACK.getrf!
using Base.LinAlg.LAPACK.getrs!
using Base.LinAlg.BLAS.gemm!
type StaticCond
    nb::Int
    ni::Int
    B::Array{Float64,2}
    iAii::Array{Float64,2}
    X::Array{Float64,2}
    Y::Array{Float64,2}
    function StaticCond(A::MatrixDecomp)
        nb = A.nb
        ni = A.ni
        B = zeros(nb, nb)
        iAii = copy(A.Aii)
        iAii, ipiv, info = getrf!(iAii)
        X = A.Abi'
        getrs!('T', iAii, ipiv, X)
        Y = copy(A.Aib)
        getrs!('N', iAii, ipiv, Y)
        
        B = copy(A.Abb)
        gemm!('T', 'N', -1.0, X, A.Aib, 1.0, B)
        
        new(nb, ni, B, iAii, X, Y)
    end
end


    

In [ ]: