二次规划:有效集法(Active-Set)

参考参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法(Active-Set),有效集法只能优化中等规模的最优化问题,如果是大规模的最优化问题,应该采用其它算法,代码如下:(代码执行环境VS2010)

QP.h
#ifndef _QP_
#define _QP_
#pragma comment (lib,"./bin/blasd.lib")
#pragma comment (lib,"./bin/clapackd.lib")
#pragma comment (lib,"./bin/libf2cd.lib")
#pragma comment (lib,"./bin/tmglibd.lib")

#include <fstream>
#include <assert.h>
#include <vector>
#include <algorithm>
#include <iostream>

#include "Lapack.h"
using namespace std;

/*
Copy right 
用拉格朗日乘子发求解等式约束的二次规划问题
求解问题为:min: 0.5*trans(x)*H*x + trans(c)*x; st: A*x = b

参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法(Active-Set)
*/
class QuadEquLag //等式约束问题
{
private:
	//
	double** rTmp;
public:
	double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b);//返回最优解
	void calculateQ(double** reverseH, int& row, int& col, double** A,double** Q);
	void calculateR(double** reverseH, int& row, int& col, double** A,double** R);
	void twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out);
	void inverseMatrix(double** a, int& aRow, int& aCol); //计算矩阵a的逆矩阵
	
	void testMatrixInverse();	
};

class QuadProg : public QuadEquLag //不等式约束问题
{
private:

public:
	//double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0);
	double optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0); //有效集算法只能处理中等规模的数据
	int rank(double** a, int& row, int& col, double tol=1e-10); //用SVD分解计算矩阵的秩
	void solveLinearEquations(double** a, int& row, int& col, double*);
	bool greaterThanZero(double* in, int& n);
	bool equalZero(double* in, int& n);
	void updateActiveSet(vector<int>& activeSet, double* mu);
	double findAfak(vector<int>& activeSet, double* dk, double** A , double* b,double* xk, int& minIndex, int& m, int& n);
	void updataXk(double* xk,double& afak, double* dk, int& n);
};

#endif

 

#include "QP.h"

void QuadEquLag::testMatrixInverse()
{
	int N = 500;
	double* A = new double[9];
	
	int LDA = 500;
	int* IPIV = new int[N];
	int LWORK = N;
	double* WORK = new double[N];
	int INFO ;
	ifstream readin;
	readin.open("./data/m.txt",ios_base::in);
	int k = 0;
	if(readin.is_open())
	{
		for(int i = 0; i < N; i++)
		{
			for(int j = 0; j < N; j++)
			{
				double tmp = 0;
				readin>>tmp;
				A[k++] = tmp;
			}
		}
	}
	int m = N; int n = m; 
	dgetrf_(&m, &n, A, &m, IPIV, &INFO);
	dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
}

