用Intel MKL解SVD和最小二乘
2010年4月2日星期五
这几天赶着移植算法，要把MATLAB程序用C++改写。涉及6000多维的方阵运算，幸亏笔记本上2008R2是64位的，64位MATLAB还能算出数来，32位电脑上的MATLAB直接就存储空间不足罢工了。
研究了几个C++下矩阵运算的库，boost带的BLAS太矬，lapack，clapack，cpplapack全是在Linux上的设置多，讲windows下配置的太少，都得用Makefile慢慢编译……搞科学计算的牛人都是用Linux，可对我这种半路出家的矬人实在是过于高深，再说我也懒得去研究这些，只想找个立刻就能用起来的。那么就Intel的MKL了。于是就装个Intel C Plus Plus Compiler Professional 11.1.060，和VS2008集成得很好。
啃路径下自带的userguide.pdf配置路径和库文件，然后按照讲解好第一个例子，增强了信心（我最腻味那种为了实现一个功能让我配置半天甚至几天的玩意）；
然后随着移植算法，啃3000页的mklman.pdf，虽然也实现了矩阵乘除、解线性方程组之类简单功能，但是毕竟都是针对Fortran语言的为主，看起来跌跌撞撞的。连随软件的例子也净是.f的，你这好像是C++编译器啊，Intel大哥-_-!
当然咯，我能看懂这些功能，少不了古狗别人的代码示例，这里特别感谢http://blog.sina.com.cn/s/blog_4840fe2a0100gjxz.html，刚开始起步时看到能用的代码对我这样的菜鸟是在是太重要了）。
后面麻烦就来了，就想把用SVD分解来解最小二乘问题的函数dgelsd用起来，就费了劲了，10几个参数mklman.pdf里讲的可以说是乱七八糟，云山雾罩，为了设置一个参数要调另外的函数，而且一会是上文，一会是后面附录的。昨晚看到上午就是没用起来。试图看看clapack里的例子，也因为愚钝没看懂。结果古狗到Intel自己的网站，结果发现美观清晰的C代码，大喜，copy下来一运行，还真用起来了，那些之前不明白的参数一下就懂了。唉后来一看，原来只提供在线版浏览，没随软件安装。唉，这比那个3000页的pdf实用多了！MD，真正的好东西咋不提供pdf啊……
Intel官网：
http://software.intel.com/en-us/articles/intel-math-kernel-library-documentation/
放到下面other里面了-.-!，我当时也路过过这里，结果看了上面的，也看了下面的，就独独忘了中间这个……
Intel Math Kernel Library LAPACK Examples [HTML]
下面就是我“千辛万苦”（其实不到半天）找到的dgelsd用法。清晰，规范……不转帖不足以表示我的感激之情！
http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/index.htm
开头加一下”mkl.h”，再在项目属性的链接器里把几个依赖的lib加上，就可以运行了。
/*******************************************************************************
* Copyright (C) 2009-2010 Intel Corporation. All Rights Reserved.
* The information and material ("Material") provided below is owned by Intel
* Corporation or its suppliers or licensors, and title to such Material remains
* with Intel Corporation or its suppliers or licensors. The Material contains
* proprietary information of Intel or its suppliers and licensors. The Material
* is protected by worldwide copyright laws and treaty provisions. No part of
* the Material may be copied, reproduced, published, uploaded, posted,
* transmitted, or distributed in any way without Intel’s prior express written
* permission. No license under any patent, copyright or other intellectual
* property rights in the Material is granted to or conferred upon you, either
* expressly, by implication, inducement, estoppel or otherwise. Any license
* under such intellectual property rights must be express and approved by Intel
* in writing.
*
********************************************************************************
*/
/*
DGELSD Example.
==============
Program computes the minimum norm-solution to a real linear least squares
problem using the singular value decomposition of A,
where A is the coefficient matrix:
0.12 -8.19 7.69 -2.26 -4.71
-6.91 2.22 -5.12 -9.08 9.96
-3.33 -8.94 -6.72 -4.40 -9.98
3.97 3.33 -2.74 -7.92 -3.20
and B is the right-hand side matrix:
7.30 0.47 -6.28
1.33 6.58 -3.42
2.68 -1.71 3.46
-9.62 -0.79 0.41
Description.
============
The routine computes the minimum-norm solution to a real linear least
squares problem: minimize ||b – A*x|| using the singular value
decomposition (SVD) of A. A is an m-by-n matrix which may be rank-deficient.
Several right hand side vectors b and solution vectors x can be handled
in a single call; they are stored as the columns of the m-by-nrhs right
hand side matrix B and the n-by-nrhs solution matrix X.
The effective rank of A is determined by treating as zero those singular
values which are less than rcond times the largest singular value.
Example Program Results.
========================
DGELSD Example Program Results
Minimum norm solution
-0.69 -0.24 0.06
-0.80 -0.08 0.21
0.38 0.12 -0.65
0.29 -0.24 0.42
0.29 0.35 -0.30
Effective rank = 4
Singular values
18.66 15.99 10.01 8.51
*/
#include <stdlib.h>
#include <stdio.h>
/* DGELSD prototype */
extern void dgelsd( int* m, int* n, int* nrhs, double* a, int* lda,
double* b, int* ldb, double* s, double* rcond, int* rank,
double* work, int* lwork, int* iwork, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, double* a, int lda );
/* Parameters */
#define M 4
#define N 5
#define NRHS 3
#define LDA M
#define LDB N
/* Main program */
int main() {
/* Locals */
int m = M, n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info, lwork, rank;
/* Negative rcond means using default (machine precision) value */
double rcond = -1.0;
double wkopt;
double* work;
/* Local arrays */
/* iwork dimension should be at least 3*min(m,n)*nlvl + 11*min(m,n),
where nlvl = max( 0, int( log_2( min(m,n)/(smlsiz+1) ) )+1 )
and smlsiz = 25 */
int iwork[3*M*0+11*M];
double s[M];
double a[LDA*N] = {
0.12, -6.91, -3.33, 3.97,
-8.19, 2.22, -8.94, 3.33,
7.69, -5.12, -6.72, -2.74,
-2.26, -9.08, -4.40, -7.92,
-4.71, 9.96, -9.98, -3.20
};
double b[LDB*NRHS] = {
7.30, 1.33, 2.68, -9.62, 0.00,
0.47, 6.58, -1.71, -0.79, 0.00,
-6.28, -3.42, 3.46, 0.41, 0.00
};
/* Executable statements */
printf( " DGELSD Example Program Results\n" );
/* Query and allocate the optimal workspace */
lwork = -1;
dgelsd( &m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, &wkopt, &lwork,
iwork, &info );
lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );
/* Solve the equations A*X = B */
dgelsd( &m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork,
iwork, &info );
/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm computing SVD failed to converge;\n" );
printf( "the least squares solution could not be computed.\n" );
exit( 1 );
}
/* Print minimum norm solution */
print_matrix( "Minimum norm solution", n, nrhs, b, ldb );
/* Print effective rank */
printf( "\n Effective rank = %6i\n", rank );
/* Print singular values */
print_matrix( "Singular values", 1, m, s, 1 );
/* Free workspace */
free( (void*)work );
exit( 0 );
} /* End of DGELSD Example */
/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
int i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
printf( "\n" );
}
}