24

LinearSystemsToolkit/AllTest.c

@@ -1,5 +1,8 @@
 #include "VectorTest.h"
+#include "MatrixTest.h"
 
 int main() {
     vectorAllTests();
+    matrixAllTests();
+    return 0;
 }

LinearSystemsToolkit/Matrix.c

@@ -0,0 +1,99 @@
+#include <malloc.h>
+#include <stdio.h>
+#include "Matrix.h"
+
+Matrix* matrixCreate(int n) {
+    Matrix* M = malloc(sizeof(Matrix));
+    M->n = n;
+    double** A = malloc(n * sizeof(double*) + n * n * sizeof(double));
+	char* p = A;
+	p += n * sizeof(double*);
+	for (int i = 0; i < n; i++) {
+		A[i] = p + i * n * sizeof(double);
+	}
+	for (int i = 0; i < n; i++) {
+		for (int j = 0; j < n; j++) {
+			A[i][j] = 0;
+		}
+	}
+	M->data = A;
+	return M;
+}
+
+void matrixDestroy(Matrix* M) {
+	free(M->data);
+	free(M);
+}
+
+Matrix* matrixIdentity(int n) {
+    Matrix* M = matrixCreate(n);
+    for (int i = 0; i < n; ++i) {
+        M->data[i][i] = 1;
+    }
+    return M;
+}
+
+Matrix* matrixReadFromFile(const char* name, int n) {
+    FILE* f = fopen(name, "r");
+    Matrix* A = matrixCreate(n);
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < n; ++j) {
+            fscanf(f, "%lf", &(A->data[i][j]));
+        }
+    }
+    fclose(f);
+    return A;
+}
+
+void matrixPrint(Matrix* M) {
+    for (int i = 0; i < M->n; ++i) {
+        for (int j = 0; j < M->n; ++j) {
+            printf("%lf ", M->data[i][j]);
+        }
+        printf("\n");
+    }
+}
+
+Matrix* matrixSum(Matrix* left, Matrix* right) {
+	Matrix* S = matrixCreate(left->n);
+	for (int i = 0; i < left->n; ++i) {
+		for (int j = 0; j < right->n; ++j) {
+			S->data[i][j] = left->data[i][j] + right->data[i][j];
+		}
+	}
+	return S;
+}
+
+Matrix* matrixDiff(Matrix* left, Matrix* right) {
+	Matrix* D = matrixCreate(left->n);
+	for (int i = 0; i < left->n; ++i) {
+		for (int j = 0; j < right->n; ++j) {
+			D->data[i][j] = left->data[i][j] - right->data[i][j];
+		}
+	}
+	return D;
+}
+
+Matrix* matrixScale(Matrix* A, double num) {
+	Matrix* S = matrixCreate(A->n);
+	for (int i = 0; i < A->n; ++i) {
+		for (int j = 0; j < A->n; ++j) {
+			S->data[i][j] = A->data[i][j] * num;
+		}
+	}
+	return S;
+}
+
+int matrixEqual(Matrix* left, Matrix* right) {
+	if (left->n != right->n) {
+		return 0;
+	}
+	for (int i = 0; i < left->n; ++i) {
+		for (int j = 0; j < left->n; ++j) {
+			if (left->data[i][j] != right->data[i][j]) {
+				return 0;
+			}
+		}
+	}
+	return 1;
+}

LinearSystemsToolkit/Matrix.h

