#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 29 11:18:44 2017

@author: troll
"""

"""Dies ist eine Lösung für Aufgabe 23.
Sie dürfen rueckwaertseinsetzen(R, b): verwenden um Aufgabe 27 zu lösen.
Sie dürfen KEINE weiteren Features in rueckwaertseinsetzen(R, b) einbauen durch
die Aufgabe 27 leichter werden würde.
"""

import numpy as np

#Infos gibt es mit help(rueckwaertseinsetzen
def rueckwaertseinsetzen(R, b):
    """
    Löst ein Gleichungssystem Rx=b mit oberer Dreiecksmatrix R und Vektor b
    Eingaben:
    R: Obere Dreicksmatriz als array (wird nicht überprüft)
    b: Rechte Seite des LGS (als array)
    Rückgabe:
    x: Lösung des Gleichungssystems Rx=b
    
    Sollte R singulär sein (einer der Diagonaleinträge ist 0), so wird eine
    Warnung ausgegeben. Die Rückgabe ist das ein 0-array (Größe wie b)
    """
    x = np.zeros(b.shape)
    # Gucken ob die Matrix invertierbar ist (geht natürlich auch in der
    # Schleife durch überprüfen des Diagonalelementes)
    if (R.flat[0:-1:(len(b)+1)]==0).any():
        print('FEHLER: Die Matrix R ist singulär!\nEs wird ein 0-array zurückgegeben')
        return np.zeros(b.shape)
    x[-1] = b[-1]/R[-1, -1]
    for k in range(len(b)-2, -1, -1):
        x[k] = (b[k]-np.sum(x[k+1:]*R[k, k+1::]))/R[k, k]
    return x

#Tests
#reguläre Matrix erstellen
while True:
    Areg=10*np.random.rand(10,10)-5
    Areg=np.triu(Areg)
    if np.min(np.abs(np.diag(Areg)))>10**-3:
        break
# Lösungsvektor
xex=10*np.random.rand(10)-5
#rechte Seite
breg=Areg@xex
#singuläre Matrix
Asin=10*np.random.rand(10,10)-5
Asin=np.triu(Areg)
Asin[7,7]=0
bsin=Asin@xex

#Lösugnen berechnen
print('Reguläre Matrix:')
xregmy=rueckwaertseinsetzen(Areg,breg)
xregnp=np.linalg.solve(Areg,breg)
print(xregmy)
print(xex)
print(np.linalg.norm(Areg@xregmy-breg))
print(np.linalg.norm(Areg@xregnp-breg))

print('\nSinguläre Matrix:')
xsinmy=rueckwaertseinsetzen(Asin,bsin)
print(xsinmy)
print(np.linalg.norm(Asin@xsinmy-bsin))
#Gibt eine Fehlermeldung
xsinmnp=np.linalg.solve(Asin,bsin)