gpt4 book ai didi

c++ - 复杂对称三对角矩阵的快速矩阵指数

转载 作者:塔克拉玛干 更新时间:2023-11-03 02:24:09 24 4
gpt4 key购买 nike

基本上我需要以上这些。我已经搜索了谷歌,但找不到实现它的方法。

我在这里找到了这个函数 http://www.guwi17.de/ublas/examples/但它太慢了。我什至按照 MATLAB 的例程编写了自己的 Pade Approximation,但它只比链接中的快一点点。

让我吃惊的是 Mathematica 计算矩阵指数的速度有多快(我不知道它是否关心矩阵是否为三边形)。

有人能帮忙吗?

编辑:这是我想出的,有什么意见吗?希望对 future 的读者有用

我已经离开 C++ 一段时间了,所以下面的代码可能有点乱/慢,所以如果你看到改进请赐教。

//Program will compute the matrix exponential expm( i gA ). where i = sqrt(-1) and gA is a matrix

#include <iostream>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>


int main() {

int n = 100;

unsigned i = 0;
unsigned j = 0;

gsl_complex img = gsl_complex_rect (0.,1.);

gsl_matrix * gA = gsl_matrix_alloc (n, n);


//Set gA:
for ( i = 0; i<n; i++) {
for ( j = 0; j<=i; j++) {
if( i == j ) {
if( i == 0 || i == n-1 ) {
gsl_matrix_set (gA, i, j, 1.);
} else {
gsl_matrix_set (gA, i, j, 2.);
}
} else if( j == i-1 ) {
gsl_matrix_set (gA, i, j, 1.);
gsl_matrix_set (gA, j, i, 1.);
} else {
gsl_matrix_set (gA, i, j, 0.);
gsl_matrix_set (gA, j, i, 0.);
}
}
}


//get SVD of gA
gsl_matrix * gA_t = gsl_matrix_alloc (n, n);
gsl_matrix_transpose_memcpy (gA_t, gA);

gsl_matrix * U = gsl_matrix_alloc (n, n);
gsl_matrix * V= gsl_matrix_alloc (n, n);
gsl_vector * S = gsl_vector_alloc (n);


gsl_vector * work = gsl_vector_alloc (n);
gsl_linalg_SV_decomp (gA_t, V, S, work);
gsl_vector_free(work);

gsl_matrix_memcpy (U, gA_t);


//Take exponential of S and form matrix
gsl_matrix_complex * Sp = gsl_matrix_complex_alloc (n, n);
gsl_matrix_complex_set_zero (Sp);
for (i = 0; i < n; i++) {
gsl_matrix_complex_set (Sp, i, i, gsl_complex_exp ( gsl_complex_mul_real ( img, gsl_vector_get(S, i) ) ) ); // Vector 'S' to matrix 'Sp'
}

gsl_matrix_complex * Uc = gsl_matrix_complex_alloc (n, n);


//convert U to a complex matrix for next step
for (i = 0; i < n; i++) {
for ( j = 0; j<n; j++) {
gsl_matrix_complex_set (Uc, i, j, gsl_complex_rect ( gsl_matrix_get(U, i, j), 0. ) );
}
}


//recombine U S Utranspose
gsl_matrix_complex * tA = gsl_matrix_complex_alloc (n, n);
gsl_matrix_complex * eA = gsl_matrix_complex_alloc (n, n);
gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, Uc, Sp, GSL_COMPLEX_ZERO, tA);
gsl_blas_zgemm (CblasNoTrans, CblasTrans, GSL_COMPLEX_ONE, tA, Uc, GSL_COMPLEX_ZERO, eA);


//Print result
std::cout << "eA:" << std::endl;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
printf ("%g + %g i\n", GSL_REAL(gsl_matrix_complex_get (eA, i, j)), GSL_IMAG(gsl_matrix_complex_get (eA, i, j)));
std::cout << "\n" << std::endl;


//Free up
gsl_matrix_free(gA_t);
gsl_matrix_free(U);
gsl_matrix_free(gA);
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_matrix_complex_free(Uc);
gsl_matrix_complex_free(tA);
gsl_matrix_complex_free(eA);

return 0;
}

最佳答案

我在问题中发布的上面的代码运行良好。我发现的另一个改进是通过特征值缩放 V 拷贝中的列,而不是执行完整的矩阵乘法。由于 zgemm 是此代码中最慢的部分,删除其中一个 zgemm 函数可将其速度提高近两倍。我将很快在此处发布完整的代码 list 。

关于c++ - 复杂对称三对角矩阵的快速矩阵指数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10124754/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com