// Kompilieren mit:
//
// g++ -I../hdnum/ -o rohrleitungsnetzwerk rohrleitungsnetzwerk.cc
//
// Das setzt voraus, dass ihr Programm in einem Ordner parallel zur
// entpackten HDNum Bibliothek liegt.

#include <iostream>
#include <cstdlib>
#include "hdnum.hh"


// Funktion zum Aufstellen der Matrix
template<class NumberType>
void flussMatrix( hdnum::DenseMatrix<NumberType> &A ) {
    int M( A.rowsize() );
    int N( A.colsize() );
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    // Numerierung wie auf Zettel 3 nur mit 0 beginnend
    // also v_0, ... ,v_(N^2-1)
    // der Referenzknoten v_r hat Druck 0

    // berechnung der kantenlaenge
    int n = floor(sqrt(N+1));

    for(int i = 0; i < N; i++) {
        int edges = 0;
        if ((i+1)%n != 0) { // nicht linker rand
            edges++;
            if (i-1 >= 0) A(i, i-1) = -1; // falls nicht der referenzknoten
        }
        if ((i+2)%n != 0) { // nicht rechter rand
            edges++;
            A(i, i+1) = -1;
        }
        if (i+n < N) { // nicht unterer rand
            edges++;
            A(i, i+n) = -1;
        }
        if (i-n >= -1) { // nicht oberer rand
            edges++;
            if (i-n >= 0) A(i, i-n) = -1; // falls nicht der referenzknoten
        }
        A(i, i) = edges;
    }
}


// Funktion zur Berechnung der Frobenius-Norm einer Matrix
template<class NumberType>
NumberType frobeniusNorm(const hdnum::DenseMatrix<NumberType> &A) {
    // Error checking
    int M(A.rowsize());
    int N(A.colsize());
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    NumberType result=0.0;

    // iteriere ueber alle zeilen und spalten, quadriere die elemente und summiere
    for (int i=0; i < N; i++) {
        for (int j=0; j < N; j++) {
            result += pow(A(i,j), 2);
        }
    }

    // ziehe wurzel aus summe
    return sqrt(result);
}


// Funktion zur Berechnung des betragsgroessten Eigenwertes mit Potenzmethode
template<class NumberType>
NumberType maxEigenwert(const hdnum::DenseMatrix<NumberType> &A) {
    // Error checking
    int M(A.rowsize());
    int N(A.colsize());
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    // start vektor
    hdnum::Vector<NumberType> r(N);
    r[0] = 0;
    // work copy
    hdnum::Vector<NumberType> r_tmp(N);
    // finde vektor mit Ar != 0
    for (int i = 0;;i++) {
        A.mv(r_tmp, r); // r_tmp = Ar
        if (r_tmp.two_norm() > 0) {
            // Ar != 0
            break;
        } else if (i >= 100) {
            // breche ab, wenn kein passender startvektor gefunden werden kann
            HDNUM_ERROR("Kann keinen Vektor mit Ax != 0 finden. Ist A = 0?");
        }
        // Ar = 0, modifiziere start vektor
        r[i % N] += 1;
    }
    // fuehre iterationsschritt maximal 10000 mal aus
    for (int k=0; k<10000; k++) {
        A.mv(r_tmp, r); // r_tmp = Ar
        r_tmp /= r_tmp.two_norm(); // normiere r_tmp
        if ((r-r_tmp).two_norm() == 0)
            // breche ab ab, wenn maximale genauigkeit erreicht wurde
            break;
        r = r_tmp;
    }
    // berechne eigenwert mit rayleigh quotient
    A.mv(r_tmp, r); // r_tmp = Ar
    return (r * r_tmp)/r.two_norm_2();
}

// Hauptprogramm
int main(int argc, char ** argv) {
    // Anzahl der Knoten
    const int N(10);
    std::cout << "Knotenanzahl N: " << N << std::endl;

    // Größe der Matrix
    const int n(N*N-1);

    // Datentyp für die Matrix
    typedef double REAL;

    // Matrix initialisieren
    hdnum::DenseMatrix<REAL> A(n,n);

    // Pretty-printing einmal setzen für alle Matrizen
    A.scientific(false);
    A.width(15);

    flussMatrix(A);
    if (N<=4)
        std::cout << A << std::endl;

    // Bei Schwierigkeiten mit Teilaufgabe a) können Sie Teilaufgaben b)
    // und c) mit folgender Matrix testen
    int size_b = 3;
    hdnum::DenseMatrix<REAL> B(size_b,size_b);
    for (std::size_t i=0; i<size_b; ++i) {
        for (std::size_t j=0; j<size_b; ++j) {
            B[i][j] = i+j;
        }
    }
    std::cout << "Zeilen-Summen-Norm: " << A.norm_infty() << std::endl;
    std::cout << "Spalten-Summen-Norm: " << A.norm_1() << std::endl;
    std::cout << "Frobenius-Norm: " << frobeniusNorm(A) << std::endl;
    std::cout << "Maximaler Eigenwert: " << maxEigenwert(A) << std::endl;

    return 0;
}