QP.CPP
double* QuadEquLag::optimizer(double** H,int& m, int& n, double* c, double** A, double* b)
{
	//H为n维对称阵,A为m*n矩阵,c为n维列向量,b为m维列向量
	double* res = new double[n];

	inverseMatrix(H, n, n);
	double** R = new double*[n];
	for(int i = 0; i < n; i++)
	{
		R[i] = new double[m];
	}
	
	calculateR(H, m, n, A,R);
	rTmp = R;

	double** Q = new double*[n];
	for(int i = 0; i < n; i++)
	{
		Q[i] = new double[n];
	}
	calculateQ(H, m, n, A,Q);
	
	double* Qc = new double[n];
	for(int i = 0; i < n; i++)
	{
		double tmp1 = 0;
		double tmp2 = 0;
		for(int k1 = 0; k1 < n; k1++)
		{
			tmp1 = tmp1 + Q[i][k1]*c[k1];
		}

		for(int k2 = 0; k2 < m; k2++)
		{
			tmp2 = tmp2 + R[i][k2]*b[k2];
		}

		res[i] = -tmp1 + tmp2;
	}

	for(int i = 0; i < n; i++)
	{
		delete[] R[i];
	}
	delete[] R;
	for(int i = 0; i < n; i++)
	{
		delete[] Q[i];
	}
	delete[] Q;
	return res;
}
void QuadEquLag::calculateQ(double** reverseH, int& m, int& n, double** A,double** Q)
{
	//H为n维对称阵,A为m*n矩阵,c为n维列向量,b为m维列向量
	================计算H逆乘A转置=====================
	//double** rHtA = new double*[n]; //n*m维
	//for(int i = 0; i < n; i++)
	//{
	//	rHtA[i] = new double[m];
	//}
	//twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
	=====================================================
	===============计算A乘H逆乘A转置=====================
	//double** ArH = new double*[m]; //m*n维
	//for(int i = 0; i < m; i++)
	//{
	//	ArH[i] = new double[n];
	//}
	//twoMatrixMultiply(A,m,n,"N",reverseH,m,n,"N",ArH);
	//double** ArHtA = new double*[m];
	//for(int i = 0; i < m; i++)
	//{
	//	ArHtA[i] = new double[m];
	//}
	//twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
	//inverseMatrix(ArHtA, m, m);
	释放内存空间
	//for(int i = 0; i < m; i++)
	//{
	//	delete[] ArH[i];
	//}
	//delete[] ArH;
	//double** tmp1 = new double*[n];
	//for(int i = 0; i < n; i++)
	//{
	//	tmp1[i] = new double[m];
	//}
	//calculateR(reverseH, m, n, A,tmp1);
	//=====================================================
	//=====================计算A乘H逆======================
	double** AtH = new double*[m];
	for(int i = 0; i < m; i++)
	{
		AtH[i] = new double[n];
	}
	twoMatrixMultiply(A,m,n,"N", reverseH,n,n,"N",AtH);
	//=====================================================
	释放内存
	//for(int i = 0; i < n; i++)
	//{
	//	delete[] rHtA[i];
	//}
	//delete[] rHtA;
	//for(int i = 0; i < m; i++)
	//{
	//	delete[] ArHtA[i];
	//}
	//delete[] ArHtA;

	double** tmp2 = new double*[n];
	for(int i = 0; i < n; i++)
	{
		tmp2[i] = new double[n];
	}
	twoMatrixMultiply(rTmp,n,m,"N", AtH,m,n,"N",tmp2);
	//释放内存空间
	/*for(int i =0; i < n; i++)
	{
		delete[] tmp1[i];
	}
	delete[] tmp1;*/
	for(int i =0; i < m; i++)
	{
		delete[] AtH[i];
	}
	delete[] AtH;
	//两个矩阵相减
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			Q[i][j] = reverseH[i][j] - tmp2[i][j];
		}
	}

	//=================内存清理==================
	for(int i =0; i < n; i++)
	{
		delete[] tmp2[i];
	}
	delete[] tmp2;
	//===========================================
}

