Инвертирование матриц с использованием LAPACK и синхронизации
Я использую LAPACK на C ++, чтобы инвертировать сложную матрицу. В частности, две функции, которые я использую:
zgetrf для разложения LU.
zgetri для инверсии.
Теперь, когда моя цель оптимизировать мой код, у меня есть вопрос о сроках. Используя общий метод обращения к матрице с LAPACK (если у вас есть лучшие / более быстрые функции для использования, пожалуйста, дайте мне знать), является ли синхронизация функции (ей) независимой от значений в матрице?
Например, будет ли быстрее инвертировать единичную матрицу, чем инвертировать густонаселенную матрицу?
Опять же, я хотел бы подчеркнуть, что я задаю этот вопрос в отношении общей инверсии LAPACK сложных матриц. Я знаю о различных трехдиагональных и полосовых функциях, которые можно использовать.
Я предполагаю, что все элементы в матрице являются сложными двойными числами.
Спасибо,
Kevin
Решение
Как предположил Киран Куни, LAPACK инвертирует порядки порядка единичной матрицы быстрее, чем случайная плотная матрица. Использование теста ниже дает мне следующие результаты (размер выборки = 1, но доказывает суть):
Измененный размер
Информация: 0
Общее время (случайное) = 2389 миллисекунд.
Информация: 0
Общее время (идентичность) = 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;
}
Другие решения