1、用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
2、的MKL了。于是就装个IntelCPlusPlusCompilerProfessional11.1.060，和VS2008集成得很好。啃路径下自带的userguide.pdf配置路径和库文件，然后按照讲解好第一个例子，增强了信心(我最腻味那种为了实现一个功能让我配置半天甚至几天的玩意);然后随着移植算法，啃3000页的mklman.pdf，虽然也实现了矩阵乘除、解线性方程组之类简单功能，但是毕竟都是针对Fortran语言的为主，看起来跌跌撞撞的。连随软件的例子也净是.f的，你这好像是C+编译器啊，Intel大哥-_-!当然咯，我能看懂这些功能，少不了古狗别人的代码示例，这里特别感谢，刚开始起步
3、时看到能用的代码对我这样的菜鸟是在是太重要了)。后面麻烦就来了，就想把用SVD分解来解最小二乘问题的函数dgelsd用起来，就费了劲了，10几个参数mklman.pdf里讲的可以说是乱七八糟，云山雾罩，为了设置一个参数要调另外的函数，而且一会是上文，一会是后面附录的。昨晚看到上午就是没用起来。试图看看clapack里的例子，也因为愚钝没看懂。结果古狗到Intel自己的网站，结果发现美观清晰的C代码，大喜，copy下来一运行，还真用起来了，那些之前不明白的参数一下就懂了。唉后来一看，原来只提供在线版浏览，没随软件安装。唉，这比那个3000页的pdf实用多了!MD，真正的好东西咋不提供pdf啊In
4、tel官网:放到下面other里面了-.-!，我当时也路过过这里，结果看了上面的，也看了下面的，就独独忘了中间这个IntelMathKernelLibraryLAPACKExamplesHTML下面就是我“千辛万苦“(其实不到半天)找到的dgelsd用法。清晰，规范不转帖不足以表示我的感激之情!开头加一下“mkl.h“，再在项目属性的链接器里把几个依赖的lib加上，就可以运行了。/*Copyright(C)2009-2010IntelCorporation.AllRightsReserved.*Theinformationandmaterial(“Material“)providedbelow
5、isownedbyIntel*Corporationoritssuppliersorlicensors,andtitletosuchMaterialremains*withIntelCorporationoritssuppliersorlicensors.TheMaterialcontains*proprietaryinformationofInteloritssuppliersandlicensors.TheMaterial*isprotectedbyworldwidecopyrightlawsandtreatyprovisions.Nopartof*theMaterialmaybecopi
6、ed,reproduced,published,uploaded,posted,*transmitted,ordistributedinanywaywithoutIntelspriorexpresswritten*permission.Nolicenseunderanypatent,copyrightorotherintellectual*propertyrightsintheMaterialisgrantedtoorconferreduponyou,either*expressly,byimplication,inducement,estoppelorotherwise.Anylicense
7、*undersuchintellectualpropertyrightsmustbeexpressandapprovedbyIntel*inwriting.*/*DGELSDExample.=Programcomputestheminimumnorm-solutiontoareallinearleastsquaresproblemusingthesingularvaluedecompositionofA,whereAisthecoefficientmatrix:0.12-8.197.69-2.26-4.71-6.912.22-5.12-9.089.96-3.33-8.94-6.72-4.40-
8、9.983.973.33-2.74-7.92-3.20andBistheright-handsidematrix:7.300.47-6.281.336.58-3.422.68-1.713.46-9.62-0.790.41Description.=Theroutinecomputestheminimum-normsolutiontoareallinearleastsquaresproblem:minimize|b-A*x|usingthesingularvaluedecomposition(SVD)ofA.Aisanm-by-nmatrixwhichmayberank-deficient.Sev
9、eralrighthandsidevectorsbandsolutionvectorsxcanbehandledinasinglecall;theyarestoredasthecolumnsofthem-by-nrhsrighthandsidematrixBandthen-by-nrhssolutionmatrixX.TheeffectiverankofAisdeterminedbytreatingaszerothosesingularvalueswhicharelessthanrcondtimesthelargestsingularvalue.ExampleProgramResults.=D
10、GELSDExampleProgramResultsMinimumnormsolution-0.69-0.240.06-0.80-0.080.210.380.12-0.650.29-0.240.420.290.35-0.30Effectiverank=4Singularvalues18.6615.9910.018.51*/#includestdlib.h#includestdio.h/*DGELSDprototype*/externvoiddgelsd(int*m,int*n,int*nrhs,double*a,int*lda,double*b,int*ldb,double*s,double*
11、rcond,int*rank,double*work,int*lwork,int*iwork,int*info);/*Auxiliaryroutinesprototypes*/externvoidprint_matrix(char*desc,intm,intn,double*a,intlda);/*Parameters*/#defineM4#defineN5#defineNRHS3#defineLDAM#defineLDBN/*Mainprogram*/intmain()/*Locals*/intm=M,n=N,nrhs=NRHS,lda=LDA,ldb=LDB,info,lwork,rank
12、;/*Negativercondmeansusingdefault(machineprecision)value*/doublercond=-1.0;doublewkopt;double*work;/*Localarrays*/*iworkdimensionshouldbeatleast3*min(m,n)*nlvl+11*min(m,n),wherenlvl=max(0,int(log_2(min(m,n)/(smlsiz+1)+1)andsmlsiz=25*/intiwork3*M*0+11*M;doublesM;doubleaLDA*N=0.12,-6.91,-3.33,3.97,-8.
13、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;doublebLDB*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;/*Executablestatements*/printf(“DGELSDExampleProgramResultsn“);/*Queryandallocatetheoptimalworkspace*/lwork=-1;dgel
14、sd(lwork=(int)wkopt;work=(double*)malloc(lwork*sizeof(double);/*SolvetheequationsA*X=B*/dgelsd(/*Checkforconvergence*/if(info0)printf(“ThealgorithmcomputingSVDfailedtoconverge;n“);printf(“theleastsquaressolutioncouldnotbecomputed.n“);exit(1);/*Printminimumnormsolution*/print_matrix(“Minimumnormsolut
15、ion“,n,nrhs,b,ldb);/*Printeffectiverank*/printf(“nEffectiverank=%6in“,rank);/*Printsingularvalues*/print_matrix(“Singularvalues“,1,m,s,1);/*Freeworkspace*/free(void*)work);exit(0);/*EndofDGELSDExample*/*Auxiliaryroutine:printingamatrix*/voidprint_matrix(char*desc,intm,intn,double*a,intlda)inti,j;printf(“n%sn“,desc);for(i=0;im;i+)for(j=0;jn;j+)printf(“%6.2f“,ai+j*lda);printf(“n“);后记:搜英文还得是古狗，哪怕是夹杂大量中文;搜中文还得是百度
