From 49d42cdd6895ae071ec6835301ce322e4fc91f1b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Sun, 11 Nov 2018 18:49:43 +0100 Subject: [PATCH] initial --- Projekt1_Knuettel_Daniel.c | 127 +++++++++++++++++++++++++++++++++++++ gaussian.py | 75 ++++++++++++++++++++++ 2 files changed, 202 insertions(+) create mode 100644 Projekt1_Knuettel_Daniel.c create mode 100644 gaussian.py diff --git a/Projekt1_Knuettel_Daniel.c b/Projekt1_Knuettel_Daniel.c new file mode 100644 index 0000000..ac34271 --- /dev/null +++ b/Projekt1_Knuettel_Daniel.c @@ -0,0 +1,127 @@ +// Daniel Knüttel Aufgabe2 + +/* +* Copyright(c) 2018 Daniel Knüttel +*/ + +/* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Dieses Programm ist Freie Software: Sie können es unter den Bedingungen +* der GNU General Public License, wie von der Free Software Foundation, +* Version 3 der Lizenz oder (nach Ihrer Wahl) jeder neueren +* veröffentlichten Version, weiterverbreiten und/oder modifizieren. +* +* Dieses Programm wird in der Hoffnung, dass es nützlich sein wird, aber +* OHNE JEDE GEWÄHRLEISTUNG, bereitgestellt; sogar ohne die implizite +* Gewährleistung der MARKTFÄHIGKEIT oder EIGNUNG FÜR EINEN BESTIMMTEN ZWECK. +* Siehe die GNU General Public License für weitere Details. +* +* Sie sollten eine Kopie der GNU General Public License zusammen mit diesem +* Programm erhalten haben. Wenn nicht, siehe . +*/ + +#include +#include + +/* + * \brief Computes the householder partition of A. + * Stores the generating vectors of Q in the memory passed as A, + * the diagonal elements in alpha and the rest of R in the memory + * passed as A. + * + * \param A Matrix to partition. After the partition this memory contains Q and R. + * \param alpha Will be filled with diagonal elements of R. + * \param n Length of A[i]. + * \param m Length of A. + * \return Status. + * */ +int householder(double ** A, double * alpha, int n, int m); + +/* + * \brief Computes the solution of and n times n system using a QR partition. + * Overwrites b. + * + * \param A Matrix that has been modified by the function householder. + * \param alpha Diagonal elements of R, produced by the function householder. + * \param n Length of A and A[i]. + * \param b Vector b for Ax=b. Will be overwritten with elements of x. + * \return Status. + * */ +int solve_QR(double ** A, double * alpha, int n, double * b); + +/* + * \brief Computes Q^Tb and stores the result in b. + * + * \param A Matrix that has been modified by the function householder. + * \param alpha Diagonal elements of R, produced by the function householder. + * \param n Length of A and A[i]. + * \param b Vector b for Ax=b. Will be overwritten with elements of the result. + * + * \return Status. + * */ +int qtb(double ** A, double * alpha, int n, double * b); + +/* + * + * \brief Uses backwards substitution to solve R = Q^Tb. + * + * \param A Matrix that has been modified by the function householder. + * \param alpha Diagonal elements of R, produced by the function householder. + * \param n Length of A and A[i]. + * \param b Vector b for Ax=b. Will be overwritten with elements of the result. + * + * \return Status. + * + * */ +int rw_subs(double ** A, double * alpha, int n, double * b); + + +/* + * \brief Print the matrix in a nicely formatted way to stream. + * + * \param stream The output stream. + * \param matrix The matrix to print. + * \param n Length of matrix[i]. + * \param m Length of matrix. + * + * */ +int fprintm(FILE * stream, double ** matrix, int n, int m); +/* + * \brief Print the vector in a nicely formatted way to stream. + * + * \param stream The output stream. + * \param vector The vector to print. + * \param n Length of vector. + * + * */ +int fprintv(FILE * stream, double * vector, int n); + +#define printm(matrix, n, m) fprintm(stdout, matrix, n, m) +#define printv(vector, n) fprintv(stdout, vector, n) + +int main(void) +{ + int status = 0; + return status; +} + + +int householder(double ** A, double * alpha, int n, int m) +{ + int i, j; + for(i = 0; i < n; i++) + { + // Bring this column in triangular shape. + } +} diff --git a/gaussian.py b/gaussian.py new file mode 100644 index 0000000..bffb692 --- /dev/null +++ b/gaussian.py @@ -0,0 +1,75 @@ +import copy +import pprint +import numpy + +def gaussian_with_pivot(A, b, round = lambda x: x): + A = numpy.array(A, dtype = numpy.float64) + b = numpy.array(b, dtype = numpy.float64) + L = numpy.identity(len(A), dtype = numpy.float64) + s = numpy.zeros(len(A), dtype = numpy.float64) + + for k in range(len(A)): + max_a = numpy.absolute(A[k:,k]).argmax() + A[[k, k + max_a]] = A[[k + max_a, k]] + b[[k, k + max_a]] = b[[k + max_a, k]] + s[k] = k + max_a + + for i in range(k + 1, len(A)): + L[i][k] = round(A[i][k]/A[k][k]) + b[i] = round(b[i] - round(L[i][k] * b[k])) + + for j in range(k + 1, len(A)): + A[i][j] = round(A[i][j] - round(L[i][k] * A[k][j])) + + A[i][k] = 0 + return A, b, L, s + +def gaussian_no_pivot(A, b, round = lambda x: x): + A = numpy.array(A, dtype = numpy.float64) + b = numpy.array(b, dtype = numpy.float64) + L = numpy.identity(len(A), dtype = numpy.float64) + + for k in range(len(A)): + + for i in range(k + 1, len(A)): + L[i][k] = round(A[i][k]/A[k][k]) + b[i] = round(b[i] - round(L[i][k] * b[k])) + + for j in range(k + 1, len(A)): + A[i][j] = round(A[i][j] - round(L[i][k] * A[k][j])) + + A[i][k] = 0 + return A, b, L + + +if( __name__ == "__main__"): + # For Exercise 1 + A = [[10, 1, 0], [15, 2, 1], [5, 0, 1]] + b = [12, 18, 4] + pprint.pprint(gaussian_with_pivot(A, b)) + + + # For Exercise 2 + A = [[7, 1, 1] + , [10, 1, 1] + , [1000, 0, 1]] + b = [10, 13, 1001] + pprint.pprint(A) + + def round4(x): + return float("%.4g" % x) + + pprint.pprint(gaussian_with_pivot(A, b, round4)) + pprint.pprint(gaussian_no_pivot(A, b, round4)) + + + A = numpy.zeros((10, 10)) + + for i in range(10): + if(i - 1 > 0): + A[i][i - 1] = -1 + A[i][i] = 2 + if(i + 1 < 10): + A[i][i + 1] = -1 + + pprint.pprint(gaussian_no_pivot(A, numpy.ones(10)))