void QuadEquLag::calculateR(double** reverseH, int& m, int& n, double** A,double** R)
{
	//H为n维对称阵,A为m*n矩阵,c为n维列向量,b为m维列向量
	//================计算H逆乘A转置=====================
	double** rHtA = new double*[n]; //n*m维
	for(int i = 0; i < n; i++)
	{
		rHtA[i] = new double[m];
	}
	twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
	//=====================================================
	//===============计算A乘H逆乘A转置=====================
	double** ArH = new double*[m]; //m*n维
	for(int i = 0; i < m; i++)
	{
		ArH[i] = new double[n];
	}
	twoMatrixMultiply(A,m,n,"N",reverseH,n,n,"N",ArH);
	double** ArHtA = new double*[m];
	for(int i = 0; i < m; i++)
	{
		ArHtA[i] = new double[m];
	}
	twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
	inverseMatrix(ArHtA, m, m);
	//释放内存空间
	for(int i = 0; i < m; i++)
	{
		delete[] ArH[i];
	}
	delete[] ArH;

	twoMatrixMultiply(rHtA,n,m,"N",ArHtA,m,m,"N",R);

	//释放内存
	for(int i = 0; i < n; i++)
	{
		delete[] rHtA[i];
	}
	delete[] rHtA;
	for(int i = 0; i < m; i++)
	{
		delete[] ArHtA[i];
	}
	delete[] ArHtA;
}
void QuadEquLag::twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out)
{
	if(strcmp(aType,"N") == 0 && strcmp(bType,"N") == 0)//a*b
	{
		for(int i = 0; i < aRow; i++)
		{
			for(int j = 0; j < bCol; j++)
			{
				double tmp = 0;
				for(int k = 0; k < aCol; k++)
				{
					tmp = tmp + a[i][k]*b[k][j];
				}
				out[i][j] = tmp;
			}
		}
	}
	else if(strcmp(aType,"N") == 0 && strcmp(bType,"T") == 0)//a*trans(b)
	{
		for(int i = 0; i < aRow; i++)
		{
			for(int j = 0; j < bRow; j++)
			{
				double tmp = 0.0;
				for(int k = 0; k < bCol; k++)
				{
					tmp = tmp + a[i][k]*b[j][k];
				}
				out[i][j] = tmp;
			}
		}
	}
	else if(strcmp(aType,"T") == 0 && strcmp(bType,"N") == 0)//trans(a)*b
	{

	}
	else if(strcmp(aType,"T") == 0 && strcmp(bType,"T") == 0)//trans(a)*trans(b)
	{

	}

}
void QuadEquLag::inverseMatrix(double** a, int& aRow, int& aCol)
{
	int N = aRow;
	double* A = new double[aRow * aCol];
	
	int LDA = aCol;
	int* IPIV = new int[N];
	int LWORK = N;
	double* WORK = new double[N];
	int INFO ;
	
	int k = 0;
	
	for(int i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			A[k++] = a[j][i];	
		}
	}
	int m = N; int n = m; 
	dgetrf_(&m, &n, A, &m, IPIV, &INFO);
	if(INFO != 0)
	{
		printf("dgetrf is wrong");
		assert(0);
	}
	INFO = 1;
	dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
	if(INFO != 0)
	{
		printf("dgetri is wrong");
		assert(0);
	}
	k = 0;
	for(int i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			a[j][i] = A[k++];	
		}
	}
	delete[] A;
}

