matrix - 对 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.0s 和其他地方都有 0.0s 的矩阵)。所以我认为作为上述 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, p, 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, i, 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, CblasNoTrans, 1.0, mat, inverse, 0.0, product);
// error = gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inverse, mat, 0.0, 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()。
推荐阅读
- javascript - 从两个数组中过滤重复项
- google-cloud-storage - @google-cloud/storage 使用 generateSignedPostPolicyV4 出错
- c# - C# Selenium 检查浏览器是否处于活动状态/焦点
- python - 在 Flask 中使用 LOGIN_USER() 变得很棘手
- javascript - 每个动态模式的唯一评论部分
- linux - 在 Ubuntu 19.10 或 20.04 上安装 Tizen Studio 3.7?
- cocoa - tvOS applicationDidBecomeActive 在电视关闭的情况下反复调用
- java - 在 TabLayout + ViewPager2 中选择第三页作为默认值
- sql - 在 oracle sql 中查找某个日期之前的产品订单的查询是什么?
- node.js - 如何在不同 .js 文件中的函数之间读取、更新和传递变量?