Matrix.h
5.59 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
#include <iostream>
#include <fstream>
#include <complex>
#ifndef LAPACK_COMPLEX_CUSTOM
#define LAPACK_COMPLEX_CUSTOM
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#include "lapacke.h"
#endif
namespace tira {
template <typename type>
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<int> 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<int> getrf(Matrix<float> M) {
Matrix<int> 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<int> getrf(Matrix<double> M) {
Matrix<int> 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;
}
};
}