对 gsl_blas_dgemm() 的调用未按预期工作

Question

我正在研究一个需要求解 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

关于我做错了什么的任何线索？

提前致谢

Answer 1

我发现了问题所在。对 invert_a_matrix() 的调用会修改传入的矩阵。因此，当我调用 gsl_blas_dgemm() 时，我并没有将逆矩阵乘以原始矩阵。

修复是在调用 invert_a_matrix() 之前分配原始矩阵的副本并将副本传递给 gsl_blas_dgemm()。