/** This class tests decomposition routines in Matrix.frink */

use Matrix.frink

A = new Matrix[[[  4,  12, -16],
                [ 12,  37, -43],
                [-16, -43,  98]]]

A = new Matrix[[[  25,  15, -5],
                [ 15,  18, 0],
                [-5, 0,  11]]]

println["Original matrix is:\n" + A.formatMatrix[]]


cb = A.CholeskyB[]
println[]
println["Cholesky-Banachiewicz decomposition is:"]
println[cb.formatMatrix[]]

Aprime = cb.multiply[cb.conjugateTranspose[]]

print["Cholesky-Banachiewicz reversal is: "]
if A.equals[Aprime]
   println["OK"]
else
   println["FAIL:\n"+Aprime.formatMatrix[]]


cc = A.CholeskyCrout[]
println[]
println["Cholesky-Crout decomposition is:"]
println[cc.formatMatrix[]]

Aprime = cc.multiply[cc.conjugateTranspose[]]

print["Cholesky-Crout reversal is: "]
if A.equals[Aprime]
   println["OK"]
else
   println["FAIL:\n"+Aprime.formatMatrix[]]


[L,U] = A.LUDecomposeCrout[]
println[]
println["LUDecomposeCrout decomposition is:"]
println[formatTable[[[L.formatMatrix[],U.formatMatrix[]]]]]

Aprime = L.multiply[U]

print["LUDecomposeCrout reversal is: "]
if A.equals[Aprime]
   println["OK"]
else
   println["FAIL:\n"+Aprime.formatMatrix[]]



[L,U] = A.LUDecomposeCrout2[]
println[]
println["LUDecomposeCrout2 decomposition is:"]
println[formatTable[[[L.formatMatrix[],U.formatMatrix[]]]]]

Aprime = L.multiply[U]

print["LUDecomposeCrout2 reversal is: "]
if A.equals[Aprime]
   println["OK"]
else
   println["FAIL:\n"+Aprime.formatMatrix[]]