前言

推荐两个比较好的教程:

BLAS (Basic Linear Algebra Subprograms)

LAPACK for Windows

命名规范

BLAS基本线性代数子程序的函数命令都有一定规范，便于记忆

<character> <name> <mod> ()

character

定义的是数据类型

也可结合起来，比如sc代表实数域和复数域的单精度类型，dz代表实数域和复数域的双精度类型。

name

针对BLAS level 1，也就是向量与向量间的操作的时候

针对BLAS level 2和3的情况，也就是矩阵参数类型

mod

提供操作的额外信息，三种Level情况

针对Level 1

针对Level 2

针对Level 3

函数实例

C接口的调用方法

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存储到二维数组aa中，矩阵元素

  Aij
  A_{ij}以列优先的方式存储到

  a[i+j∗lda]
  a[i+j*lda]中，或者以行优先的形式存储到

  a[j+i∗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.
意思是联系中的以为数组是将每行列的元素放入到数组单元中。代码明明是每行的元素放入到其中。这到底是按行存储还是按列存储？两篇参考博客介绍行列存储的区别： 行优先和列优先的问题,矩阵存储的两种方式——行优先与列优先。个人感觉，姑且认为是按行存储的吧，毕竟代码是这样写的，虽然与文字不一致。如果有同学有何见解，希望在评论区讨论讨论^_^

Level 1所有函数

所有函数概览

cblas_?asum

• 作用：计算向量元素和
• 定义函数

float cblas_sasum (const MKL_INT n, const float *x, const MKL_INT incx);
float cblas_scasum (const MKL_INT n, const void *x, const MKL_INT incx);
double cblas_dasum (const MKL_INT n, const double *x, const MKL_INT incx);
double cblas_dzasum (const MKL_INT n, const void *x, const MKL_INT incx);

• 运算

  计算实数向量的元素和，或者计算复数向量的实部以及虚部的和
  

  res=|Re(x1)|+|Im(x1)|+|Re(x2|+|Im(x2)|+⋯+|Re(xn)|+|Im(xn)|
  

  res=|Re (x_1)|+|Im(x_1)|+|Re (x_2|+|Im(x_2)|+\\cdots+|Re (x_n)|+|Im(x_n)|

• 输入参数

  n
  n : 向量的元素个数

  xx : 数组，大小至少是

  (1+(n−1)∗abs(incx))
  (1+(n-1)*abs(incx))

  incx
  incx : 指定索引向量

  x
  x的增量

• 返回值:向量所有元素的和

cblas_?axpy

• 作用 : 计算向量-标量的积，然后加到结果上
• 定义函数

void cblas_saxpy (const MKL_INT n, const float a, const float *x, const MKL_INT incx, float *y, const MKL_INT incy);
void cblas_daxpy (const MKL_INT n, const double a, const double *x, const MKL_INT incx, double *y, const MKL_INT incy);
void cblas_caxpy (const MKL_INT n, const void *a, const void *x, const MKL_INT incx, void *y, const MKL_INT incy);
void cblas_zaxpy (const MKL_INT n, const void *a, const void *x, const MKL_INT incx, void *y, const MKL_INT incy);

• 运算

向量与向量之间的操作

y:=a*x+y

• 输入参数

  n
  n : 指定向量x,yx,y的元素个数

  a
  a : 标量aa

  x
  x : 数组，大小至少是(1+(n−1)∗abs(incx))(1+(n-1)*abs(incx))

  incx
  incx : 指定索引向量

  x
  x的增量

  yy : 数组，大小至少是

  (1+(n−1)∗abs(incy))
  (1+(n-1)*abs(incy))

  incy
  incy : 指定索引向量

  y
  y的增量

• 返回值 : 最终更新得到的向量yy

cblas_?copy

• 作用 : 拷贝一个向量到另一个向量

• 定义函数

  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
  y=x

• 输入参数

  n
  n : 指定向量x,yx,y的元素个数

  x
  x : 数组，大小至少是(1+(n−1)∗abs(incx))(1+(n-1)*abs(incx))

  incx
  incx : 指定索引向量

  x
  x的增量

  yy : 数组，大小至少是

  (1+(n−1)∗abs(incy))
  (1+(n-1)*abs(incy))

  incy
  incy : 指定索引向量

  y
  y的增量

• 返回值 : 当nn是正数的时候，向量

  x
  x的拷贝被返回，否则参数不变。

cblas_?dot

• 作用 : 计算向量-向量的点乘

• 定义函数

  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

  res=\\sum_{i=1}^nx_i*y_i

• 输入参数

  n,x,incx,y,incy
  n,x,incx,y,incy分别代表元素个数，数组

  x
  x，数组xx的索引增量，数组

  y
  y，数组yy的索引增量

• 返回值 : 返回两个向量的点乘；如果

  n<0
  n<0，返回0

cblas_?sdot

• 作用 : 计算双精度向量-向量的点乘

• 定义函数

  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
  n : 输入向量xx和

  y
  y的维度

  sbsb : 内积的单精度缩放值(仅针对sdsdot)

  sx,sy
  sx,sy : 数组，包含单精度输入向量

  incx,incy
  incx,incy : 两个数组的索引增量

• 返回值 : 当

  n
  n为正的时候，返回两个数组的点乘(sdsdot结果加一个标量sb);若n≤0n\\leq0，对于sdsdot返回

  sb
  sb，对于dsdot返回0

cblas_?dotc

• 作用 : 计算一个共轭向量与另一个向量的点积

• 定义函数

  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
  

  res=\\sum_{i=1}^n conjg(x_i)*y_i

• 输入参数

  n,x,incx,y,incy
  n,x,incx,y,incy分别代表元素个数，数组

  x
  x及其索引增量，数组yy及其增量

• 输出 : 如果

  n>0
  n>0，输出共轭向量

  x
  x与非共轭向量yy的点乘，否则返回0

cblas_?dotu

• 作用 : 计算复数域的向量-向量的点积

• 定义函数

  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
  

  res=\\sum_{i=1}^n x_i*y_i
  
  其中

  xi
  x_i和

  yi
  y_i分别是复数向量

  x
  x和yy的元素

• 输入参数 : 同上

• 输出参数 : 如果

  n>0
  n>0，返回点积，否则返回0

cblas_?nrm2

• 作用 : 计算一个向量的欧几里得范数(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||
  

  res=||x||

• 输入参数 : 元素个数，数组，索引增量

• 返回值 : 向量

  x
  x的欧几里得范数

cblas_?rot

• 作用 : 平面上绕点旋转

• 定义函数

  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]

  [x_i,y_i]=[x_i,y_i]\\left[
  \\begin{matrix}
  c&-s\\\\
  s&c
  \\end{matrix}
  \\right]
  
  官方写法
  

  xi=c∗xi+s∗yiyi=c∗yi−s∗xi
  

  x_i = c*x_i + s*y_i\\\\
  y_i = c*y_i – s*x_i

• 输入参数 : 元素个数，两个向量，对应的索引增量，

  c,s
  c,s表示所绕点的坐标标量

• 输出参数 :

  x
  x的每个元素被c∗x+s∗yc*x+s*y代替，

  y
  y的每个元素被c∗y−s∗xc*y-s*x代替

cblas_?rotg

• 作用 : 计算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)
  (a,b),返回参数

  c,s,r,z
  c,s,r,z对应Givens旋转，

  c,s
  c,s对应的是酉阵(unitary matrix)，类似于
  

  [c−ssc][ab]=[r0]
  

  \\begin{bmatrix}
  c&s\\\\
  -s&c
  \\end{bmatrix}
  \\begin{bmatrix}
  a\\\\
  b
  \\end{bmatrix}=
  \\begin{bmatrix}
  r\\\\
  0
  \\end{bmatrix}
  
  参数

  z
  z这样定义 : 如果|a|>|b||a|>|b|,那么

  z=s
  z=s，否则如果

  c≠0
  c\\neq 0,那么

  z=1c
  z=\\frac{1}{c}，否则

  z=1
  z=1

• 输入参数 :

  a,b
  a,b分别提供点

  p
  p的横xx坐标和纵

  y
  y坐标

• 输出参数 : Givens的四个参数

cblas_?rotm

• 作用 : 修正绕平面是一点的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,yx,y，这些向量的每个元素都会被替代为
  

  [xiyi]=H[xiyi]
  

  \\begin{bmatrix}
  x_i\\\\
  y_i
  \\end{bmatrix}=
  H\\begin{bmatrix}
  x_i\\\\
  y_i
  \\end{bmatrix}
  
  其中

  i=1⋯n
  i=1\\cdots n，H是修正Givens旋转矩阵，值存储在

  parm[1]
  parm[1]至

  parm[4]
  parm[4]中

• 输入参数

  n,x,incx,y,incy,
  n,x,incx,y,incy,分别表示元素个数，两个向量，对应增量索引

  param : 包含五个参数，

  param[0]
  param[0]代表切换标志，

  param[1−4]
  param[1- 4]均包含

  h11,h12,h21,h22
  h_{11},h_{12},h_{21},h_{22}，分别对应数组

  H
  H的成分，依据标识符，可以得到如下数组HH
  

  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]
  

  flag = -1.0: H =\\begin{bmatrix}
  h_{11} &h_{12}\\\\
  h_{21} &h_{22}
  \\end{bmatrix}\\\\
  flag = 0.0: H =\\begin{bmatrix}
  1.0& h_{12}\\\\
  h_{21} &1.0
  \\end{bmatrix}\\\\
  flag = 1.0: H =\\begin{bmatrix}
  h_{11}& 1.0\\\\
  -1.0 &h_{22}
  \\end{bmatrix}\\\\
  flag = -2.0: H =
  \\begin{bmatrix}
  1.0 &0.0\\\\
  0.0 & 1.0
  \\end{bmatrix}
  
  后三种情况，矩阵

  H
  H中的−0.1,−1.0,0.0-0.1,-1.0,0.0可以由标志指定，无需在向量中写出来

• 输出 :

  x
  x的每一个向量被h11∗x[i]+h12∗y[i]h_{11}*x[i] + h_{12}*y[i] 替换，

  y
  y的每一个向量被h21∗x[i]+h22∗y[i]h_{21}*x[i] + h_{22}*y[i] 替换

cblas_?rotmg

• 作用 : 计算修正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)
  (x_1,y_1),计算将结果向量的y置零时候，计算得到的Givens旋转矩阵

  H
  H
  

  [x10]=H[x1d1−−√y1d2−−√]

  \\begin{bmatrix}
  x_1\\\\
  0
  \\end{bmatrix}=
  H\\begin{bmatrix}
  x_1\\sqrt{d_1}\\\\
  y_1\\sqrt{d_2}
  \\end{bmatrix}

