test_chol

// 定义控制台应用程序的入口点。
//


#include "stdafx.h"
#include <mat.h>  
#include <iostream>
#include <iomanip>
#include <Eigen/Dense>
#include <Eigen/LU>
#include <Eigen/Cholesky>
using Eigen::MatrixXd;
using Eigen::VectorXd;
using namespace std;
using namespace Eigen;
void chol(MatrixXd P,double c,MatrixXd &A)
{

LLT<MatrixXd> lltOfP(P);
MatrixXd L=lltOfP.matrixL();        //产生下三角矩阵
A<<c*L;


}
void mat_read(MatrixXd &P,MatrixXd &c,MatrixXd &A)
{


MATFile *pmatFileP = NULL,*pmatFilec = NULL,*pmatFileA = NULL;     
mxArray *pMxArrayP = NULL,*pMxArrayc = NULL,*pMxArrayA = NULL;    
double *inputP,*inputc,*outputA; 




pmatFileP = matOpen("inputP.mat","r");    
pMxArrayP = matGetVariable(pmatFileP,"inputP");    
inputP = (double*) mxGetData(pMxArrayP);  
mwSize M_P = mxGetM(pMxArrayP);  
mwSize N_P = mxGetN(pMxArrayP);    
MatrixXd tempP(M_P,N_P);    
for (int i=0; i<M_P; i++)    
for (int j=0; j<N_P; j++)    
tempP(i,j) = inputP[M_P*j+i];   
P=tempP;
matClose(pmatFileP);    
mxFree(inputP);  
//cout<<"tempP"<<endl<<fixed<<setprecision(5)<<tempP<<endl;


pmatFilec = matOpen("inputc.mat","r");    
pMxArrayc = matGetVariable(pmatFilec,"inputc");    
inputc = (double*) mxGetData(pMxArrayc);    
mwSize M_c = mxGetM(pMxArrayc);  
mwSize N_c = mxGetN(pMxArrayc);    
MatrixXd tempc(M_c,N_c);    
for (int i=0; i<M_c; i++)    
for (int j=0; j<N_c; j++)    
tempc(i,j) = inputc[M_c*j+i];  
c=tempc;
matClose(pmatFilec);    
mxFree(inputc);  
//cout<<"tempc"<<endl<<fixed<<setprecision(5)<<tempc<<endl;


pmatFileA = matOpen("outputA.mat","r");    
pMxArrayA = matGetVariable(pmatFileA,"outputA");    
outputA = (double*) mxGetData(pMxArrayA);  
mwSize M_X = mxGetM(pMxArrayA);  
mwSize N_X = mxGetN(pMxArrayA);    
MatrixXd tempX(M_X,N_X);    
for (int i=0; i<M_X; i++)    
for (int j=0; j<N_X; j++)    
tempX(i,j) = outputA[M_X*j+i];  
A=tempX;
matClose(pmatFileA);    
mxFree(outputA);  
//cout<<"tempX"<<endl<<fixed<<setprecision(5)<<tempX<<endl;




}
void test_chol()
{
MatrixXd P,c,A;
mat_read(P,c,A);
int length=400;
int wideth=60;

MatrixXd inputP(3,3),outputA(3,3);
double inputc;
for(int ut=0;ut<length;++ut)
for(int uk=0;uk<wideth;++uk)
{
int index=ut*wideth+uk;
for(int i=0;i<3;++i)   /从P中提取输入数据
for(int j=0;j<3;++j)
{
inputP(i,j)=P(index*3+i,j);
}
inputc=c(index,0);      //从c中提取输入数据

chol(inputP,inputc,outputA);


/检验outputX与matlab中X的输出的差异，如果大于1e-5，报错
double varience;
double permissibleError=1e-18;
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
{


varience=outputA(i,j)-A(index*3+i,j);
if(varience>=permissibleError)
{
printf("sigmas call error!\n");
printf("(%d,%d)th call error\n",ut,uk);
printf("outputA(%d,%d) error\n",i,j);
printf("outputA=%.6f,X=%.6f,varience=%.6f\n",outputA(i,j),A(index*3+i,j),varience);
return ;
}
}
}
printf("sigmas call success!\n");
}


int _tmain(int argc, _TCHAR* argv[])
{
test_chol();
return 0;
}