1. 程式人生 > >使用Eigen解Ax=b線性方程組

使用Eigen解Ax=b線性方程組

#include <iostream>
using namespace std;
#include <ctime>
// Eigen 部分
#include <Eigen/Core>
// 稠密矩陣的代數運算(逆,特徵值等)
#include <Eigen/Dense>

#define MATRIX_SIZE 100

/****************************
* 本程式演示了 Eigen 基本型別的使用
****************************/

int main( int argc, char** argv )
{
    // 解方程
    // 我們求解 A * x = b 這個方程
    // 直接求逆自然是最直接的,但是求逆運算量大

    Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > A1;
    A1 = Eigen::MatrixXd::Random( MATRIX_SIZE, MATRIX_SIZE );

    Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > b1;
    b1 = Eigen::MatrixXd::Random( MATRIX_SIZE, 1 );

    clock_t time_stt = clock(); // 計時
    // 直接求逆
    Eigen::Matrix<double,MATRIX_SIZE,1> x = A1.inverse()*b1;
    cout <<"time use in normal inverse is " << 1000* (clock() - time_stt)/(double)CLOCKS_PER_SEC << "ms"<< endl;
    cout<<x<<endl;
    // QR分解colPivHouseholderQr()
    time_stt = clock();
    x = A1.colPivHouseholderQr().solve(b1);
    cout <<"time use in Qr decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;
    //QR分解fullPivHouseholderQr()
    time_stt = clock();
    x = A1.fullPivHouseholderQr().solve(b1);
    cout <<"time use in Qr decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;
    /* //llt分解 要求矩陣A正定
    time_stt = clock();
    x = A1.llt().solve(b1);
    cout <<"time use in llt decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;*/
    /*//ldlt分解  要求矩陣A正或負半定
    time_stt = clock();
    x = A1.ldlt().solve(b1);
    cout <<"time use in ldlt decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;*/
    //lu分解 partialPivLu()
    time_stt = clock();
    x = A1.partialPivLu().solve(b1);
    cout <<"time use in lu decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;
    //lu分解(fullPivLu()
    time_stt = clock();
    x = A1.fullPivLu().solve(b1);
    cout <<"time use in lu decomposition is " <<1000*(clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;
    cout <<x<<endl;

    return 0;

}