Fabio Bobrow
Cubli library
Diff: Matrix/Matrix.cpp
- Revision:
- 1:085840a3d767
- Parent:
- 0:431ee55036ca
- Child:
- 2:dc7840a67f77
--- a/Matrix/Matrix.cpp Wed Feb 13 13:07:37 2019 +0000
+++ b/Matrix/Matrix.cpp Wed Feb 13 17:06:28 2019 +0000
@@ -2,50 +2,74 @@
Matrix::Matrix() : _rows(1), _cols(1)
{
- // Allocate memory
- _data = (float **)malloc(sizeof(float *));
- _data[0] = (float *)malloc(sizeof(float));
- // Initialize values
- _data[0][0] = 0.0;
+ // Allocate memory
+ _data = (float **)malloc(sizeof(float *));
+ _data[0] = (float *)malloc(sizeof(float));
+ // Initialize values
+ _data[0][0] = 0.0;
}
Matrix::Matrix(int rows, int cols) : _rows(rows), _cols(cols)
{
- // Allocate memory
- _data = (float **)malloc(_rows * sizeof(float *)); // = new float*[_rows];
- for (int i = 1; i <= _rows; i++)
- {
- _data[i-1] = (float *)malloc(_cols * sizeof(float)); // = new float[_cols];
- }
- // Initialize values
- for (int i = 1; i <= _rows; i++)
- {
- for (int j = 1; j <= _cols; j++)
- {
- _data[i-1][j-1] = 0.0;
+ // Allocate memory
+ _data = (float **)malloc(_rows * sizeof(float *)); // = new float*[_rows];
+ for (int i = 1; i <= _rows; i++) {
+ _data[i-1] = (float *)malloc(_cols * sizeof(float)); // = new float[_cols];
}
- }
+ // Initialize values
+ for (int i = 1; i <= _rows; i++) {
+ for (int j = 1; j <= _cols; j++) {
+ _data[i-1][j-1] = 0.0;
+ }
+ }
}
Matrix::~Matrix()
{
// Free memory
- for (int i = 0; i < _rows; ++i) {
- free(_data[i]);
+ for (int i = 1; i <= _rows; i++) {
+ free(_data[i-1]);
}
free(_data);
}
+Matrix& Matrix::operator=(const Matrix& m)
+{
+ for (int i = 1; i <= _rows; i++) {
+ for (int j = 1; j <= _cols; j++) {
+ _data[i-1][j-1] = m._data[i-1][j-1];
+ }
+ }
+ return *this;
+}
+
float& Matrix::operator()(int row, int col)
{
- // Return cell data
- return _data[row-1][col-1];
+ // Return cell data
+ return _data[row-1][col-1];
}
Matrix operator+(const Matrix& A, const Matrix& B)
{
- // Auxiliary matrix where C = A+B
+ /*// Auxiliary matrix where C = A+B
+ Matrix C(A._rows,A._cols);
+ for (int i = 1; i <= A._rows; i++)
+ {
+ for (int j = 1; j <= A._cols; j++)
+ {
+ C._data[i-1][j-1] += A._data[i-1][j-1]+B._data[i-1][j-1];
+ }
+ }
+ return C;
+ */
+
+ /*Matrix temp(m1);
+ return (temp += m2);
+ */
+
+ //Matrix A(m1);
+ //Matrix B(m2);
Matrix C(A._rows,A._cols);
for (int i = 1; i <= A._rows; i++)
{
@@ -61,12 +85,10 @@
{
// Auxiliary matrix where C = A-B
Matrix C(A._rows,A._cols);
- for (int i = 1; i <= A._rows; i++)
- {
- for (int j = 1; j <= A._cols; j++)
- {
- C._data[i-1][j-1] += A._data[i-1][j-1]-B._data[i-1][j-1];
- }
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._cols; j++) {
+ C._data[i-1][j-1] += A._data[i-1][j-1]-B._data[i-1][j-1];
+ }
}
return C;
}
@@ -75,15 +97,12 @@
{
// Auxiliary matrix where C = A*B
Matrix C(A._rows,B._cols);
- for (int i = 1; i <= A._rows; i++)
- {
- for (int j = 1; j <= B._cols; j++)
- {
- for (int k = 1; k <= A._cols; k++)
- {
- C._data[i-1][j-1] += A._data[i-1][k-1]*B._data[k-1][j-1];
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= B._cols; j++) {
+ for (int k = 1; k <= A._cols; k++) {
+ C._data[i-1][j-1] += A._data[i-1][k-1]*B._data[k-1][j-1];
+ }
}
- }
}
return C;
}
@@ -92,12 +111,34 @@
{
// Auxiliary matrix where C = scalar*A
Matrix C(A._rows,A._cols);
- for (int i = 1; i <= A._rows; i++)
- {
- for (int j = 1; j <= A._cols; j++)
- {
- C._data[i-1][j-1] += scalar*A._data[i-1][j-1];
- }
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._cols; j++) {
+ C._data[i-1][j-1] += scalar*A._data[i-1][j-1];
+ }
+ }
+ return C;
+}
+
+Matrix operator*(const Matrix& A, float scalar)
+{
+ // Auxiliary matrix where C = scalar*A
+ Matrix C(A._rows,A._cols);
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._cols; j++) {
+ C._data[i-1][j-1] += A._data[i-1][j-1]*scalar;
+ }
+ }
+ return C;
+}
+
+Matrix operator/(const Matrix& A, float scalar)
+{
+ // Auxiliary matrix where C = scalar*A
+ Matrix C(A._rows,A._cols);
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._cols; j++) {
+ C._data[i-1][j-1] += A._data[i-1][j-1]/scalar;
+ }
}
return C;
}
@@ -106,10 +147,8 @@
{
// Auxiliary matrix where C = A'
Matrix C(A._cols, A._rows);
- for (int i = 1; i <= C._rows; i++)
- {
- for (int j = 1; j <= C._cols; j++)
- {
+ for (int i = 1; i <= C._rows; i++) {
+ for (int j = 1; j <= C._cols; j++) {
C._data[i-1][j-1] = A._data[j-1][i-1];
}
}
@@ -118,33 +157,25 @@
Matrix inverse(const Matrix& A)
{
- // Apply A = LDL' factorization where L is a lower triangular matrix and D
+ // Apply A = LDL' factorization where L is a lower triangular matrix and D
// is a block diagonal matrix
Matrix L(A._cols, A._rows);
Matrix D(A._cols, A._rows);
float L_sum;
float D_sum;
- for (int i = 1; i <= A._rows; i++)
- {
- for (int j = 1; j <= A._rows; j++)
- {
- if (i >= j)
- {
- if (i == j)
- {
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._rows; j++) {
+ if (i >= j) {
+ if (i == j) {
L._data[i-1][j-1] = 1.0;
D_sum = 0.0;
- for (int k = 1; k <= (i-1); k++)
- {
+ for (int k = 1; k <= (i-1); k++) {
D_sum += L._data[i-1][k-1]*L._data[i-1][k-1]*D._data[k-1][k-1];
}
D._data[i-1][j-1] = A._data[i-1][j-1] - D_sum;
- }
- else
- {
+ } else {
L_sum = 0.0;
- for (int k = 1; k <= (j-1); k++)
- {
+ for (int k = 1; k <= (j-1); k++) {
L_sum += L._data[i-1][k-1]*L._data[j-1][k-1]*D._data[k-1][k-1];
}
L._data[i-1][j-1] = (1.0/D._data[j-1][j-1])*(A._data[i-1][j-1]-L_sum);
@@ -152,26 +183,19 @@
}
}
}
- // Compute the inverse of L and D matrices
+ // Compute the inverse of L and D matrices
Matrix L_inv(A._cols, A._rows);
Matrix D_inv(A._cols, A._rows);
float L_inv_sum;
- for (int i = 1; i <= A._rows; i++)
- {
- for (int j = 1; j <= A._rows; j++)
- {
- if (i >= j)
- {
- if (i == j)
- {
+ for (int i = 1; i <= A._rows; i++) {
+ for (int j = 1; j <= A._rows; j++) {
+ if (i >= j) {
+ if (i == j) {
L_inv._data[i-1][j-1] = 1.0/L._data[i-1][j-1];
D_inv._data[i-1][j-1] = 1.0/D._data[i-1][j-1];
- }
- else
- {
+ } else {
L_inv_sum = 0.0;
- for (int k = j; k <= (i-1); k++)
- {
+ for (int k = j; k <= (i-1); k++) {
L_inv_sum += L._data[i-1][k-1]*L_inv._data[k-1][j-1];
}
L_inv._data[i-1][j-1] = -L_inv_sum;
@@ -179,7 +203,18 @@
}
}
}
- // Compute the inverse of A matrix in terms of the inverse of L and D
- // matrices
+ // Compute the inverse of A matrix in terms of the inverse of L and D
+ // matrices
return transpose(L_inv)*D_inv*L_inv;
}
+
+float norm(const Matrix& A)
+{
+ // Auxiliary matrix where C = A'
+ float n;
+ n = 0.0;
+ for (int i = 1; i <= A._rows; i++) {
+ n += A._data[i-1][0]*A._data[i-1][0];
+ }
+ return sqrt(n);
+}