//double* QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
//{
//	//H为n维对称阵,A为m*n矩阵,c为n维列向量,b为m维列向量
//
//	
//
//	//vector<int> J0 ;
//	//
//	1)计算有效集J0,构造矩阵Aa0,设k=0
//	//for(int i = 0; i < m; i++)
//	//{
//	//		
//	//	double tmp = 0.0;
//	//	for(int j = 0; j < n; j++)
//	//	{
//	//		tmp = tmp + A[i][j] * x0[j];
//	//	}
//	//	if(tmp == b[i])
//	//	{
//	//		J0.push_back(i);
//	//	}
//	//}
//	//int nn=0;
//	//while(true)
//	//{
//	//	//2)计算gk,并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
//	//	double* gk = new double[n];
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		double tmp = 0.0;
//	//		for(int j = 0; j < n; j++)
//	//		{
//	//			tmp = tmp + H[i][j] * x0[j];
//	//		}
//	//		gk[i] = c[i] + tmp;
//	//	}
//	//	double** Aakgk = new double*[n];
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		Aakgk[i] = new double[J0.size() + 1];
//	//	}
//	//	
//	//	for(unsigned int i = 0; i < J0.size(); i++)
//	//	{
//	//		for(int j = 0; j < n; j++)
//	//		{
//	//			Aakgk[j][i] = A[J0[i]][j];
//	//		}
//	//	}
//	//	
//	//	int last = J0.size();
//	//	for(int j = 0; j < n; j++)
//	//	{
//	//		Aakgk[j][last] = gk[j];
//	//	}
//	//	
//	//	int col = J0.size() + 1;
//	//	int rankAakgk = rank(Aakgk,n,col);
//	//	col--;
//	//	int rankAak = rank(Aakgk, n, col);
//
//	//	if(rankAakgk == rankAak)
//	//	{
//	//	//3)如果gk中的值全大于零,则迭代停止,找到最优值,gk相当于罚参数
//	//		solveLinearEquations(Aakgk, n, col, gk);//gk作为输入参数和输出参数,最终的方程组的解存储在此
//	//	
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			delete[] Aakgk[i];
//	//		}
//	//		delete[] Aakgk;
//	//	
//	//		/*if(greaterThanZero(gk,n) == true)
//	//		{
//	//			delete[] gk;
//	//			break;
//	//		}*/
//	//		if(greaterThanZero(gk,n) == false)
//	//		{
//	//			updateActiveSet(J0, gk );//更新有效集,将gk中的负值对应的有效集去掉,如果gk中全部为负值,去掉最小的负值所对应的有效集
//	//			/*delete[] gk;*/
//	//		}
//	//	}
//	//
//	//	//4)求解dk
//
//	//	double** Aak = new double*[n];
//	//	int mTmp ;
//	//	double* bTmp = NULL;
//	//	if(J0.size() == 0)
//	//	{
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			Aak[i] = new double[1];
//	//			Aak[i][0] = 0.0;
//	//		}
//	//		mTmp = 1;
//	//		bTmp = new double[1];
//	//		for(unsigned int i = 0; i < 1; i++)
//	//		{
//	//			bTmp[i] = 0.0;
//	//		}
//	//	}
//	//	else
//	//	{
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			Aak[i] = new double[J0.size()];
//	//		}
//	//		for(unsigned int i = 0; i < J0.size(); i++)
//	//		{
//	//			for(int j = 0; j < n; j++)
//	//			{
//	//				Aak[j][i] = A[J0[i]][j];
//	//			}
//	//		}
//	//	
//	//		mTmp =(unsigned) J0.size();
//	//		bTmp = new double[J0.size()];
//	//		for(unsigned int i = 0; i < J0.size(); i++)
//	//		{
//	//			bTmp[i] = 0.0;
//	//		}
//	//	}
//	//	cout<<nn++<<endl;
//
//	//	double* dk = QuadEquLag::optimizer(H, mTmp, n, gk, Aak, bTmp);
//	//	if(equalZero(dk,n) == true && greaterThanZero(gk,n) == true)
//	//	{
//	//		delete[] gk;
//	//		break;
//	//	}
//	//	{
//	//		delete[] gk;
//	//	}
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		delete[] Aak[i];
//	//	}
//	//	delete[] Aak;
//	//	delete[] bTmp;
//	//	//5)求解afak,并设置x(k+1) = x(k) + afak*dk;
//	//	double* xk = x0;
//	//	int minIndex;
//	//	double afak = findAfak(J0, dk,A, b, xk,minIndex, m,n);
//	//	updataXk(xk,afak,dk,n);
//
//	//	//6)
//	//	if(afak < 1)
//	//	{
//	//		J0.push_back(minIndex);
//	//	}
//	//	sort(J0.begin(), J0.end());
//	//}
//	
//	return x0;
//}

double QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
{
	//H为n维对称阵,A为m*n矩阵,c为n维列向量,b为m维列向量
	vector<int> J0 ;
	bool come5 = true;
	//0)计算有效集J0,构造矩阵Aa0,设k=0
	for(int i = 0; i < m; i++)
	{
			
		double tmp = 0.0;
		for(int j = 0; j < n; j++)
		{
			tmp = tmp + A[i][j] * x0[j];
		}
		if(tmp == b[i])
		{
			J0.push_back(i);
		}
	}
	int nn=0;
	double* gk = new double[n];
	int nm = 0;
	while(true)
	{
		//1)计算gk,并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
		cout << nm++ <<endl;
		if(come5 == true)
		{
			come5 = false;

			
			for(int i = 0; i < n; i++)
			{
				double tmp = 0.0;
				for(int j = 0; j < n; j++)
				{
					tmp = tmp + H[i][j] * x0[j];
				}
				gk[i] = c[i] + tmp;
			}
		}
		//2)求解线性方程组[Q A'; A 0]*[dk v]' = [-gk 0]
		double* gkk = new double[n+J0.size()];
		double** QA = new double*[(n+J0.size())];
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			QA[i] = new double[n+J0.size()];
		}
		for(int i = 0; i < n; i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[i][j] = H[i][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[i+n][j] = A[J0[i]][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[j][i+n] = A[J0[i]][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(unsigned int j = 0; j < J0.size(); j++)
			{
				QA[i+n][j+n] = 0.0;
			}
		}
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			if(i < n)
				gkk[i] = -gk[i];
			else
				gkk[i] = 0.0;
		}
		int row = n + J0.size();
		int col = n + J0.size();

		double* gkkk = new double[n+J0.size()];
		for(unsigned int i = 0; i < n + J0.size(); i++)
		{
			gkkk[i] = gkk[i];
		}
		solveLinearEquations(QA, row, col, gkkk);
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			delete[] QA[i];
		}
		delete[] QA;

		//3)如果dk=0,转到6)
		if(equalZero(gkkk,n) != true)
		{
		//4)计算步长
			int minIndex;
			double afak =findAfak(J0, gkkk, A ,  b, x0, minIndex,  m, n);
			if(afak < 1)
			{
				J0.push_back(minIndex);
			}
			sort(J0.begin(), J0.end());
		//5)更新xk
			updataXk(x0,afak,gkkk,n);
			come5 = true;
			delete[] gkk;
			delete[] gkkk;
			continue;
		}
		else if(equalZero(gkkk,n) == true)
		{
		//6)
			double* v = new double[J0.size()];
			for(unsigned int i = 0; i < J0.size(); i++)
			{
				v[i] = gkkk[n+i];
			}
			int num = J0.size();
			if(greaterThanZero(v,num) == true)
			{
				break;
			}
			else
			{
				updateActiveSet(J0, v);
			}
			delete[] v;
			
			delete[] gkkk;
			continue;

		}
	}
	delete[] gk;
	
	double obj = 0;
	double* o = new double[n];
	for(int i = 0; i < n; i++)
	{
		double tmp = 0.0;
		for(int j = 0; j < n; j++)
		{
			tmp = tmp + x0[j]*H[i][j];
		}
		o[i] = tmp;
	}
	double cx = 0.0;
	for(int i = 0; i < n; i++)
	{
		obj = obj + o[i] * x0[i] ;
		cx = cx + c[i] * x0[i];
	}
	delete[] o;

	return obj * 0.5 + cx;
		
		;
}

bool QuadProg::equalZero(double* in, int& n)
{
	for(int i = 0; i < n; i++)
	{
		if(in[i] > 0.0000000001)
			return false;
	}
	return true;
}

