推荐两个比较好的教程:
BLAS (Basic Linear Algebra Subprograms)
LAPACK for Windows
BLAS基本线性代数子程序的函数命令都有一定规范,便于记忆
<character> <name> <mod> ()定义的是数据类型
s实数域,单精度c复数域,单精度d实数域,双精度z复数域,双精度也可结合起来,比如sc代表实数域和复数域的单精度类型,dz代表实数域和复数域的双精度类型。
针对BLAS level 1,也就是向量与向量间的操作的时候
dot向量点乘rot向量旋转swap向量交换针对BLAS level 2和3的情况,也就是矩阵参数类型
ge一般矩阵gb一般带状矩阵sy对称矩阵sp对称矩阵(压缩存储)sb对称带状矩阵he埃尔米特矩阵 Hermitian matrixhp埃尔米特矩阵(压缩存储)hb埃尔米特带状矩阵tr三角矩阵tp三角矩阵(压缩存储)tb三角带状矩阵提供操作的额外信息,三种Level情况
针对Level 1
c共轭向量u非共轭向量gGivens 旋转结构m修正 Givens 旋转mg修正 Givens 旋转结构针对Level 2
mv矩阵-向量乘积sv求解具有一个未知向量的线性方程组r矩阵的一阶更新r2矩阵的二阶更新针对Level 3
mm矩阵-矩阵乘积sm求解具有多个未知向量的线性方程组rk矩阵的k阶更新r2k矩阵的2k阶更新C接口相对于Fortran有一个优势就是,可以指定是行优先还是列优先。调用方法就是在正常的命令规范前面加个前缀cblas_。对于复数函数?dotc和?dotu还需要加一个后缀_sub,通过指针返回复数结果,作为后一个参数添加进来。
使用CALBAS需要遵循一些规则 :
输入参数需要使用const修饰非复数标量输入参数直接使用值传递复数标量参数使用void指针传递矩阵参数使用地址传递BLAS类型参数使用合适的枚举类型代替Level 2和3中有一个CLABS_LAYOUT类型的额外参数,作为第一个输入参数。这个参数指定二维数组是行优先CblasRowMajor还是列优先CblasColMajor定义的枚举类型:
enum CBLAS_LAYOUT { CblasRowMajor=101, /* row-major arrays */ CblasColMajor=102}; /* column-major arrays */ enum CBLAS_TRANSPOSE { CblasNoTrans=111, /* trans='N' */ CblasTrans=112, /* trans='T' */ CblasConjTrans=113}; /* trans='C' */ enum CBLAS_UPLO { CblasUpper=121, /* uplo ='U' */ CblasLower=122}; /* uplo ='L' */ enum CBLAS_DIAG { CblasNonUnit=131, /* diag ='N' */ CblasUnit=132}; /* diag ='U' */ enum CBLAS_SIDE { CblasLeft=141, /* side ='L' */ CblasRight=142}; /* side ='R' */三种存储方案:
全部存储:比如将矩阵 A 存储到二维数组a中,矩阵元素 Aij 以列优先的方式存储到 a[i+j∗lda] 中,或者以行优先的形式存储到 a[j+i∗lda] 中。这里的lda就是数组的引导维度。压缩存储:用于存储对称阵,埃尔米特矩阵,三角矩阵。对于列优先的布局,上三角和下三角按照列存储到一个一维数组;或者按照行优先的布局,上三角和下三角按行存储到一个一维数组中。带状存储:带状矩阵被压缩存储到一个二维数组中。对于列优先布局,矩阵的列被存储到对应的数组列中,矩阵对角部分被存储到数组的特殊行中。对于行优先的布局,矩阵的行被存储到对应数组的行中,矩阵对角部分被存储到数组的指定行中。Fortran中以列优先的方式存储二维数组;在C中,需要数组为行优先的格式,这是C的约定。讲道理的话,这句话应该是这样写,但是官方文档给出的英文是 The BLAS routines follow the Fortran convention of storing two-dimensional arrays using column-major layout. When calling BLAS routines from C, remember that they require arrays to be in column-major format, not the row-major format that is the convention for C. Unless otherwise specified, the psuedo-code examples for the BLAS routines illustrate matrices stored using column-major layout. 这里却写着是按照列存储的。但是在Intel® Math Kernel Library Getting Started Tutorial: Using the Intel® Math Kernel Library for Matrix Multiplication中的实例却是这样 代码是 这个for循环明显是将一行一行的赋值,说明连续存储的单元是行中相邻的元素。但是后面的英文又是 The one-dimensional arrays in the exercises store the matrices by placing the elements of each column in successive cells of the arrays. 意思是联系中的以为数组是将每行列的元素放入到数组单元中。代码明明是每行的元素放入到其中。这到底是按行存储还是按列存储?两篇参考博客介绍行列存储的区别: 行优先和列优先的问题,矩阵存储的两种方式——行优先与列优先。个人感觉,姑且认为是按行存储的吧,毕竟代码是这样写的,虽然与文字不一致。如果有同学有何见解,希望在评论区讨论讨论^_^
运算
计算实数向量的元素和,或者计算复数向量的实部以及虚部的和
res=|Re(x1)|+|Im(x1)|+|Re(x2|+|Im(x2)|+⋯+|Re(xn)|+|Im(xn)|输入参数
n : 向量的元素个数
x : 数组,大小至少是 (1+(n−1)∗abs(incx))
incx : 指定索引向量 x 的增量
返回值:向量所有元素的和
向量与向量之间的操作 y:=a∗x y
输入参数
n : 指定向量x,y的元素个数
a : 标量a
x : 数组,大小至少是(1 (n−1)∗abs(incx))
incx : 指定索引向量 x 的增量
y : 数组,大小至少是 (1+(n−1)∗abs(incy))
incy : 指定索引向量 y 的增量
返回值 : 最终更新得到的向量y
作用 : 拷贝一个向量到另一个向量
定义函数
void cblas_scopy (const MKL_INT n, const float *x, const MKL_INT incx, float *y, const MKL_INT incy); void cblas_dcopy (const MKL_INT n, const double *x, const MKL_INT incx, double *y, const MKL_INT incy); void cblas_ccopy (const MKL_INT n, const void *x, const MKL_INT incx, void *y, const MKL_INT incy); void cblas_zcopy (const MKL_INT n, const void *x, const MKL_INT incx, void *y, const MKL_INT incy);运算
拷贝向量 y=x
输入参数
n : 指定向量x,y的元素个数
x : 数组,大小至少是(1 (n−1)∗abs(incx))
incx : 指定索引向量 x 的增量
y : 数组,大小至少是 (1+(n−1)∗abs(incy))
incy : 指定索引向量 y 的增量
返回值 : 当n是正数的时候,向量 x 的拷贝被返回,否则参数不变。
作用 : 计算向量-向量的点乘
定义函数
float cblas_sdot (const MKL_INT n, const float *x, const MKL_INT incx, const float *y, const MKL_INT incy); double cblas_ddot (const MKL_INT n, const double *x, const MKL_INT incx, const double *y, const MKL_INT incy);运算 : res=∑i=1nxi∗yi
输入参数
n,x,incx,y,incy 分别代表元素个数,数组 x ,数组x的索引增量,数组 y ,数组y的索引增量
返回值 : 返回两个向量的点乘;如果 n<0 ,返回0
作用 : 计算双精度向量-向量的点乘
定义函数
float cblas_sdsdot (const MKL_INT n, const float sb, const float *sx, const MKL_INT incx, const float *sy, const MKL_INT incy); double cblas_dsdot (const MKL_INT n, const float *sx, const MKL_INT incx, const float *sy, const MKL_INT incy);运算
?sdot计算的是两个双精度向量的内积(点积),中间结果的累积是双精度的,但是sdsdot返回的结果是单精度的,使用dsdot可以输出双精度的结果。其中sdsdot为点积结果加一个标量值sb
输入参数
n : 输入向量x和 y 的维度
sb : 内积的单精度缩放值(仅针对sdsdot)
sx,sy : 数组,包含单精度输入向量
incx,incy : 两个数组的索引增量
返回值 : 当 n 为正的时候,返回两个数组的点乘(sdsdot结果加一个标量sb);若n≤0,对于sdsdot返回 sb ,对于dsdot返回0
作用 : 计算一个共轭向量与另一个向量的点积
定义函数
void cblas_cdotc_sub (const MKL_INT n, const void *x, const MKL_INT incx, const void *y, const MKL_INT incy, void *dotc); void cblas_zdotc_sub (const MKL_INT n, const void *x, const MKL_INT incx, const void *y, const MKL_INT incy, void *dotc);运算
res=∑i=1nconjg(xi)∗yi输入参数
n,x,incx,y,incy 分别代表元素个数,数组 x 及其索引增量,数组y及其增量
输出 : 如果 n>0 ,输出共轭向量 x 与非共轭向量y的点乘,否则返回0
作用 : 计算复数域的向量-向量的点积
定义函数
void cblas_cdotu_sub (const MKL_INT n, const void *x, const MKL_INT incx, const void *y, const MKL_INT incy, void *dotu); void cblas_zdotu_sub (const MKL_INT n, const void *x, const MKL_INT incx, const void *y, const MKL_INT incy, void *dotu);运算
res=∑i=1nxi∗yi 其中 xi 和 yi 分别是复数向量 x 和y的元素输入参数 : 同上
输出参数 : 如果 n>0 ,返回点积,否则返回0
作用 : 计算一个向量的欧几里得范数(Euclidean norm)
定义函数
float cblas_snrm2 (const MKL_INT n, const float *x, const MKL_INT incx); double cblas_dnrm2 (const MKL_INT n, const double *x, const MKL_INT incx); float cblas_scnrm2 (const MKL_INT n, const void *x, const MKL_INT incx); double cblas_dznrm2 (const MKL_INT n, const void *x, const MKL_INT incx);运算
res=||x||输入参数 : 元素个数,数组,索引增量
返回值 : 向量 x 的欧几里得范数
作用 : 平面上绕点旋转
定义函数
void cblas_srot (const MKL_INT n, float *x, const MKL_INT incx, float *y, const MKL_INT incy, const float c, const float s); void cblas_drot (const MKL_INT n, double *x, const MKL_INT incx, double *y, const MKL_INT incy, const double c, const double s); void cblas_csrot (const MKL_INT n, void *x, const MKL_INT incx, void *y, const MKL_INT incy, const float c, const float s); void cblas_zdrot (const MKL_INT n, void *x, const MKL_INT incx, void *y, const MKL_INT incy, const double c, const double s);运算
通俗写法 [xi,yi]=[xi,yi][cs−sc] 官方写法
xi=c∗xi+s∗yiyi=c∗yi−s∗xi输入参数 : 元素个数,两个向量,对应的索引增量, c,s 表示所绕点的坐标标量
输出参数 : x 的每个元素被c∗x s∗y代替, y 的每个元素被c∗y−s∗x代替
作用 : 计算Givens旋转参数
定义函数
void cblas_srotg (float *a, float *b, float *c, float *s); void cblas_drotg (double *a, double *b, double *c, double *s); void cblas_crotg (void *a, const void *b, float *c, void *s); void cblas_zrotg (void *a, const void *b, double *c, void *s);运算
给定一个点的笛卡尔(Cartesian)坐标 (a,b) ,返回参数 c,s,r,z 对应Givens旋转, c,s 对应的是酉阵(unitary matrix),类似于
[c−ssc][ab]=[r0] 参数 z 这样定义 : 如果|a|>|b|,那么 z=s ,否则如果 c≠0 ,那么 z=1c ,否则 z=1输入参数 : a,b 分别提供点 p 的横x坐标和纵 y 坐标
输出参数 : Givens的四个参数
作用 : 修正绕平面是一点的Givens旋转
定义函数
void cblas_srotm (const MKL_INT n, float *x, const MKL_INT incx, float *y, const MKL_INT incy, const float *param); void cblas_drotm (const MKL_INT n, double *x, const MKL_INT incx, double *y, const MKL_INT incy, const double *param);运算
给定两个向量x,y,这些向量的每个元素都会被替代为
[xiyi]=H[xiyi] 其中 i=1⋯n ,H是修正Givens旋转矩阵,值存储在 parm[1] 至 parm[4] 中输入参数
n,x,incx,y,incy, 分别表示元素个数,两个向量,对应增量索引
param : 包含五个参数, param[0] 代表切换标志, param[1−4] 均包含 h11,h12,h21,h22 ,分别对应数组 H 的成分,依据标识符,可以得到如下数组H
flag=−1.0:H=[h11h21h12h22]flag=0.0:H=[1.0h21h121.0]flag=1.0:H=[h11−1.01.0h22]flag=−2.0:H=[1.00.00.01.0] 后三种情况,矩阵 H 中的−0.1,−1.0,0.0可以由标志指定,无需在向量中写出来输出 : x 的每一个向量被h11∗x[i] h12∗y[i]替换, y 的每一个向量被h21∗x[i] h22∗y[i]替换
作用 : 计算修正Givens旋转的参数
定义函数
void cblas_srotmg (float *d1, float *d2, float *x1, const float y1, float *param); void cblas_drotmg (double *d1, double *d2, double *x1, const double y1, double *param);运算
给定输入向量的笛卡尔坐标 (x1,y1) ,计算将结果向量的y置零时候,计算得到的Givens旋转矩阵 H [x10]=H[x1d1−−√y1d2−−√]
输入参数
d1 代表向量的 x 轴缩放因子,d2代表向量的y轴缩放因子, x1,y1 是输入向量的横纵坐标
输出
d1 是更新矩阵的第一对角元
d2 是更新矩阵的第二对角元
x1 是缩放之前旋转向量的x坐标
param 同上cblas_?rotmg的输入 parma 一样
作用 : 计算向量和标量的乘积
定义函数
void cblas_sscal (const MKL_INT n, const float a, float *x, const MKL_INT incx); void cblas_dscal (const MKL_INT n, const double a, double *x, const MKL_INT incx); void cblas_cscal (const MKL_INT n, const void *a, void *x, const MKL_INT incx); void cblas_zscal (const MKL_INT n, const void *a, void *x, const MKL_INT incx); void cblas_csscal (const MKL_INT n, const float a, void *x, const MKL_INT incx); void cblas_zdscal (const MKL_INT n, const double a, void *x, const MKL_INT incx);运算
x=a∗x输入参数
n 向量的元素个数,a标量, x 数组,incx向量x的索引增量
输出 : 更新后的向量 x
作用 : 交换向量值
定义函数
void cblas_sswap (const MKL_INT n, float *x, const MKL_INT incx, float *y, const MKL_INT incy); void cblas_dswap (const MKL_INT n, double *x, const MKL_INT incx, double *y, const MKL_INT incy); void cblas_cswap (const MKL_INT n, void *x, const MKL_INT incx, void *y, const MKL_INT incy); void cblas_zswap (const MKL_INT n, void *x, const MKL_INT incx, void *y, const MKL_INT incy);计算:交换向量x和 y ,互相代替元素值
输入参数 : 正常的五个输入,元素个数,数组,增量索引
输出 : 交换后的向量x,y
作用 : 找到绝对值最大的元素的索引
定义函数
CBLAS_INDEX cblas_isamax (const MKL_INT n, const float *x, const MKL_INT incx); CBLAS_INDEX cblas_idamax (const MKL_INT n, const double *x, const MKL_INT incx); CBLAS_INDEX cblas_icamax (const MKL_INT n, const void *x, const MKL_INT incx); CBLAS_INDEX cblas_izamax (const MKL_INT n, const void *x, const MKL_INT incx);计算 : 给定向量 x ,函数i?amax返回实数域中绝对值最大的向量元素x[i]的位置,或者返回复数域 |Re(x[i])|+|Im(x[i])| 和的最大值
如果 n 非正,则返回
如果向量中有好几个位置的值等于最大元素,第一个位置将被返回
输入参数 : n,x,incx分别代表向量元素个数,数组,索引增量
返回值 : 返回最大值元素的位置,比如 x[index−1] 具有最大的绝对值
作用 : 返回最小值的位置
定义函数
CBLAS_INDEX cblas_isamin (const MKL_INT n, const float *x, const MKL_INT incx); CBLAS_INDEX cblas_idamin (const MKL_INT n, const double *x, const MKL_INT incx); CBLAS_INDEX cblas_icamin (const MKL_INT n, const void *x, const MKL_INT incx); CBLAS_INDEX cblas_izamin (const MKL_INT n, const void *x, const MKL_INT incx);计算 : i?amin返回的是向量的最小绝对值的位置。与i?amax类似
输入参数: 同i?amax
返回值: 绝对值最小元素位置的索引,比如 x[index−1] 具有最小的绝对值
作用: 计算复数的绝对值,是一个辅助函数,用于辅助其它函数的实现
定义函数
float cblas_scabs1 (const void *z); double cblas_dcabs1 (const void *z);计算
res=|Re(z)|+|Im(z)|输入: 标量 z
返回值: 复数z的绝对值
