Diff: Matrix.cpp

Revision:

2:4a6b46653abf

Parent:

1:2eca9b376580

--- a/Matrix.cpp	Thu Apr 16 08:51:04 2015 +0000
+++ b/Matrix.cpp	Mon Apr 20 14:54:55 2015 +0000
@@ -32,6 +32,16 @@
     memcpy(components, m.GetpComponents(), sizeof(float)*row*col);
 }
 
+Matrix Matrix::operator-() const{
+    Matrix retMat(*this);
+
+    for (int i = 0; i < row * col; i++) {
+        retMat.components[i] = - this->components[i];
+    }
+
+    return retMat;
+}
+
 Matrix& Matrix::operator=(const Matrix& m) {
     if (this == &m) return *this;
     row = m.row;
@@ -139,6 +149,61 @@
     return retVal*decSign;
 }
 
+float Matrix::det() const {
+    if (row != col) AbortWithMsg("failed to calculate det : matrix is not square");
+    
+    Matrix temp(*this);
+    int decSign = 1;
+
+    for (int j = 0; j < col - 1; j++) {
+
+        // 列内のみで最大の要素を探す
+        int maxNo = j;
+        for (int k = j; k < row; k++) {
+            if (temp.components[maxNo * col + j] < temp.components[k * col + j]) maxNo = k;
+        }
+        if (maxNo != j) {
+            temp.SwapRow(j + 1, maxNo + 1);
+            decSign *= -1;
+        }
+        // 列内の最大要素が小さ過ぎる場合、行内の最大要素も探す
+        if (fabs(temp.components[j * col + j]) < NEARLY_ZERO) {
+            maxNo = j;
+            for (int k = j; k < col; k++) {
+                if (temp.components[j * col + maxNo] < temp.components[j * col + k])maxNo = k;
+            }
+            if (maxNo != j) {
+                temp.SwapCol(j + 1, maxNo + 1);
+                decSign *= -1;
+            }
+
+            // 列内、行内の最大要素を選んでも小さすぎる場合はエラー
+            if (fabs(temp.components[j * col + j]) < NEARLY_ZERO) {
+                if (row != col) AbortWithMsg("failed to calculate det : Division by Zero");
+            }
+        }
+
+        float c1 = 1.0f / temp.components[j * col + j];
+
+        for (int i = j + 1; i < row; i++) {
+            float c2 = temp.components[i * col + j] * c1;
+            for (int k = j; k < col; k++) {
+                temp.components[i * col + k] = temp.components[i * col + k] - c2 * temp.components[j * col + k];
+            }
+        }
+        
+    }
+
+    if (fabs(temp.components[(row - 1) * col + (col - 1)]) < NEARLY_ZERO) return 0.0f;
+
+    float retVal = 1.0f;
+    for (int i = 0; i < row; i++) {
+        retVal *= temp.components[i * col + i];
+    }
+
+    return retVal * decSign;
+}
+
 Matrix Matrix::LU_Decompose(int* sign, Matrix* p) const{
     if (row != col) AbortWithMsg("failed to LU decomposition: matrix is not square");
     if (sign != 0) *sign = 1;
@@ -295,4 +360,15 @@
     memcpy(components + (rowNo2 - 1) * col, temp, sizeof(float) * col);
 
     delete[] temp;
+}
+
+void Matrix::SwapCol(int colNo1, int colNo2) {
+    if (colNo1 > col || colNo2 > col) AbortWithMsg("Index Out of Bounds Error !!");
+    float temp = 0.0f;
+
+    for (int i = 0; i < row; i++) {
+        temp = components[i * col + colNo1];
+        components[i * col + colNo1] = components[i * col + colNo2];
+        components[i * col + colNo2] = temp;
+    }
 }
\ No newline at end of file