3

C++ で LAPACK を使用して、複雑な行列を反転しています。具体的には、私が使用している2つの機能は次のとおりです。

zgetrfLU分解用。

zgetri反転のために。

コードを最適化するという私の目標として、タイミングについて質問があります。LAPACK で一般的な行列反転法を使用する場合 (使用する関数がより適切で迅速な場合はお知らせください)、関数のタイミングは行列の値とは無関係ですか?

たとえば、人口密度の高い行列を反転するよりも単位行列を反転する方が速いでしょうか?

繰り返しますが、複素行列の一般的な LAPACK 反転に関してこの質問をしていることを強調したいと思います。私は、使用できるさまざまな三重対角関数とバンド関数について知っています。

行列のすべての要素が複雑な double であると想定しています。

ありがとう、ケビン

4

1 に答える 1

2

Kieran Cooney が仮定したように、LAPACK は単位行列をランダムな密行列よりも桁違いに高速に反転します。以下のテストを使用すると、次の結果が得られます (サンプル サイズ = 1 ですが、要点を証明しています)。

サイズ変更
情報: 0
合計時間 (ランダム) = 2389 ミリ秒。
情報: 0
合計時間 (ID) = 14 ミリ秒。 

#include "lapacke.h"

#include <iostream>
#include <vector>
#include <Eigen/Core>
#include <chrono>

lapack_int getSize(lapack_int n, lapack_complex_double* a,
    const lapack_int* ipiv, lapack_complex_double* work)
{
    lapack_complex_double resizetome;
    lapack_int hello = -1;
    lapack_int info = -1;

    LAPACK_zgetri(&n, a, &n, ipiv, &resizetome, &hello, &info);

    return lapack_int(resizetome.real());

}
void invert(lapack_int n, lapack_complex_double* a,
    lapack_int* ipiv, lapack_complex_double* work, lapack_int lwork, lapack_int *info)
{
    // LU factor
    LAPACK_zgetrf(&n, &n, a, &n, ipiv, info);

    // Invert
    LAPACK_zgetri(&n, a, &n, ipiv, work, &lwork, info);
}

int main(int argc, char* argv[]) {

    int sz = 1000;

    int ln = sz;
    int llda = sz;
    int lipiv = 1;
    int llwork = -1;
    int linfo = 0;

    srand(time(NULL));

    typedef Eigen::MatrixXcd lapackMat;
    lapackMat ident = lapackMat::Identity(sz, sz).eval();
    lapackMat randm = lapackMat::Random(sz, sz);
    lapackMat work = lapackMat::Zero(1, 1);
    Eigen::VectorXi ipvt(sz);
    randm;

    work.resize(1,
        getSize(ln, randm.data(), ipvt.data(), work.data())
        );

    std::cout << "Resized\n";

    // Timing for random matrix
    {
        auto startTime = std::chrono::high_resolution_clock::now();

        invert(ln, randm.data(), ipvt.data(), work.data(), llwork, &linfo);

        auto endTime = std::chrono::high_resolution_clock::now();

        std::cout << "Info: " << linfo << "\n";

        std::cout << "Total Time (random) = " <<
            std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count()
            << " milliseconds.\n";
    }

    // Timing for identity matrix
    {
        auto startTime = std::chrono::high_resolution_clock::now();

        invert(ln, ident.data(), ipvt.data(), work.data(), llwork, &linfo);

        auto endTime = std::chrono::high_resolution_clock::now();

        std::cout << "Info: " << linfo << "\n";

        std::cout << "Total Time (identity) = " <<
            std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count()
            << " milliseconds.\n";

    }

    return 0;
}
于 2015-08-10T13:52:48.887 に答える