• 输入参数

  d1
  d1代表向量的

  x
  x轴缩放因子，d2d2代表向量的y轴缩放因子,

  x1,y1
  x1,y1是输入向量的横纵坐标

• 输出

  d1
  d1是更新矩阵的第一对角元

  d2
  d2是更新矩阵的第二对角元

  x1
  x1是缩放之前旋转向量的x坐标

  param
  param同上cblas_?rotmg的输入

  parma
  parma一样

cblas_?scal

• 作用 : 计算向量和标量的乘积

• 定义函数

  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
  

  x=a*x

• 输入参数

  n
  n向量的元素个数，aa标量，

  x
  x数组，incxincx向量x的索引增量

• 输出 : 更新后的向量

  x
  x

cblas_?swap

• 作用 : 交换向量值

• 定义函数

  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);

• 计算：交换向量xx和

  y
  y，互相代替元素值

• 输入参数 : 正常的五个输入，元素个数，数组，增量索引

• 输出 : 交换后的向量x,yx,y

cblas_i?amax

• 作用 : 找到绝对值最大的元素的索引

• 定义函数

  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
  x，函数i?amax返回实数域中绝对值最大的向量元素x[i]x[i]的位置，或者返回复数域

  |Re(x[i])|+|Im(x[i])|
  |Re(x[i])|+|Im(x[i])|和的最大值

  如果

  n
  n非正，则返回

  如果向量中有好几个位置的值等于最大元素，第一个位置将被返回

• 输入参数 : n,x,incxn,x,incx分别代表向量元素个数，数组，索引增量

• 返回值 : 返回最大值元素的位置，比如

  x[index−1]
  x[index-1]具有最大的绝对值

cblas_i?amin

• 作用 : 返回最小值的位置

• 定义函数

  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]
  x[index-1] 具有最小的绝对值

cblas_?cabs1

• 作用: 计算复数的绝对值，是一个辅助函数，用于辅助其它函数的实现

• 定义函数

  float cblas_scabs1 (const void *z);
  double cblas_dcabs1 (const void *z);

• 计算
  

  res=|Re(z)|+|Im(z)|
  

  res=|Re(z)|+|Im(z)|

• 输入: 标量

  z
  z

• 返回值: 复数zz的绝对值