如何解决对 gsl_blas_dgemm() 的调用未按预期工作
我正在研究一个需要求解 3 个变量的 3 个联立方程的问题。阅读:
https://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html
我发现这可以通过确定 3 个方程的矩阵表示的逆矩阵来解决。我发现这个 C 代码使用 github 上的 GSL 库来计算逆矩阵:
https://gist.github.com/bjd2385/7f4685e703f7437e513608f41c65bbd7
(非常感谢它的作者 Doyle 先生。)
我读过,如果一个矩阵乘以其逆矩阵应该得到一个单位矩阵(一个矩阵,对角线上有 1.0 秒,其他地方都是 0.0 秒)。所以我认为作为对上述 github 代码的健全性检查,我可以修改它以将结果逆乘以原始矩阵并显示结果。如果结果是单位矩阵,则证明逆计算。
我发现,至少对于简单的 2X2 矩阵情况,虽然逆计算的结果看起来是正确的,但随后的矩阵乘法不会产生单位矩阵。
我是这个 GSL 库的新手,所以也许我只是没有正确调用 gsl_blas_dgemm() 库函数来执行矩阵乘法。
/* A simple example of inverting a matrix using the gsl */
/* 1-26-2021,modified to sanity check result */
/* code doesn't seem to work */
#define HAVE_INLINE
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_linalg.h>
gsl_matrix *invert_a_matrix(gsl_matrix *matrix);
void print_mat_contents(gsl_matrix *matrix);
void randomize_mat_contents(gsl_matrix *matrix);
static size_t size = 2;
/************************************************************
* PROCEDURE: invert_a_matrix
*
* DESCRIPTION: Invert a matrix using GSL.
*
* RETURNS:
* gsl_matrix pointer
*/
gsl_matrix *
invert_a_matrix(gsl_matrix *matrix)
{
gsl_permutation *p = gsl_permutation_alloc(size);
int s;
// Compute the LU decomposition of this matrix
gsl_linalg_LU_decomp(matrix,p,&s);
// Compute the inverse of the LU decomposition
gsl_matrix *inv = gsl_matrix_alloc(size,size);
gsl_linalg_LU_invert(matrix,inv);
gsl_permutation_free(p);
return inv;
}
/************************************************************
* PROCEDURE: print_mat_contents
*
* DESCRIPTION: Print the contents of a gsl-allocated matrix
*
* RETURNS:
* None.
*/
void
print_mat_contents(gsl_matrix *matrix)
{
size_t i,j;
double element;
for (i = 0; i < size; ++i) {
for (j = 0; j < size; ++j) {
element = gsl_matrix_get(matrix,i,j);
printf("%f ",element);
}
printf("\n");
}
}
/************************************************************
* PROCEDURE: randomize_mat_contents
*
* DESCRIPTION: Overwrite entries in matrix with randomly
* generated values.
*
* RETURNS:
* None.
*/
void
randomize_mat_contents(gsl_matrix *matrix)
{
size_t i,j;
double random_value;
double range = 1.0 * RAND_MAX;
for (i = 0; i < size; ++i) {
for (j = 0; j < size; ++j) {
// generate a random value
random_value = rand() / range;
// set entry at i,j to random_value
gsl_matrix_set(matrix,j,random_value);
}
}
}
int
main(void)
{
srand(time(NULL));
gsl_matrix *mat = gsl_matrix_alloc(size,size);
// fill this matrix with random doubles
randomize_mat_contents(mat);
// let's see the original Now
printf("Original matrix:\n");
print_mat_contents(mat);
printf("\n");
// compute the matrix inverse
gsl_matrix *inverse = invert_a_matrix(mat);
printf("Inverted matrix:\n");
print_mat_contents(inverse);
printf("\n");
gsl_matrix *product = gsl_matrix_calloc(size,size);
// if inverse is truly the inverse of mat,then mat * inverse should = identity matrix
printf("product before:\n");
print_mat_contents(product);
printf("\n");
int error;
// neither of these results in an identity matrix,8^(
error = gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,inverse,mat,0.0,product);
// error = gsl_blas_dgemm(CblasNoTrans,product);
if (error) {
fprintf(stderr,"gsl_blas_dgemm returned %d\n",error);
}
printf("inverse * mat:\n");
print_mat_contents(product);
gsl_matrix_free(mat);
gsl_matrix_free(inverse);
gsl_matrix_free(product);
return 0;
}
总结:
创建一个包含随机内容的矩阵,打印出来。 计算它的逆,打印逆。 调用 gsl_blas_dgemm() 将矩阵乘以其逆矩阵,打印应该是单位矩阵的内容。
在我的 ubuntu 20.04 笔记本电脑上编译、链接和运行程序时,我得到了这个:
~/opengl/matrix_code/2by2_inverter$ gcc inverter.c -lgsl -lgslcblas -lm
~/opengl/matrix_code/2by2_inverter$ ./a.out
Original matrix:
0.317588 0.113800
0.280836 0.190114
Inverted matrix:
6.689703 -4.004360
-9.882010 11.175224
product before:
0.000000 0.000000
0.000000 0.000000
inverse * mat:
-1.416400 0.402961
6.743603 -0.124569
~/opengl/matrix_code/2by2_inverter$
现在,如果我“手动”将原始矩阵乘以它的逆矩阵,我发现结果是一个 2X2 单位矩阵:
Upper left corner: 0.317588*6.689703+0.113800*-9.882010 = .999996658 ~= 1.0
Upper right corner: 0.317588*-4.004360+0.113800*11.175224 = .000003808 ~= 0.0
Lower left corner: 0.280836*6.689703+0.190114*-9.882010 = .000000982 ~= 0.0
Lower right corner: 0.280836*-4.004360+0.190114*11.175224 = .999998091 ~= 1.0
当然,单位矩阵的坐标不完全是 1.0 和 0.0,但预计会有一些误差。由此我得出结论,函数 invert_a_matrix() 正在做正确的事情,至少对于 2X2 矩阵。
但尽我所能,我无法弄清楚如何调用 gsl_blas_dgemm() 来生成单位矩阵。
请注意,我通过以下方式从 Ubuntu 存储库安装了 GSL 库:
~$ sudo apt-get install libgsl-dev
任何关于我做错了什么的线索?
提前致谢
解决方法
我想出了问题。调用 invert_a_matrix() 修改传入的矩阵。所以当我调用 gsl_blas_dgemm() 时,我并没有将逆矩阵乘以原始矩阵。
修复是在调用 invert_a_matrix() 之前分配原始矩阵的副本并将副本传递给 gsl_blas_dgemm()。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。