c语言使用英特尔MKL函数库进行SVD分解

英特尔官方下载的mkl说明文档写的过于简略，反而在线的开发人员参考文档中写的很详细，这里指一下路：Singular Value Decomposition (intel.com)

SVD分解就是奇异值分解，具体计算方法可以参考这篇文章：奇异值分解（SVD） - 知乎 (zhihu.com)。

下面直接上代码：

#include <stdio.h>
#include <stdlib.h>
#include "mkl.h"

int main()
{
	/***********************************************************输入参数*******************************************************************/
	int matrix_layout = LAPACK_ROW_MAJOR;//指定矩阵存储布局以行为主还是列为主
	char jobu = 'A';
	/*如果jobu = 'A'，则U的所有m列返回到数组u中；
	如果 jobu = 'S'，则 U 的前 min(m, n) 列（左奇异向量）在数组 u 中返回；
	如果 jobu = 'O'，则 U 的前 min(m, n) 列（左奇异向量）在数组 a 上被覆盖；
	如果 jobu = 'N'，则不计算 U 列（无左奇异向量）*/
	char jobvt = 'A';
	/*如果jobvt = 'A'，则VT/VH的所有n行都返回到数组vt中；
	如果 jobvt = 'S'，则在数组 vt 中返回 VT/VH（右奇异向量）的前 min(m,n) 行；
	如果 jobvt = 'O'，则在数组 a 上覆盖 VT/VH)（右奇异向量）的前 min(m,n) 行；
	如果 jobvt = 'N'，则不计算 VT/VH（无右奇异向量）行 。*/
	lapack_int m = 3;//矩阵行数
	lapack_int n = 2;//矩阵列数
	lapack_int lda = 3;//矩阵行为主则为最大行数，列为主则为最大列数
	double a[9]= { 0,1.0,0,1.0,1.0,0,1.0,0,0 };
	//a（列主布局的大小至少为 max(1, lda*n)，行主布局的大小至少为 max(1, lda*m)）是一个包含 m×n 矩阵 A 的数组。
	lapack_int ldu=3;
	lapack_int ldvt=2;
	/*ldu、ldvt分别是输出数组 u 和 vt 的前导维度。
	约束：
	ldu≥1； ldvt≥1。
	如果jobu = 'A'，ldu≥m；
	如果 jobu = 'S'，则 ldu≥m 为列主布局，ldu≥ min(m, n) 为行主布局；
	如果jobvt = 'A'，ldvt≥n；
	如果 jobvt = 'S'， ldvt≥ min(m, n) 列主要布局和 ldvt≥n 行主要布局。*/


	/***********************************************************输出参数*******************************************************************/
	double s[2];//数组，大小至少为 max(1, min(m,n))。包含经过排序的 A 的奇异值，使得 s[i] ≥ s[i + 1]。
	double u[9];
	/*如果 jobu = 'A'，则 u 包含 m×m 正交 / 酉矩阵 U。
	如果 jobu = 'S'，则 u 包含 U 的前 min(m, n) 列（左奇异向量按列存储）。
	如果 jobu = 'N' 或 'O'，则不引用 u。*/
	double vt[4];
	/*如果 jobvt = 'A'，则包含正交 / 酉矩阵 vtnnVT / VH。
	如果 jobvt = 'S'，则包含 vtmnVT / VH 的前 min(, ) 行（按行存储的右奇异向量）。
	如果 jobvt = 'N' 或 'O'，则不引用 vt。*/
	double superb[9];
	//如果 ?bdsqr 不收敛（返回值信息 > 0 ），则在退出时 super(0:min(m,n)-2) 包含上对角矩阵的未收敛超对角元素，其对角线在）。满足 BsBA = u*B*VT或 A = u*B*VH，因此它具有相同的奇异值，以及和 .Auvt 关联的奇异向量
	int info = LAPACKE_dgesvd(matrix_layout, jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, superb);
	//getchar();

	printf("U矩阵：\n");
	for (int i = 0; i < 9; i++) {
		printf( "%6.4lf ", u[i]);
	}
	printf( "\n");
	printf("VT矩阵：\n");
	for (int i = 0; i < 4; i++) {
		printf( "%6.4lf ", vt[i]);
	}
	printf( "\n");
	printf("奇异值：\n");
	for (int i = 0; i < 2; i++) {
		printf( "%6.4lf ", s[i]);
	}
	
	return 0;
}

上面代码中计算的矩阵是文章开头的知乎教程链接中的矩阵，即：

首先在变量matrix_layout中指定矩阵是按列存储还是按行存储，这个会直接影响到输入矩阵和输出矩阵的存储格式。
比如本文使用的是按行存储，而MKL要求输入的a矩阵为一个包含A矩阵的方阵，因此需要声明a矩阵为：a[9]= { 0,1,0,1,1,0,1,0,0 };（从第一行开始顺序输入所有行）即：

若是按列存储，需要声明a矩阵为：a[9]= { 0,1,1,1,1,0 };（从第一列开始顺序输入前两列）即可。
因为我需要返回完整的U矩阵和VT矩阵，因此我指定jobu = ‘A’;和 jobvt = ‘A’;并指定两个矩阵的维数：double u[9]; double vt[4];
MKL的svd函数不会直接生成一个3×2的奇异值矩阵，因此只需声明s[2];用来存储两个奇异值。
测试结果：

来和matlab计算结果对比下：




matlab做数学计算真的方便啊呜呜呜