Show toolbar

2013年7月18日 星期四

Solving Matrix Inverse using Pointer

標題:使用指標求解n維矩陣
VC++ (main.cpp):
#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <string>
#include <regex>
#include <iomanip>
using namespace std;

void print(double *Arr, int m)
{
    int i = 0, j = 0;
    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            cout << *(Arr+i*m+j) << setprecision(2) << "\t";
        }
        cout << endl;
    }
}
void printInverse(double *Matrix, double *Inverse, int m)
{
    int i = 0, j = 0;
    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            cout << *(Matrix+i*m+j) << setprecision(2) << "\t";
        }
        cout << "|\t";
        for(j=0;j<m;j++)
        {
            cout << *(Inverse+i*m+j) << setprecision(2) << "\t";
        }
        cout << endl;
    }
}

void readfile(string filePath, double *Matrix)
{
    int m = 0, n = 0, i = 0;
    string lines = "";
    string::const_iterator start;
    match_results<string::const_iterator> matchMatrix; //迭代器
    regex regMatrix("([0-9]+)+");
    ifstream outputFile(filePath);
    if (outputFile.is_open()) //Check File
    {
        while(outputFile.good())
        {
            getline(outputFile, lines);
            start = lines.begin();
            while (regex_search(start, lines.cend(), matchMatrix, regMatrix))
            {
                *(Matrix+i) = atof(matchMatrix[1].str().c_str());
                i++;
                start = matchMatrix[0].second;
            }
            m++;
        }
        outputFile.close();
    }
}
double plusnum(double num)
{
    if(num < 0) return -1 * num;
    return num;
}

void rowSwap(double *Matrix, int ma, int mb, int m)
{
    int i = 0;
    double *Temp = new double[m];
    for(i=0;i<m;i++)
    {
        *(Temp+i) = *(Matrix+ma*m+i);
        *(Matrix+ma*m+i) = *(Matrix+mb*m+i);
        *(Matrix+mb*m+i) = *(Temp+i);
    }
}

void inverse(double *Matrix, double *Inverse, int m)
{
    int i = 0, j = 0, k = 0;
    double temp = 0;
    //*(Matrix + rows * m + columns) 看成 Matrix[rows][columns]
    for(int i=0; i<m; i++)
    {
        for(int j=i+1; j<m; j++)
        {
            if(plusnum(*(Matrix + j * m + i)) > plusnum(*(Matrix + i * m + i)))
            {
                rowSwap(Matrix, i, j, m); // 互換Matrix矩陣兩行
                rowSwap(Inverse, i, j, m); // 互換Inverse矩陣兩行
                //cout << "[" << i << "][" << j << "]Swap:" << endl;
                //printInverse(Matrix, Inverse, m);
            }
        }
        // 消去
        for(int j=0; j<m; j++) //第 i 次消去運算,以第 i 行作為主要行,將其上下各行的第 i 列元素化為零
        {
            if(j != i)
            {
                // 進行消去運算,將第 i 列上下方的元素化為零
                double temp = *(Matrix + j * m + i) / *(Matrix + i * m + i);
                
                // 注意此處 k 可從 i+1 開始
                //for(int k=i+1; k<m; k++)
                for(int k=0; k<m; k++)
                {
                    *(Matrix + j * m + k) -= (*(Matrix + i * m + k) * temp);
                }
                // 注意此處 k 要從 0 開始
                for(int k=0; k<m; k++)
                {
                    *(Inverse + j * m + k) -= (*(Inverse + i * m + k) * temp);
                }
                //cout << "[" << i << "][" << j << "]Matrix:" << endl;
                //printInverse(Matrix, Inverse, m);
            }
        }
    }

    // Matrix 已經是對角矩陣,只要在將其對角線元素化為 1 即可
    // 此時僅需對 Inverse 矩陣作運算
    //cout << "[" << i << "][" << j << "]Inverse:" << endl;
    //printInverse(Matrix, Inverse, m);
    for(int i=0; i<m; i++)
    {
        for(int j=0; j<m; j++)
        {
            *(Inverse + i * m + j) /= *(Matrix + i * m + i);
        }
        //cout << "[" << i << "][" << j << "]Inverse:" << endl;
        //printInverse(Matrix, Inverse, m);
    }
}
 
int main(int argc, char *argv[])
{
    const int m = 4;
    int i = 0, j = 0;
    double Matrix[m][m] = {};
    double Inverse[m][m] = {};

    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            if(i == j)
                Inverse[i][j] = 1;
            else
                Inverse[i][j] = 0;
        }
    }
    
    readfile("matrix.txt", &Matrix[0][0]);

    cout << "Matrix:" << endl;
    print(&Matrix[0][0], m);

    inverse(&Matrix[0][0], &Inverse[0][0], m);
    cout << endl;
 
    cout << "Inverse:" << endl;
    print(&Inverse[0][0], m);

    system("pause");
    return 0;
}
TEXT (matrix.txt):

說明:
先讀取矩陣數據再利用指標的方式求得n維度反矩陣的解。

沒有留言:

張貼留言