#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# easycnm - Easy Computer Numerics Library
# Copyright (C) 2005 Josef Spillner <js177634@inf.tu-dresden.de>
# Published under GNU GPL conditions

import Numeric
import cnmlib

a = Numeric.array(((5, 2), (2, 4)), Numeric.Float)
b = Numeric.array((4, 16))
b = Numeric.reshape(b, (2, 1))

print "Gegeben sei die Matrix A (Gleichung Ax=b):"
print a

""" ============================================= """

print "Gaußsches Eliminationsverfahren und Invertierungsverfahren"

g = Numeric.concatenate((a, b), 1)
gstaffel = cnmlib.matrix_gauss(g)

print "Matrix (A|b) in der gestaffelten Form:"
print gstaffel

x = cnmlib.matrix_solve(gstaffel)

print "Der Lösungsvektor x ist:"
print x

i = Numeric.identity(2)
g = Numeric.concatenate((a, i), 1)
gstaffel = cnmlib.matrix_gauss(g)

print "Matrix A (gestaffelt) nach Gauss-Vorstufe zur Inversen:"
print gstaffel

ainv = cnmlib.matrix_invert(gstaffel)

print "Inverse Matrix zu A:"
print ainv

""" ============================================= """

print "LU-Zerlegung"
(l, u) = cnmlib.matrix_lu(a)

print "L und U:"
print l
print u

x = cnmlib.matrix_solve_lu(l, u, b)

print "Lösungsvektor x:"
print x

""" ============================================= """

print "Cholesky-Zerlegung"
(l, d, lt) = cnmlib.matrix_cholesky(a)

print "L und D und LT:"
print l
print d
print lt

x = cnmlib.matrix_solve_cholesky(l, d, lt, b)

print "Lösungsvektor x"
print x

""" ============================================= """

print "Jacobi-Verfahren (Gesamtschrittverfahren)"

x = Numeric.array((0, 3))
x = Numeric.reshape(x, (2, 1))

print "Startvektor x:"
print x

x = cnmlib.matrix_jacobi_start(a, b, x)

print "Lösung x nach 1. Schritt:"
print x

x = cnmlib.matrix_jacobi_step()

print "Lösung x nach 2. Schritt:"
print x

x = cnmlib.matrix_jacobi_step()

print "Lösung x nach 3. Schritt:"
print x

x = cnmlib.matrix_jacobi_step()

print "Lösung x nach 4. Schritt:"
print x

""" ============================================= """

print "Gauss-Seidel-Verfahren (Einzelschrittverfahren)"

x = Numeric.array((0, 3))
x = Numeric.reshape(x, (2, 1))

print "Startvektor x:"
print x

x = cnmlib.matrix_gaussseidel_start(a, b, x)

print "Lösung x nach 1. Schritt:"
print x

x = cnmlib.matrix_gaussseidel_step()

print "Lösung x nach 2. Schritt:"
print x

x = cnmlib.matrix_gaussseidel_step()

print "Lösung x nach 3. Schritt:"
print x

""" ============================================= """

print "Methode der kleinsten Quadrate (Normalgleichungsverfahren)"

x = Numeric.array((0, 3, 4))
x = Numeric.reshape(x, (3, 1))
y = Numeric.array((1, 8, 10))
y = Numeric.reshape(y, (3, 1))

ab = cnmlib.matrix_ngv(x, y)

print "Lösungsvektor (a, b)^T:"
print ab

""" ============================================= """

print "QR-Zerlegung (Householder-Algorithmus)"

a = Numeric.array(((12, -51, 4), (6, 167, -68), (-4, 24, -41)))

(q, r) = cnmlib.matrix_householder(a)

print "Lösungsmatrizen Q, R:"
print q
print r