#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Musterlösung für Aufgabe 35

import numpy as np

def QRZerlegung(A):
    ''' Berechnet die QR-Zerlegung der Matrix A aus R^{mxn}
    Input: Matrix A
    Output: Q, R
    '''
    R = A.astype('float') #oder R = A*1.
    m, n = A.shape
    Q = np.eye(m)
    for k in range(n):
        u = R[k:, [k]]
        #eckige Klammern um k sorgen dafür, dass der Vektor u eindeutig eine
        #Spalte bleibt und @ genau das macht, was man erwartet/möchte.
        #Ohne die Klammern wird einer flacher Vektor daraus.
        alpha = -np.sign(R[k, k])*np.linalg.norm(u)
        # A hat den Rang n, d.h. ||v|| ist immer ungleich 0
        if alpha == 0:
            #Ist R[k,k] = 0, ergibt sich alpha = 0, da np.sign(0)=0
            #Den Spezialfall decken wir hier ab.
            alpha = -np.linalg.norm(u)
        e1 = np.eye(m-k, 1)
        u -= alpha*e1
        Qk = np.eye(m)
        Qk[k:, k:] = np.eye(m-k) - (u @ u.T)/(alpha*(alpha - R[k, k]))
        for j in range(k+1, n):
            R[k:, [j]] -= (R[k:, [j]].T @ u)*u / (alpha*(alpha-R[k, k]))
        R[k:, [k]] = alpha*e1
        Q = Q @ Qk
    return Q, R