Pygsl 的基本使用方法
这几天因为要用 python 处理一些有关矩阵的运算,NumPy 和 SciPy 都很合适。不过,因为我在 C 程序中比较中意 gsl,看到 pythons 有 gsl 的绑定—— pygsl。就尝试了一下。
最开始,让我手足无措的是不知道如何在 pygsl 接口中使用向量或矩阵这样基本的运算元素。我尝试了用 python 的列表,结果出错。认真看了看 pygsl 的文档,才知道目前他们推荐使用 NumPy 的 array 来产生向量和矩阵,然后才可以调用 pygsl 所封装的 gsl 接口。下面,我将 gsl 文档中一个 blas 的示例改写成 pygsl 程序来演示 pygsl 的基本用法。
原 C 程序 gsl-dgemm.c 如下:
#include <stdio.h>
#include <gsl/gsl_blas.h>
int
main (void)
{
double a[] = { 0.11, 0.12, 0.13, 0.21, 0.22, 0.23 };
double b[] = { 1011, 1012, 1021, 1022, 1031, 1032 };
double c[] = { 0.00, 0.00, 0.00, 0.00 };
gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);
/* Compute C = A B */
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
1.0, &A.matrix, &B.matrix,
0.0, &C.matrix);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
return 0;
}
#include <gsl/gsl_blas.h>
int
main (void)
{
double a[] = { 0.11, 0.12, 0.13, 0.21, 0.22, 0.23 };
double b[] = { 1011, 1012, 1021, 1022, 1031, 1032 };
double c[] = { 0.00, 0.00, 0.00, 0.00 };
gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);
/* Compute C = A B */
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
1.0, &A.matrix, &B.matrix,
0.0, &C.matrix);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
return 0;
}
可采用以下方式编译该程序:
$ gcc gsl-dgemm.c -o gsl-dgemm -lgsl -lgslcblas -lm
若采用 pygsl,对应的程序可写为:
#! /usr/bin/env python
# -*- coding:utf-8 -*-
from numpy import *
from pygsl.blas import *
A = array ([[0.11, 0.12, 0.13], [0.21, 0.22, 0.23]])
B = array ([[1011, 1012], [1021, 1022], [1031, 1032]])
C = array ([[0.00, 0.00], [0.00, 0.00]])
C = dgemm (1.0, A, B, 0.0, C)
print (C)
# -*- coding:utf-8 -*-
from numpy import *
from pygsl.blas import *
A = array ([[0.11, 0.12, 0.13], [0.21, 0.22, 0.23]])
B = array ([[1011, 1012], [1021, 1022], [1031, 1032]])
C = array ([[0.00, 0.00], [0.00, 0.00]])
C = dgemm (1.0, A, B, 0.0, C)
print (C)
在 python 交互模式下,可以很方便地查看 pygsl 文档。不过,现在 pygsl 的文档还不是很健全。许多 pygsl 的接口的用法最好还是参照 gsl 的 C 接口文档来理解。