From 35b5b79cca17f3191edf1acc7062d98a5a4abec0 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Mon, 18 Oct 2021 10:18:45 -0500 Subject: [PATCH] converting to TIRA --- tira/Matrix.h | 165 --------------------------------------------------------------------------------------------------------------------------------------------------------------------- 1 file changed, 0 insertions(+), 165 deletions(-) delete mode 100644 tira/Matrix.h diff --git a/tira/Matrix.h b/tira/Matrix.h deleted file mode 100644 index f154f7c..0000000 --- a/tira/Matrix.h +++ /dev/null @@ -1,165 +0,0 @@ -#include -#include -#include - -#ifndef LAPACK_COMPLEX_CUSTOM -#define LAPACK_COMPLEX_CUSTOM -#define lapack_complex_float std::complex -#define lapack_complex_double std::complex -#include "lapacke.h" -#endif - -namespace tira { - - template - class Matrix { - type* m_ptr; - int m_rows, m_cols; - - void allocate(int rows, int cols) { - if (m_ptr != NULL) - free(m_ptr); - m_rows = rows; - m_cols = cols; - m_ptr = (type*)malloc(m_rows * m_cols * sizeof(type)); - } - - public: - /// Default constructor, creates an empty matrix (generally filled as a parameter) - Matrix() { - m_rows = m_cols = 0; - m_ptr = NULL; - } - - /// Contructor, creates a rows x cols matrix with undefined values - Matrix(int rows, int cols) : Matrix() { - allocate(rows, cols); - } - - /// De-allocate data and clear the matrix for garbage collection - void del() { - m_rows = 0; - m_cols = 0; - free(m_ptr); - m_ptr = NULL; - } - - void create(int rows, int cols) { - allocate(rows, cols); - } - - int rows() { return m_rows; } - int cols() { return m_cols; } - type* get_ptr() { return m_ptr; } - - void force_square() { - if (m_rows != m_cols) - throw "Matrix must be square!"; - if (m_rows == 0 || m_cols == 0) - throw "Matrix undefined!"; - } - - inline void set(int row, int col, type val) { - m_ptr[col * m_rows + row] = val; - } - - inline type get(int row, int col) { - return m_ptr[col * m_rows + row]; - } - - /// Return the number of bytes in the matrix array - size_t bytes() { return m_rows * m_cols * sizeof(type); } - - /// Create a deep copy of the current matrix and return it - Matrix copy() { - Matrix result(m_rows, m_cols); - memcpy(result.m_ptr, m_ptr, bytes()); - return result; - } - - /// Calculate the determinant of the matrix - type det() { - force_square(); //throw an exception if the matrix isn't square - Matrix tmp = copy(); //copy the current matrix to create space for LUP decomposition - Matrix idx = getrf(tmp); //perform LUP decomposition - - // multiply all elements along the diagonal - type determinant = get(0, 0); //initialize the determinant - for (int i = 1; i < m_rows; i++) //for each remaining element of the diagonal - determinant *= tmp.get(i, i); //calculate the running factor - tmp.del(); //delete memory for the temporary matrix - - // calculate the number of row swaps - int swaps = 0; //initialize the number of row swaps - for (int i = 0; i < m_rows; i++) //iterate through each element of the swap vector - if (idx.get(i, 0) != i + 1) //if the element does not equal the row number - swaps++; //a row was swapped, so increment the swap counter - if (swaps % 2 != 0) //if there was an odd number of row swaps - determinant *= (-1); //multiply the determinant by -1 - return determinant; //return the determinant - } - - /// Calculate the trace of the matrix - type trace() { - force_square(); //throw an exception if the matrix isn't square - type trace = 0; //initialize the trace summation - for (int i = 0; i < m_rows; i++) //for each element along the diagonal - trace += get(i, i); //sum the trace - return trace; //return the matrix trace result - } - - /// Output the matrix to a CSV file - void csv(std::string filename) { - std::ofstream outfile(filename.c_str()); //open the file for writing - for (int r = 0; r < m_rows; r++) { //for each row - for (int c = 0; c < m_cols; c++) { //for each column - outfile << get(r, c); //output the value to a file - if (c != m_cols - 1) outfile << ","; //if this isn't the last column, add a comma - } - if (r != m_rows - 1) outfile << std::endl; //if this isn't the last row, end the line - } - } - - /* Static functions used to generate specialized matrices */ - - /// Static function for generating an identity matrix - static Matrix I(int dim) { - Matrix result(dim, dim); //create and allocate space for the matrix - memset(result.m_ptr, 0, dim * dim * sizeof(type)); //set all elements to zero - for (int i = 0; i < dim; i++) //for each element on the diagonal - result.set(i, i, (type)1); //set that element to 1 - return result; //return the completed matrix - } - - /// Static function for generating a geometric matrix for robustness testing - static Matrix geometric_matrix(int dim) { - Matrix result(dim, dim); //create and allocate space for the matrix - for (int r = 0; r < dim; r++) //for each row - for (int c = 0; c < dim; c++) //for each column - result.set(r, c, pow((type)(2 + r), c)); //set the geometric coefficient - return result; //return the completed matrix - } - - /// Static function analytically generates a solution (right hand side) vector for a geometric matrix - static Matrix geometric_vector(int dim) { - Matrix result(dim, 1); //generate and allocate space for the vector - for (int r = 0; r < dim; r++) - result.set(r, 0, (pow((type)(2 + r), dim) - 1) / (type)(r + 1)); - return result; //return the completed vector - } - - /* Static functions used for LAPACK and BLAS */ - - /// Calls the appropriate LAPACK GETRF function for performing LUP decomposition - static Matrix getrf(Matrix M) { - Matrix idx(M.m_rows, 1); - LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr()); - return idx; - } - static Matrix getrf(Matrix M) { - Matrix idx(M.m_rows, 1); - LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr()); - return idx; - } - }; -} \ No newline at end of file -- libgit2 0.21.4