diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..efed2b5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +# Ignore Mac DS_Store files in all directories +.DS_Store + + +#Ignore EAGLE board and schematic auto save files (.b#x and .s#x) in all directories +*.b#? +*.s#? + + +#Development file for making local changes to Arduino libraries +#Ignore in all directories +.development \ No newline at end of file diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 6a02661..0000000 --- a/.travis.yml +++ /dev/null @@ -1,27 +0,0 @@ -language: c -sudo: false - -# Blacklist -branches: - except: - - gh-pages - -env: - global: - - PRETTYNAME="MatrixMath" -# Optional, will default to "$TRAVIS_BUILD_DIR/Doxyfile" -# - DOXYFILE: $TRAVIS_BUILD_DIR/Doxyfile - -before_install: - - source <(curl -SLs https://raw.githubusercontent.com/adafruit/travis-ci-arduino/master/install.sh) - -install: -# - arduino --install-library "AUnit" - -script: - - build_main_platforms - -# Generate and deploy documentation -after_success: - - source <(curl -SLs https://raw.githubusercontent.com/adafruit/travis-ci-arduino/master/library_check.sh) - - source <(curl -SLs https://raw.githubusercontent.com/adafruit/travis-ci-arduino/master/doxy_gen_and_deploy.sh) diff --git a/MatrixMath.cpp b/MatrixMath.cpp deleted file mode 100644 index 6d7ff6d..0000000 --- a/MatrixMath.cpp +++ /dev/null @@ -1,203 +0,0 @@ -/* - * MatrixMath.cpp Library for Matrix Math - * - * Created by Charlie Matlack on 12/18/10. - * Modified from code by RobH45345 on Arduino Forums, algorithm from - * NUMERICAL RECIPES: The Art of Scientific Computing. - */ - -#include "MatrixMath.h" - -#define NR_END 1 - -MatrixMath Matrix; // Pre-instantiate - -// Matrix Printing Routine -// Uses tabs to separate numbers under assumption printed mtx_type width won't cause problems -void MatrixMath::Print(mtx_type* A, int m, int n, String label) -{ - // A = input matrix (m x n) - int i, j; - Serial.println(); - Serial.println(label); - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - Serial.print(A[n * i + j]); - Serial.print("\t"); - } - Serial.println(); - } -} - -void MatrixMath::Copy(mtx_type* A, int n, int m, mtx_type* B) -{ - int i, j; - for (i = 0; i < m; i++) - for(j = 0; j < n; j++) - { - B[n * i + j] = A[n * i + j]; - } -} - -//Matrix Multiplication Routine -// C = A*B -void MatrixMath::Multiply(mtx_type* A, mtx_type* B, int m, int p, int n, mtx_type* C) -{ - // A = input matrix (m x p) - // B = input matrix (p x n) - // m = number of rows in A - // p = number of columns in A = number of rows in B - // n = number of columns in B - // C = output matrix = A*B (m x n) - int i, j, k; - for (i = 0; i < m; i++) - for(j = 0; j < n; j++) - { - C[n * i + j] = 0; - for (k = 0; k < p; k++) - C[n * i + j] = C[n * i + j] + A[p * i + k] * B[n * k + j]; - } -} - - -//Matrix Addition Routine -void MatrixMath::Add(mtx_type* A, mtx_type* B, int m, int n, mtx_type* C) -{ - // A = input matrix (m x n) - // B = input matrix (m x n) - // m = number of rows in A = number of rows in B - // n = number of columns in A = number of columns in B - // C = output matrix = A+B (m x n) - int i, j; - for (i = 0; i < m; i++) - for(j = 0; j < n; j++) - C[n * i + j] = A[n * i + j] + B[n * i + j]; -} - - -//Matrix Subtraction Routine -void MatrixMath::Subtract(mtx_type* A, mtx_type* B, int m, int n, mtx_type* C) -{ - // A = input matrix (m x n) - // B = input matrix (m x n) - // m = number of rows in A = number of rows in B - // n = number of columns in A = number of columns in B - // C = output matrix = A-B (m x n) - int i, j; - for (i = 0; i < m; i++) - for(j = 0; j < n; j++) - C[n * i + j] = A[n * i + j] - B[n * i + j]; -} - - -//Matrix Transpose Routine -void MatrixMath::Transpose(mtx_type* A, int m, int n, mtx_type* C) -{ - // A = input matrix (m x n) - // m = number of rows in A - // n = number of columns in A - // C = output matrix = the transpose of A (n x m) - int i, j; - for (i = 0; i < m; i++) - for(j = 0; j < n; j++) - C[m * j + i] = A[n * i + j]; -} - -void MatrixMath::Scale(mtx_type* A, int m, int n, mtx_type k) -{ - for (int i = 0; i < m; i++) - for (int j = 0; j < n; j++) - A[n * i + j] = A[n * i + j] * k; -} - - -//Matrix Inversion Routine -// * This function inverts a matrix based on the Gauss Jordan method. -// * Specifically, it uses partial pivoting to improve numeric stability. -// * The algorithm is drawn from those presented in -// NUMERICAL RECIPES: The Art of Scientific Computing. -// * The function returns 1 on success, 0 on failure. -// * NOTE: The argument is ALSO the result matrix, meaning the input matrix is REPLACED -int MatrixMath::Invert(mtx_type* A, int n) -{ - // A = input matrix AND result matrix - // n = number of rows = number of columns in A (n x n) - int pivrow = 0; // keeps track of current pivot row - int k, i, j; // k: overall index along diagonal; i: row index; j: col index - int pivrows[n]; // keeps track of rows swaps to undo at end - mtx_type tmp; // used for finding max value and making column swaps - - for (k = 0; k < n; k++) - { - // find pivot row, the row with biggest entry in current column - tmp = 0; - for (i = k; i < n; i++) - { - if (abs(A[i * n + k]) >= tmp) // 'Avoid using other functions inside abs()?' - { - tmp = abs(A[i * n + k]); - pivrow = i; - } - } - - // check for singular matrix - if (A[pivrow * n + k] == 0.0f) - { - Serial.println("Inversion failed due to singular matrix"); - return 0; - } - - // Execute pivot (row swap) if needed - if (pivrow != k) - { - // swap row k with pivrow - for (j = 0; j < n; j++) - { - tmp = A[k * n + j]; - A[k * n + j] = A[pivrow * n + j]; - A[pivrow * n + j] = tmp; - } - } - pivrows[k] = pivrow; // record row swap (even if no swap happened) - - tmp = 1.0f / A[k * n + k]; // invert pivot element - A[k * n + k] = 1.0f; // This element of input matrix becomes result matrix - - // Perform row reduction (divide every element by pivot) - for (j = 0; j < n; j++) - { - A[k * n + j] = A[k * n + j] * tmp; - } - - // Now eliminate all other entries in this column - for (i = 0; i < n; i++) - { - if (i != k) - { - tmp = A[i * n + k]; - A[i * n + k] = 0.0f; // The other place where in matrix becomes result mat - for (j = 0; j < n; j++) - { - A[i * n + j] = A[i * n + j] - A[k * n + j] * tmp; - } - } - } - } - - // Done, now need to undo pivot row swaps by doing column swaps in reverse order - for (k = n - 1; k >= 0; k--) - { - if (pivrows[k] != k) - { - for (i = 0; i < n; i++) - { - tmp = A[i * n + k]; - A[i * n + k] = A[i * n + pivrows[k]]; - A[i * n + pivrows[k]] = tmp; - } - } - } - return 1; -} diff --git a/MatrixMath.h b/MatrixMath.h index 08234ad..b8c32c8 100644 --- a/MatrixMath.h +++ b/MatrixMath.h @@ -1,38 +1,298 @@ /* - * MatrixMath.h Library for Matrix Math - * - * Created by Charlie Matlack on 12/18/10. - * Modified from code by RobH45345 on Arduino Forums, algorithm from - * NUMERICAL RECIPES: The Art of Scientific Computing. - * Modified to work with Arduino 1.0/1.5 by randomvibe & robtillaart - * Made into a real library on GitHub by Vasilis Georgitzikis (tzikis) - * so that it's easy to use and install (March 2015) + + FILENAME: MatrixMath.h + + + MatrixMath.h Library for Matrix Math + Created by Charlie Matlack on 12/18/10. + Modified from code by RobH45345 on Arduino Forums, algorithm from + NUMERICAL RECIPES: The Art of Scientific Computing. + Modified to work with Arduino 1.0/1.5 by randomvibe & robtillaart + Made into a real library on GitHub by Vasilis Georgitzikis (tzikis) + so that it's easy to use and install (March 2015) + + 2020/08/31:1444 (UTC-5) -- Author: Orlando S. Hoilett + - moved to a templated class to make the class more usable across + different data types. This also allows the user to avoid floats and + doubles in order to save space and increase processing speed + - changed function names to have first letter be lower case instead + of upper cass + - added a print function without the "String label" parameter + - added a scale function that edits doesn't modify input array + - added const where sensible + */ + #ifndef MatrixMath_h #define MatrixMath_h -#define mtx_type double - #if defined(ARDUINO) && ARDUINO >= 100 #include "Arduino.h" #else #include "WProgram.h" #endif + +//Ensures print function will work across different architectures. Some Arduino +//variants use SerialUSB, not Serial, for printing to the Serial Monitor/Plotter. +#if defined(ARDUINO_ARCH_SAMD) + #define MatrixMathSerial SerialUSB +#else + #define MatrixMathSerial Serial +#endif + + +template + class MatrixMath { + public: //MatrixMath(); - void Print(mtx_type* A, int m, int n, String label); - void Copy(mtx_type* A, int n, int m, mtx_type* B); - void Multiply(mtx_type* A, mtx_type* B, int m, int p, int n, mtx_type* C); - void Add(mtx_type* A, mtx_type* B, int m, int n, mtx_type* C); - void Subtract(mtx_type* A, mtx_type* B, int m, int n, mtx_type* C); - void Transpose(mtx_type* A, int m, int n, mtx_type* C); - void Scale(mtx_type* A, int m, int n, mtx_type k); - int Invert(mtx_type* A, int n); + static void print(const Type* A, int m, int n); + static void print(const Type* A, int m, int n, String label); + + static void copy(const Type* A, int n, int m, Type* B); + + static void multiply(const Type* A, const Type* B, int m, int p, int n, Type* C); + static void add(const Type* A, const Type* B, int m, int n, Type* C); + static void subtract(const Type* A, const Type* B, int m, int n, Type* C); + + static void transpose(const Type* A, int m, int n, Type* C); + static void scale(Type* A, int m, int n, Type k); + static void scale(const Type* A, int m, int n, Type k, Type* C); + + static int invert(float* A, int n); + }; -extern MatrixMath Matrix; + +// Matrix Printing Routine +// Uses tabs to separate numbers under assumption printed float width won't cause problems +template +void MatrixMath::print(const Type* A, int m, int n) +{ + // A = input matrix (m x n) + int i, j; + + for (i = 0; i < m; i++) + { + for (j = 0; j < n; j++) + { + MatrixMathSerial.print(A[n * i + j]); + MatrixMathSerial.print("\t"); + } + MatrixMathSerial.println(); + } +} + +// Matrix Printing Routine +// Uses tabs to separate numbers under assumption printed float width won't cause problems +template +void MatrixMath::print(const Type* A, int m, int n, String label) +{ + // A = input matrix (m x n) + int i, j; + + MatrixMathSerial.println(label); + for (i = 0; i < m; i++) + { + for (j = 0; j < n; j++) + { + MatrixMathSerial.print(A[n * i + j]); + MatrixMathSerial.print("\t"); + } + MatrixMathSerial.println(); + } +} + +template +void MatrixMath::copy(const Type* A, int n, int m, Type* B) +{ + int i, j; + for (i = 0; i < m; i++) + for(j = 0; j < n; j++) + { + B[n * i + j] = A[n * i + j]; + } +} + +//Matrix Multiplication Routine +// C = A*B +template +void MatrixMath::multiply(const Type* A, const Type* B, int m, int p, int n, Type* C) +{ + // A = input matrix (m x p) + // B = input matrix (p x n) + // m = number of rows in A + // p = number of columns in A = number of rows in B + // n = number of columns in B + // C = output matrix = A*B (m x n) + int i, j, k; + for (i = 0; i < m; i++) + for(j = 0; j < n; j++) + { + C[n * i + j] = 0; + for (k = 0; k < p; k++) + C[n * i + j] = C[n * i + j] + A[p * i + k] * B[n * k + j]; + } +} + + +//Matrix Addition Routine +template +void MatrixMath::add(const Type* A, const Type* B, int m, int n, Type* C) +{ + // A = input matrix (m x n) + // B = input matrix (m x n) + // m = number of rows in A = number of rows in B + // n = number of columns in A = number of columns in B + // C = output matrix = A+B (m x n) + int i, j; + for (i = 0; i < m; i++) + for(j = 0; j < n; j++) + C[n * i + j] = A[n * i + j] + B[n * i + j]; +} + + +//Matrix Subtraction Routine +template +void MatrixMath::subtract(const Type* A, const Type* B, int m, int n, Type* C) +{ + // A = input matrix (m x n) + // B = input matrix (m x n) + // m = number of rows in A = number of rows in B + // n = number of columns in A = number of columns in B + // C = output matrix = A-B (m x n) + int i, j; + for (i = 0; i < m; i++) + for(j = 0; j < n; j++) + C[n * i + j] = A[n * i + j] - B[n * i + j]; +} + + +//Matrix Transpose Routine +template +void MatrixMath::transpose(const Type* A, int m, int n, Type* C) +{ + // A = input matrix (m x n) + // m = number of rows in A + // n = number of columns in A + // C = output matrix = the transpose of A (n x m) + int i, j; + for (i = 0; i < m; i++) + for(j = 0; j < n; j++) + C[m * j + i] = A[n * i + j]; +} + +template +void MatrixMath::scale(Type* A, int m, int n, Type k) +{ + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + A[n * i + j] = A[n * i + j] * k; +} + + +template +void MatrixMath::scale(const Type* A, int m, int n, Type k, Type* C) +{ + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + C[n * i + j] = A[n * i + j] * k; +} + + +//Matrix Inversion Routine +// * This function inverts a matrix based on the Gauss Jordan method. +// * Specifically, it uses partial pivoting to improve numeric stability. +// * The algorithm is drawn from those presented in +// NUMERICAL RECIPES: The Art of Scientific Computing. +// * The function returns 1 on success, 0 on failure. +// * NOTE: The argument is ALSO the result matrix, meaning the input matrix is REPLACED +template +int MatrixMath::invert(float* A, int n) +{ + // A = input matrix AND result matrix + // n = number of rows = number of columns in A (n x n) + int pivrow; // keeps track of current pivot row + int k, i, j; // k: overall index along diagonal; i: row index; j: col index + int pivrows[n]; // keeps track of rows swaps to undo at end + float tmp; // used for finding max value and making column swaps + + for (k = 0; k < n; k++) + { + // find pivot row, the row with biggest entry in current column + tmp = 0; + for (i = k; i < n; i++) + { + if (abs(A[i * n + k]) >= tmp) // 'Avoid using other functions inside abs()?' + { + tmp = abs(A[i * n + k]); + pivrow = i; + } + } + + // check for singular matrix + if (A[pivrow * n + k] == 0.0f) + { + Serial.println("Inversion failed due to singular matrix"); + return 0; + } + + // Execute pivot (row swap) if needed + if (pivrow != k) + { + // swap row k with pivrow + for (j = 0; j < n; j++) + { + tmp = A[k * n + j]; + A[k * n + j] = A[pivrow * n + j]; + A[pivrow * n + j] = tmp; + } + } + pivrows[k] = pivrow; // record row swap (even if no swap happened) + + tmp = 1.0f / A[k * n + k]; // invert pivot element + A[k * n + k] = 1.0f; // This element of input matrix becomes result matrix + + // Perform row reduction (divide every element by pivot) + for (j = 0; j < n; j++) + { + A[k * n + j] = A[k * n + j] * tmp; + } + + // Now eliminate all other entries in this column + for (i = 0; i < n; i++) + { + if (i != k) + { + tmp = A[i * n + k]; + A[i * n + k] = 0.0f; // The other place where in matrix becomes result mat + for (j = 0; j < n; j++) + { + A[i * n + j] = A[i * n + j] - A[k * n + j] * tmp; + } + } + } + } + + // Done, now need to undo pivot row swaps by doing column swaps in reverse order + for (k = n - 1; k >= 0; k--) + { + if (pivrows[k] != k) + { + for (i = 0; i < n; i++) + { + tmp = A[i * n + k]; + A[i * n + k] = A[i * n + pivrows[k]]; + A[i * n + pivrows[k]] = tmp; + } + } + } + return 1; +} + + #endif diff --git a/examples/EXAMPLE01_basic/EXAMPLE01_basic.ino b/examples/EXAMPLE01_basic/EXAMPLE01_basic.ino new file mode 100644 index 0000000..8c2c9fb --- /dev/null +++ b/examples/EXAMPLE01_basic/EXAMPLE01_basic.ino @@ -0,0 +1,100 @@ +/* + * FILENAME: EXAMPLE01_basic.ino + * AUTHOR: Orlando S. Hoilett + * CONTACT: orlandohoilett@gmail.com + * VERSION: 1.0.0 + * + * + * AFFILIATIONS + * Linnes Lab, Weldon School of Biomedical Engineering, + * Purdue University, West Lafayette, IN 47907 + * + * + * DESCRIPTION + * Basic test of the MatrixMath library. + * + * + * UPDATES + * Version 1.0.0 + * 2020/08/31:1411> + * - Initialized + * + */ + + +#include + + +const uint16_t N = 2; +const uint16_t M = 3; + +uint16_t A[N][N]; +uint16_t B[N][N]; +uint16_t C[N][N]; + +uint16_t D[M][M]; +uint16_t E[M][M]; + + +void setup() +{ + Serial.begin(9600); + randomSeed(analogRead(A0)); //initialize random number generator + + + //Initialize matrix "A" + A[0][0] = 1; + A[0][1] = 2; + A[1][0] = 3; + A[1][1] = 4; + Serial.println("Initializing matrix, A..."); + MatrixMath::print((uint16_t*)A, N, N, "A ="); + Serial.println(); + + + //Scale matrix "A" + Serial.println("Scaling matrix, A..."); + MatrixMath::scale((uint16_t*)A, N, N, 3); + MatrixMath::print((uint16_t*)A, N, N, "A = "); + Serial.println(); + + + //Transpose matrix "A" + Serial.println("Transposing matrix, A into matrix, C..."); + MatrixMath::transpose((uint16_t*)A, N, N, (uint16_t*)C); + MatrixMath::print((uint16_t*)C, N, N, "C ="); + Serial.println(); + + + //Initialize matrix "D" + for(uint16_t i = 0; i < M; i++) + { + for(uint16_t j = 0; j < M; j++) + { + D[i][j] = random(0,10); + } + } + Serial.println("Initializing matrix, D..."); + MatrixMath::print((uint16_t*)D, N, N, "D = "); + Serial.println(); + + + //Transposing matrix, D into matrix, E + Serial.println("Transposing matrix, D into matrix, E..."); + MatrixMath::transpose((uint16_t*)D, N, N, (uint16_t*)E); + MatrixMath::print((uint16_t*)E, N, N, "E ="); + Serial.println(); + + + //Scaling matrix, A into matrix, B + Serial.println("Scale matrix, A into matrix, B..."); + MatrixMath::scale((uint16_t*)A, N, N, 10, (uint16_t*)B); + MatrixMath::print((uint16_t*)A, N, N, "A = "); + Serial.println(); + MatrixMath::print((uint16_t*)B, N, N, "B = "); + Serial.println(); +} + +void loop() +{ +} diff --git a/examples/MatrixMathExample/.due.test.skip b/examples/MatrixMathExample/.due.test.skip deleted file mode 100644 index e69de29..0000000 diff --git a/examples/MatrixMathExample/MatrixMathExample.ino b/examples/MatrixMathExample/MatrixMathExample.ino deleted file mode 100644 index ca82e6d..0000000 --- a/examples/MatrixMathExample/MatrixMathExample.ino +++ /dev/null @@ -1,76 +0,0 @@ -#include - - -#define N (2) - -mtx_type A[N][N]; -mtx_type B[N][N]; -mtx_type C[N][N]; -mtx_type v[N]; // This is a row vector -mtx_type w[N]; - -mtx_type maxVal = 10; // maxValimum random matrix entry range - -void setup() -{ - Serial.begin(9600); - - // Initialize matrices - for (int i = 0; i < N; i++) - { - v[i] = i + 1; // vector of sequential numbers - for (int j = 0; j < N; j++) - { - A[i][j] = random(maxVal) - maxVal / 2.0f; // A is random - if (i == j) - { - B[i][j] = 1.0f; // B is identity - } - else - { - B[i][j] = 0.0f; - } - } - } - -} - -void loop() -{ - - Matrix.Multiply((mtx_type*)A, (mtx_type*)B, N, N, N, (mtx_type*)C); - - Serial.println("\nAfter multiplying C = A*B:"); - Matrix.Print((mtx_type*)A, N, N, "A"); - - Matrix.Print((mtx_type*)B, N, N, "B"); - Matrix.Print((mtx_type*)C, N, N, "C"); - Matrix.Print((mtx_type*)v, N, 1, "v"); - - Matrix.Add((mtx_type*) B, (mtx_type*) C, N, N, (mtx_type*) C); - Serial.println("\nC = B+C (addition in-place)"); - Matrix.Print((mtx_type*)C, N, N, "C"); - Matrix.Print((mtx_type*)B, N, N, "B"); - - Matrix.Copy((mtx_type*)A, N, N, (mtx_type*)B); - Serial.println("\nCopied A to B:"); - Matrix.Print((mtx_type*)B, N, N, "B"); - - Matrix.Invert((mtx_type*)A, N); - Serial.println("\nInverted A:"); - Matrix.Print((mtx_type*)A, N, N, "A"); - - Matrix.Multiply((mtx_type*)A, (mtx_type*)B, N, N, N, (mtx_type*)C); - Serial.println("\nC = A*B"); - Matrix.Print((mtx_type*)C, N, N, "C"); - - // Because the library uses pointers and DIY indexing, - // a 1D vector can be smoothly handled as either a row or col vector - // depending on the dimensions we specify when calling a function - Matrix.Multiply((mtx_type*)C, (mtx_type*)v, N, N, 1, (mtx_type*)w); - Serial.println("\n C*v = w:"); - Matrix.Print((mtx_type*)v, N, 1, "v"); - Matrix.Print((mtx_type*)w, N, 1, "w"); - - while(1); -} diff --git a/keywords.txt b/keywords.txt new file mode 100644 index 0000000..f143a4d --- /dev/null +++ b/keywords.txt @@ -0,0 +1,30 @@ +####################################### +# Syntax Coloring Map for MatrixMath +####################################### + +####################################### +# Datatypes (KEYWORD1) +####################################### + +MatrixMath KEYWORD1 + + + +####################################### +# Methods and Functions (KEYWORD2) +####################################### + +print KEYWORD2 +copy KEYWORD2 +multiply KEYWORD2 +add KEYWORD2 +subtract KEYWORD2 +transpose KEYWORD2 +scale KEYWORD2 +invert KEYWORD2 + + + +####################################### +# Constants (LITERAL1) +####################################### \ No newline at end of file diff --git a/library.properties b/library.properties index de6bd4d..05e60cf 100644 --- a/library.properties +++ b/library.properties @@ -1,9 +1,9 @@ name=MatrixMath -version=1.0 +version=2.0.0 author=Charlie Matlack maintainer=Charlie Matlack sentence=Minimal linear algebra library -paragraph=A minimal linear algebra library for Arduino. This gives you all the basics in a lean package, up to in-place matrix inversion. Matrices are represented as simple 2D arrays, so you need to check dimension agreement manually. A far more capable, testable, and friendly linear algebra library for Arduino is BasicLinearAlgebra +paragraph=A minimal linear algebra library for Arduino. This gives you all the basics in a lean package, up to in-place matrix inversion. Matrices are represented as simple 2D arrays, so you need to check dimension agreement manually. category=Data Processing url=https://github.com/eecharlie/MatrixMath architectures=*