gistfile1.cpp

#include <cmath> // sqrt
#include <cstdio>
#include <cstdlib>
#include <stack>
#include <vector>
#include <iostream>

using namespace std;

const unsigned int size = 5;
const unsigned int iterations = 10;

double euclideanNorm(const unsigned int size, double *x)
{
    double sumOfSquares = .0;
    for (unsigned int i = 0; i < size; i++)
        sumOfSquares += x[i] * x[i];

    return sqrt (sumOfSquares);
}

double* normalizeVector(const unsigned int size, double *x)
{
    double norm = euclideanNorm(size, x);
    double *result = new double[size];

    for (unsigned int i = 0; i < size; i++)
        result[i] = x[i] / norm;

    return result;
}

double* multiplicateMatrix(const unsigned int size, double matrix[][size], double *x)
{
    double *result = new double[size];

    for (unsigned int i = 0; i < size; i++)
    {
        result[i] = .0;
        for (unsigned int k = 0; k < size; k++)
            result[i] = matrix[i][k] * x[k];
    }

    return result;
}

void outputEigenvalues(const unsigned int size, double *yActual, double *yPrevious)
{
    for (unsigned int i = 0; i < size; i++)
        fprintf(stdout, "%14.4lf, ", yActual[i] / yPrevious[i]);

    cout << endl;
}

void outputVector(const unsigned int size, double *x)
{
    for (unsigned int i = 0; i < size; i++)
        fprintf(stdout, "%14.4lf ", x[i]);

    cout << endl;
}

int main(void)
{
    /*
     * Variablendeklaration.
     *
     * Speichere nur die letzten beiden Iterationen,
     * alles andere wird nicht benötigt.
     */

    double* xPrevious = new double[size];
    double* yPrevious = new double[size];

    double* xActual = new double[size];
    double* yActual = new double[size];

    /*
     * Matrix und Startvektor.
     */

    double matrix[size][size] = {
        { 5.0,  4.0,  4.0,  5.0,  6.0 },
        { 0.0,  8.0,  5.0,  6.0,  7.0 },
        { 0.0,  0.0,  6.0,  7.0,  8.0 },
        { 0.0,  0.0,  0.0, -4.0,  9.0 },
        { 0.0,  0.0,  0.0,  0.0, -2.0 }
    };

    double startVector[size] = {1.0, 1.0, 1.0, 1.0, 1.0};

    /*
     * Algorithmus: Direkte Vektoriteration.
     */

    xActual = normalizeVector(size, startVector);
    for (int i = 1; i < iterations; i++)
    {
        xPrevious = xActual;
        delete[] xActual;

        yPrevious = yActual;
        delete[] yActual;

        yActual = multiplicateMatrix(size, matrix, xActual);
        xActual = normalizeVector(size, yActual);

        outputVector(size, xActual);
        outputEigenvalues(size, yActual, yPrevious);
        cout << endl;
    }

    // Ausgabe

    //outputEigenvalues(size, yActual, yPrevious);
}