codebase / stimlib

Commit 35b5b79cca17f3191edf1acc7062d98a5a4abec0

Authored by David Mayerich
1 parent ed60ee52

converting to TIRA

Showing 1 changed file with 0 additions and 165 deletions   Show diff stats
Inline Side-by-side
tira/Matrix.h deleted
View file @ed60ee5
1   -#include <iostream>
2   -#include <fstream>
3   -#include <complex>
4   -
5   -#ifndef LAPACK_COMPLEX_CUSTOM
6   -#define LAPACK_COMPLEX_CUSTOM
7   -#define lapack_complex_float std::complex<float>
8   -#define lapack_complex_double std::complex<double>
9   -#include "lapacke.h"
10   -#endif
11   -
12   -namespace tira {
13   -
14   - template <typename type>
15   - class Matrix {
16   - type* m_ptr;
17   - int m_rows, m_cols;
18   -
19   - void allocate(int rows, int cols) {
20   - if (m_ptr != NULL)
21   - free(m_ptr);
22   - m_rows = rows;
23   - m_cols = cols;
24   - m_ptr = (type*)malloc(m_rows * m_cols * sizeof(type));
25   - }
26   -
27   - public:
28   - /// Default constructor, creates an empty matrix (generally filled as a parameter)
29   - Matrix() {
30   - m_rows = m_cols = 0;
31   - m_ptr = NULL;
32   - }
33   -
34   - /// Contructor, creates a rows x cols matrix with undefined values
35   - Matrix(int rows, int cols) : Matrix() {
36   - allocate(rows, cols);
37   - }
38   -
39   - /// De-allocate data and clear the matrix for garbage collection
40   - void del() {
41   - m_rows = 0;
42   - m_cols = 0;
43   - free(m_ptr);
44   - m_ptr = NULL;
45   - }
46   -
47   - void create(int rows, int cols) {
48   - allocate(rows, cols);
49   - }
50   -
51   - int rows() { return m_rows; }
52   - int cols() { return m_cols; }
53   - type* get_ptr() { return m_ptr; }
54   -
55   - void force_square() {
56   - if (m_rows != m_cols)
57   - throw "Matrix must be square!";
58   - if (m_rows == 0 || m_cols == 0)
59   - throw "Matrix undefined!";
60   - }
61   -
62   - inline void set(int row, int col, type val) {
63   - m_ptr[col * m_rows + row] = val;
64   - }
65   -
66   - inline type get(int row, int col) {
67   - return m_ptr[col * m_rows + row];
68   - }
69   -
70   - /// Return the number of bytes in the matrix array
71   - size_t bytes() { return m_rows * m_cols * sizeof(type); }
72   -
73   - /// Create a deep copy of the current matrix and return it
74   - Matrix copy() {
75   - Matrix result(m_rows, m_cols);
76   - memcpy(result.m_ptr, m_ptr, bytes());
77   - return result;
78   - }
79   -
80   - /// Calculate the determinant of the matrix
81   - type det() {
82   - force_square(); //throw an exception if the matrix isn't square
83   - Matrix tmp = copy(); //copy the current matrix to create space for LUP decomposition
84   - Matrix<int> idx = getrf(tmp); //perform LUP decomposition
85   -
86   - // multiply all elements along the diagonal
87   - type determinant = get(0, 0); //initialize the determinant
88   - for (int i = 1; i < m_rows; i++) //for each remaining element of the diagonal
89   - determinant *= tmp.get(i, i); //calculate the running factor
90   - tmp.del(); //delete memory for the temporary matrix
91   -
92   - // calculate the number of row swaps
93   - int swaps = 0; //initialize the number of row swaps
94   - for (int i = 0; i < m_rows; i++) //iterate through each element of the swap vector
95   - if (idx.get(i, 0) != i + 1) //if the element does not equal the row number
96   - swaps++; //a row was swapped, so increment the swap counter
97   - if (swaps % 2 != 0) //if there was an odd number of row swaps
98   - determinant *= (-1); //multiply the determinant by -1
99   - return determinant; //return the determinant
100   - }
101   -
102   - /// Calculate the trace of the matrix
103   - type trace() {
104   - force_square(); //throw an exception if the matrix isn't square
105   - type trace = 0; //initialize the trace summation
106   - for (int i = 0; i < m_rows; i++) //for each element along the diagonal
107   - trace += get(i, i); //sum the trace
108   - return trace; //return the matrix trace result
109   - }
110   -
111   - /// Output the matrix to a CSV file
112   - void csv(std::string filename) {
113   - std::ofstream outfile(filename.c_str()); //open the file for writing
114   - for (int r = 0; r < m_rows; r++) { //for each row
115   - for (int c = 0; c < m_cols; c++) { //for each column
116   - outfile << get(r, c); //output the value to a file
117   - if (c != m_cols - 1) outfile << ","; //if this isn't the last column, add a comma
118   - }
119   - if (r != m_rows - 1) outfile << std::endl; //if this isn't the last row, end the line
120   - }
121   - }
122   -
123   - /* Static functions used to generate specialized matrices */
124   -
125   - /// Static function for generating an identity matrix
126   - static Matrix I(int dim) {
127   - Matrix result(dim, dim); //create and allocate space for the matrix
128   - memset(result.m_ptr, 0, dim * dim * sizeof(type)); //set all elements to zero
129   - for (int i = 0; i < dim; i++) //for each element on the diagonal
130   - result.set(i, i, (type)1); //set that element to 1
131   - return result; //return the completed matrix
132   - }
133   -
134   - /// Static function for generating a geometric matrix for robustness testing
135   - static Matrix geometric_matrix(int dim) {
136   - Matrix result(dim, dim); //create and allocate space for the matrix
137   - for (int r = 0; r < dim; r++) //for each row
138   - for (int c = 0; c < dim; c++) //for each column
139   - result.set(r, c, pow((type)(2 + r), c)); //set the geometric coefficient
140   - return result; //return the completed matrix
141   - }
142   -
143   - /// Static function analytically generates a solution (right hand side) vector for a geometric matrix
144   - static Matrix geometric_vector(int dim) {
145   - Matrix result(dim, 1); //generate and allocate space for the vector
146   - for (int r = 0; r < dim; r++)
147   - result.set(r, 0, (pow((type)(2 + r), dim) - 1) / (type)(r + 1));
148   - return result; //return the completed vector
149   - }
150   -
151   - /* Static functions used for LAPACK and BLAS */
152   -
153   - /// Calls the appropriate LAPACK GETRF function for performing LUP decomposition
154   - static Matrix<int> getrf(Matrix<float> M) {
155   - Matrix<int> idx(M.m_rows, 1);
156   - LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
157   - return idx;
158   - }
159   - static Matrix<int> getrf(Matrix<double> M) {
160   - Matrix<int> idx(M.m_rows, 1);
161   - LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
162   - return idx;
163   - }
164   - };
165   -}
166 0 \ No newline at end of file