Shape preservation for blocked cholesky

ring.jl

import Base: ctranspose, *, -, *, \, zero, showlimited, print, show, writemime, showcompact, showcompact_lim, inv, one
import Base.display
import Base.LinAlg.chol!
import Base.LinAlg.chol


const Dadd = Dict(
      ('L','L') => 'L', 
      ('U','U') => 'U',
      ('L','U') => 'A',
 
      ('A','L') => 'A',
      ('A','U') => 'A',
      ('A','A') => 'A',

      ('0','L') => 'L',
      ('0','U') => 'U',
      ('0','A') => 'A',
      ('0','F') => 'F',
      ('0','0') => '0',

      ('A', 'F') => 'F',
      ('F', 'F') => 'F'


      )

const Dmul = Dict(
      ('L','L') => 'L', 
      ('U','U') => 'U', 

      ('L','U') => 'A', 
      ('L','A') => 'A',
      ('A','U') => 'A',

      ('A', '0') => '0',
      ('F', '0') => '0',
      ('L', '0') => '0',
      ('U', '0') => '0',
      ('0', '0') => '0',
      ('0', 'A') => '0',
      ('0', 'L') => '0',
      ('0', 'F') => '0',
      ('0', 'U') => '0',
      ('F', 'F') => 'F',
      ('A', 'A') => 'F',



      )

const Dldiv = Dict(
      ('L','A') => 'A', 
      ('L','0') => '0',
      ('L','L') => 'L',
      ('L','F') => 'F',
      ('F','F') => 'F',
      ('F','0') => '0',
      ('F','1') => 'F',
      ('L','1') => 'L',
      ('A','1') => 'F',
      ('U','1') => 'U',
      ('A','A') => 'A', 

      )
     
     
type Ring
   x::Float64
   c::Char
end

zero(::Type{Ring}) = Ring(0.0, '0')
one(::Type{Ring}) = Ring(1, 'A')

print(io::IO, r::Ring) = print(io,r.c)

showlimited(io::IO, r::Ring) = print(io, r)
writemime(io::IO, ::MIME"text/plain", r::Ring) = print(io, r)
showcompact_lim(io::IO, r::Ring) = print(io, r)
showcompact(io::IO, r::Ring) = print(io, r)


ctranspose(r::Ring) = Ring(r.x, Dict('L'=>'U', 'U'=>'L', 'A'=>'A', '0'=>'0', 'F'=>'F')[r.c])

+(r1::Ring, r2::Ring) = Ring(r1.x + r2.x, Dadd[minmax(r1.c, r2.c)])
-(r1::Ring, r2::Ring) = Ring(r1.x - r2.x, Dadd[minmax(r1.c, r2.c)])
*(r1::Ring, r2::Ring) = Ring(r1.x * r2.x, Dmul[r1.c, r2.c])
*(f1::Float64, r2::Ring) = Ring(f1 * r2.x, r2.c)
\(r1::Ring, r2::Ring) = Ring(r1.x \ r2.x, Dldiv[r1.c, r2.c])
inv(r::Ring) = Ring(inv(r.x), Dldiv[r.c, '1'])


function root(r::Ring) 
#    assert(r.c == 'A' || r.c == '0')
    Ring(sqrt(r.x), Dict('A'=>'U', '0'=>'0', 'F'=>'F')[r.c])      
end

chol(r::Ring, v) = chol!(r, v)
chol(r::Ring) = chol!(r, Val{:L})

function chol!(r::Ring, ::Type{Val{:U}})
    Ring(sqrt(r.x), Dict('A'=>'U', '0'=>'0', 'F'=>'F')[r.c])      
end

function chol!(r::Ring, ::Type{Val{:L}})
    Ring(sqrt(r.x), Dict('A'=>'L', '0'=>'0', 'F'=>'F')[r.c])
end

 
L = Ring[Ring(3.,'L') Ring(1.,'L')
         Ring(1.,'L') Ring(3.,'L')]
A = L*(L')
chol!(copy(A))

F = Float64[
1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  3  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  2  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  3  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  2  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0
1  2  1  1  1  1  1  1  1  2  1  1  2  1  9  0  0
2  1  1  1  1  1  1  1  2  1  1  2  1  2  7  9  0
1  1  1  1  1  1  2  1  1  2  1  9  1  1  7  4  9
]

n = size(F,1)
R = randn(n,n)
F = R .* F

U = F * F' + n*eye(n)

A2 = map(x-> Ring(x, x != 0 ? 'A':'0'), U)
L2 = chol!(copy(A2))
display(A2)
display(L2)

UU = rand(n*n, n*n)

1(abs(chol!(UU*UU' .* kron(U, U) + kron(U,U) + 100*n*n*eye(n*n))) .> eps() )

Thanks for the link. I like the example and surprisingly enough, in this week, I've actually looked at the shape preserving property of the arrow matrix under Cholesky/LDLt factorizations. We might need the backward Cholesky at some point such that we can get the same property when the arrowhead is to the upper left.