@@ -0,0 +1,141 @@
+#pragma once
+
+#include "Vector.h"
+
+/**
+ * @brief A matrix of regular values.
+ * 
+ * This struct represents a square matrix.
+*/
+typedef struct Matrix {
+	double** data; // matrix values
+	int      n;    // matrix size
+} Matrix;
+
+/**
+ * @brief Initialize a matrix.
+ * 
+ * @param n The size of the matrix.
+ * 
+ * @return A pointer to the matrix.
+*/
+Matrix* matrixCreate(int n);
+
+/**
+ * @brief Free a matrix.
+ * 
+ * @param A The matrix to free.
+*/
+void matrixDestroy(Matrix* A);
+
+/**
+ * @brief Create an identity matrix.
+ * 
+ * @param n The size of the matrix.
+ * 
+ * @return A pointer to the matrix.
+*/
+Matrix* matrixIdentity(int n);
+
+/**
+ * @brief Read a matrix from a file.
+ * 
+ * @param name The name of the file.
+ * @param n    The size of the matrix.
+ * 
+ * @return A pointer to the matrix.
+*/
+Matrix* matrixReadFromFile(const char* name, int n);
+
+/**
+ * @brief Print a matrix to console.
+ * 
+ * @param A The matrix to print.
+*/
+void matrixPrint(Matrix* A);
+
+/**
+ * @brief Sums two matrices.
+ * 
+ * @param left  Left term.
+ * @param right Right term.
+ * 
+ * @return The sum.
+*/
+Matrix* matrixSum(Matrix* left, Matrix* right);
+
+/**
+ * @brief Subtracts two matrices.
+ * 
+ * @param left  Left term.
+ * @param right Right term.
+ * 
+ * @return The difference.
+*/
+Matrix* matrixDiff(Matrix* left, Matrix* right);
+
+/**
+ * @brief Scales a matrix by some factor.
+ * 
+ * @param A   The matrix to scale.
+ * @param num The factor by which the matrix is scaled.
+ * 
+ * @return The scaled matrix.
+*/
+Matrix* matrixScale(Matrix* A, double num);
+
+/**
+ * @brief Multiplies a matrix by a vector.
+ * 
+ * Matrix and vector dimensions must match. The result of such multiplication 
+ * is b = A * vec.
+ * 
+ * @param A   The matrix to multiply.
+ * @param vec The vector to multiply.
+ * 
+ * @return A vector containing the product.
+*/
+Vector* matrixMultiplyByVector(Matrix* A, Vector* vec);
+
+/**
+ * @brief Multiplies two matrices.
+ * 
+ * @param left  The left term.
+ * @param right The right term.
+ * 
+ * @return The product of the two matrices' multiplication.
+*/
+Matrix* matrixMultiply(Matrix* left, Matrix* right);
+
+/**
+ * @brief Creates a normalized orthogonal matrix.
+ * 
+ * Orthogonal matrix is calculated according to the formula:
+ * Q = I - 2 * w * w' / (w' * w), where
+ * I is an identity matrix,
+ * w is a given vector seed
+ * 
+ * @param seed The seed vector.
+ * 
+ * @return A normalized orthogonal matrix.
+*/
+Matrix* matrixOrthogonal(Vector* seed);
+
+/**
+ * @brief Creates a diagonal matrix from a vector.
+ * 
+ * @param diag Vector containing the diagonal elements.
+ * 
+ * @return A diagonal matrix.
+*/
+Matrix* matrixDiagonal(Vector* diag);
+
+/**
+ * @brief Checks if two matrices are equal.
+ * 
+ * @param left  The left term.
+ * @param right The right term.
+ * 
+ * @return 1 if the two matrices are equal, 0 otherwise.
+*/
+int matrixEqual(Matrix* left, Matrix* right);

LinearSystemsToolkit/MatrixTest.c

@@ -0,0 +1,81 @@
+#include <assert.h>
+#include "MatrixTest.h"
+#include "Matrix.h"
+
+/**
+ * matrixIdentity can return an identity matrix.
+*/
+void matrixIdentityReturnsIdentityMatrix() {
+    int size  = 3;
+    Matrix* A = matrixIdentity(size);
+    for (int i = 0; i < size; ++i) {
+        assert(1 == A->data[i][i]);
+    }
+    matrixDestroy(A);
+}
+
+/**
+ * matrixSum can return the sum of two matrices.
+*/
+void matrixSumSumsTwoMatrices() {
+    int size  = 3;
+    Matrix* A = matrixCreate(size);
+    Matrix* B = matrixCreate(size);
+    Matrix* S = matrixCreate(size);
+    for (int i = 0; i < size; ++i) {
+        for (int j = 0; j < size; ++j) {
+            A->data[i][j] = i;
+            B->data[i][j] = j;
+            S->data[i][j] = i + j;
+        }
+    }
+    assert(1 == matrixEqual(matrixSum(A, B), S));
+    matrixDestroy(A);
+    matrixDestroy(B);
+    matrixDestroy(S);
+}
+
+/**
+ * matrixDiff can compute the difference of two matrices.
+*/
+void matrixDiffComputesDifferenceOfTwoMatrices() {
+    int size  = 3;
+    Matrix* A = matrixCreate(size);
+    Matrix* B = matrixCreate(size);
+    Matrix* D = matrixCreate(size);
+    for (int i = 0; i < size; ++i) {
+        for (int j = 0; j < size; ++j) {
+            A->data[i][j] = i + j;
+            B->data[i][j] = j;
+            D->data[i][j] = i;
+        }
+    }
+    assert(1 == matrixEqual(matrixDiff(A, B), D));
+    matrixDestroy(A);
+    matrixDestroy(B);
+    matrixDestroy(D);
+}
+
+/**
+ * matrixScale can scale a matrix by some factor.
+*/
+void matrixScaleScalesMatrixBySomeFactor() {
+    int size  = 3;
+    double factor = 3;
+    Matrix* A = matrixCreate(size);
+    Matrix* S = matrixCreate(size);
+    for (int i = 0; i < size; ++i) {
+        for (int j = 0; j < size; ++j) {
+            A->data[i][j] = i + j;
+            S->data[i][j] = (i + j) * factor;
+        }
+    }
+    assert(1 == matrixEqual(matrixScale(A, factor), S));
+    matrixDestroy(A);
+    matrixDestroy(S);
+}
+
+void matrixAllTests() {
+    matrixIdentityReturnsIdentityMatrix();
+    matrixSumSumsTwoMatrices();
+}

LinearSystemsToolkit/MatrixTest.h

@@ -0,0 +1,3 @@
+#pragma once
+
+void matrixAllTests();