当前位置：首页 > news >正文

C/C++语言使用技巧（三十三):强大的数学计算库gsl与算法实现

news 来源：原创 2024/9/21 14:17:11

gsl（GNU Scientific Library）是专门为应用数学和科学技术领域的数值计算提供支持的软件库。gsl使用C语言编写，同时也为其他语言做了相应的封装。gsl在GNU通用公共许可下是免费的。该函数库提供了广泛的数学算法的实现函数，包括随机数生成器，特殊函数和最小二乘拟合等等。目前该函数库提供有超过1000个函数，这些函数包含的范围有：复数计算、多项式求根、特殊函数、向量和矩阵运算、排列、组合、排序、线性代数、特征值和特征向量、快速傅里叶变换(FFT)、数值积分、随机数生成、随机数分布、统计、蒙特卡洛积分、模拟退火、常微分方程组、插值、数值微分、方程求根、最小二乘拟合、小波变换等。更详细的介绍请参考维基百科http://en.wikipedia.org/wiki/GNU_Scientific_Library。
官网：https://www.gnu.org/software/gsl/
教程文档：https://www.gnu.org/software/gsl/doc/html/index.html
下载文档：https://www.gnu.org/software/gsl/doc/latex/gsl-ref.pdf
源码：https://github.com/ampl/gsl
ubuntu一键安装命令：sudo apt install gsl-bin libgsl2 libjson-c2 libgsl0-dev alsa-utils
在这里插入图片描述

dem测试1：矩阵


#include <stdio.h>
#include <gsl/gsl_matrix.h>
int main (void)
{
        int i, j;
        gsl_matrix * m = gsl_matrix_alloc (10, 3);
        for (i = 0; i < 10; i++)
                for (j = 0; j < 3; j++)
                        gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
        for (i = 0; i < 100; i++)
                for (j = 0; j < 3; j++)
                        printf ("m(%d,%d) = %g\n", i, j,
                                        gsl_matrix_get (m, i, j));
        gsl_matrix_free (m);
        return 0;
}

执行结果：

gcc testgsl.c -lm -lgsl -lgslcblas -o main && ./main
m(0,0) = 0.23
m(0,1) = 1.23
m(0,2) = 2.23
m(1,0) = 100.23
m(1,1) = 101.23
m(1,2) = 102.23
m(2,0) = 200.23
m(2,1) = 201.23
m(2,2) = 202.23
m(3,0) = 300.23
m(3,1) = 301.23
m(3,2) = 302.23
m(4,0) = 400.23
m(4,1) = 401.23
m(4,2) = 402.23
m(5,0) = 500.23
m(5,1) = 501.23
m(5,2) = 502.23
m(6,0) = 600.23
m(6,1) = 601.23
m(6,2) = 602.23
m(7,0) = 700.23
m(7,1) = 701.23
m(7,2) = 702.23
m(8,0) = 800.23
m(8,1) = 801.23
m(8,2) = 802.23
m(9,0) = 900.23
m(9,1) = 901.23
m(9,2) = 902.23
gsl: ../gsl/gsl_matrix_double.h:275: ERROR: first index out of range
Default GSL error handler invoked.
Aborted (core dumped)

dem测试2：矩阵计算 ——>逆矩阵

#include <stdio.h>
#include <gsl/gsl_linalg.h>

int main (void)
{
        double a_data[] = { 1.0, 0.6, 0.0,
                0.0, 1.5, 1.0,
                0.0, 1.0, 1.0 };
        /*
         * Inverse is
         *    1  -1.2   1.2
         *    0   2.0  -2.0
         *    0  -2.0   3.0
         */
        double inva[9];
        int s, i, j;

        gsl_matrix_view m
                = gsl_matrix_view_array(a_data, 3, 3);
        gsl_matrix_view inv
                = gsl_matrix_view_array(inva,3,3);
        gsl_permutation * p = gsl_permutation_alloc (3);

        printf("The matrix is\n");
        for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j)
                        printf(j==2?"%6.3f\n":"%6.3f ", gsl_matrix_get(&m.matrix,i,j));

        gsl_linalg_LU_decomp (&m.matrix, p, &s);    
        gsl_linalg_LU_invert (&m.matrix, p, &inv.matrix);

        printf("The inverse is\n");
        for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j)
                        printf(j==2?"%6.3f\n":"%6.3f ",gsl_matrix_get(&inv.matrix,i,j));
        gsl_permutation_free (p);
        return 0;
}

执行结果：

gcc testgsl.c -lm -lgsl -lgslcblas -o main && ./main
The matrix is
 1.000  0.600  0.000
 0.000  1.500  1.000
 0.000  1.000  1.000
The inverse is
 1.000 -1.200  1.200
 0.000  2.000 -2.000
 0.000 -2.000  3.000