void QuadProg::updataXk(double* xk,double& afak, double* dk, int& n)
{
	for(int i = 0; i < n; i++)
	{
		xk[i] = xk[i] + afak * dk[i];
	}
}
double QuadProg::findAfak(vector<int>& activeSet, double* dk, double** A, double* b, double* xk, int& minIndex, int& m, int& n)
{
	
	double minVal = 0.0;
	vector<double> minVals;
	vector<int> minIndexs;
	
	for(int i = 0; i < m; i++)
	{
		vector<int>::iterator iter = find(activeSet.begin(),activeSet.end(),i);
		if(iter == activeSet.end()) //i不在有效集中
		{
			double aixk = 0.0;
			double aidk = 0.0;
			for(int j = 0; j < n; j++)
			{
				aidk = aidk + A[i][j] * dk[j];
				aixk = aixk + A[i][j] * xk[j];
			}
			if(aidk > 0.0)
			{
				double rating = (-aixk + b[i])/(aidk);
				minVals.push_back(rating);
				minIndexs.push_back(i);
			}
		}
	}
	minVal = minVals[0];
	minIndex = minIndexs[0];
	for(unsigned int i = 1; i< minVals.size(); i++)
	{
		if(minVals[i] < minVal)
		{
			minVal = minVals[i];
			minIndex = minIndexs[i];
		}
	}
	if(1 < minVal)
		return 1.0;
	else
		return minVal;
}
void QuadProg::updateActiveSet(vector<int>& activeSet, double* mu)
{
	int m =(int) activeSet.size();
	int minIndex = 0;
	double minVal = mu[0];
	for(int i = 1; i < m; i++)
	{
		if(mu[i] < minVal)
		{
			minVal = mu[i];
			minIndex = i;
		}
	}

	if(minVal < 0.0)
	{
		activeSet.erase(activeSet.begin() + minIndex);
	}
	
}
bool QuadProg::greaterThanZero(double* in, int& n)
{
	for(int i = 0; i < n; i++)
	{
		if(in[i] < 0.0)
			return false;
	}
	return true;
}
int QuadProg::rank(double** a,int& row, int& col, double tol)
{
	/*s = svd(A);rank = sum(s > tol);*/
	//===========SVD=============
	char jobu[]="A";
	char jobvt[]="A";
	int m = row;
	int n = col;
	double *b=new double[m*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<m;j++)
		{
			b[i*m+j]=a[j][i];
		}
	}
	int lda=m;

	int minVal = min(m,n);
	double *s=new double[minVal];
	
	int ldu=m;
	double *u=new double[ldu*m];
	
	int ldvt=n;
	double *vt=new double[ldvt*n];
	int lwork=max(1,5*min(m,n));
	double *work=new double[lwork];
	
	int info;

	dgesvd_(jobu, jobvt, &m, &n, b, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,&info);
	//===========================
	delete[] b;
	delete[] u;
	delete[] vt;
	int rank = 0;
	for(int i = 0; i < minVal; i++)
	{
		if(s[i] > tol)
			rank++;
	}

	delete[] s;
	return rank;
}

void QuadProg::solveLinearEquations(double** a, int& row, int& col, double* b)
{
	int N = row;
	int NRHS = 1;
	double* A = new double[(row) * (col)];

	int index = 0;
	for(int i = 0; i < row; i++)
	{
		for(int j = 0; j < col; j++)
		{
			A[index++] = a[j][i];
		}
	}

	int LDA = row;
	int* IPIV = new int[N];
	double* B = b;
	int LDB = max(1, N);
	int INFO;
	dgesv_(&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO );

	delete[] A;
	delete[] IPIV;
}

测试主程序为:

#include "src\Optimizer\QP.h"
#include <fstream>
using namespace std;

int main(int argv, char** argc)
{
	double** H = new double*[400];
	double** A = new double*[400];
	for(int i = 0; i < 400; i++)
	{
		A[i] = new double[400];
		H[i] = new double[400];
	}

	ifstream readin;
	readin.open("./data/H.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			for(int j = 0; j < 400; j++)
			{
				double tmp;
				readin >> tmp;
				H[i][j] = tmp;
			}
		}
	}
	readin.close();

	double* f = new double[400];
	readin.open("./data/f.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			double tmp;
			readin >> tmp;
			f[i] = tmp;
		}
	}
	readin.close();

	for(int i = 0; i < 400; i++)
	{
		for(int j = 0; j < 400; j++)
		{
			A[i][j] = 0.0;
		}
		A[i][i] = 1.0;
	}

	double* b = new double[400];
	readin.open("./data/b.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			double tmp;
			readin >> tmp;
			b[i] = tmp;
		}
	}
	b[399] = 1000000000;
	readin.close();
	double* x0 = new double[400];
	for(int i = 0; i < 400; i++)
	{
		x0[i] = 0.5;
	}
	
	int m = 400;
	int n = 400;
	QuadProg* qu = new QuadProg();
	double obj = qu->optimizer( H,m,n,  f, A, b, x0);
	return 1;
}
本程序中用到了Lapack库和测试数据(Lapack是进行矩阵计算的库,Matlab中也采用了它), 可以到此处下载,漏掉的Lapack.h文件见这里: 点击打开链接
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值