lapackdriver.cpp

图像分割算法

	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of ZGELSS is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

// Gelsd
Matrix<float> Lapack::Gelsd(Matrix<float>& A, Matrix<float>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<float> F = A.Clone();
	Matrix<float> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<float>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int smlsiz = 25;
	int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
	nlvl = (nlvl>0)?nlvl:0;
	int bs = 0;
	if(m>=n)
	{
		bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	else
	{
		bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	Vector<float> work(bs);
	int iwork_size = 3*min*nlvl + 11*min;
	Vector<int> iwork(iwork_size);
	Vector<float> s(min);
	float rcond = max*numeric_limits<float>::epsilon();

	SGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);

	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SGELSD is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<double> Lapack::Gelsd(Matrix<double>& A, Matrix<double>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<double> F = A.Clone();
	Matrix<double> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<double>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int smlsiz = 25;
	int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
	nlvl = (nlvl>0)?nlvl:0;
	int bs = 0;
	if(m>=n)
	{
		bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	else
	{
		bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	Vector<double> work(bs);
	int iwork_size = 3*min*nlvl + 11*min;
	Vector<int> iwork(iwork_size);
	Vector<double> s(min);
	double rcond = max*numeric_limits<double>::epsilon();

	DGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);

	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DGELSD is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexFloat> Lapack::Gelsd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexFloat> F = A.Clone();
	Matrix<ComplexFloat> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int smlsiz = 25;
	int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
	nlvl = (nlvl>0)?nlvl:0;
	int bs = -1;

	Vector<ComplexFloat> work(1);
	int lrwork = 0;
	if(m>=n)
	{
		lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	else
	{
		lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	Vector<float> rwork(lrwork);
	int iwork_size = 3*min*nlvl + 11*min;
	Vector<int> iwork(iwork_size);
	Vector<float> s(min);
	float rcond = max*numeric_limits<float>::epsilon();

	CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
	bs = (int)real(work[0]);
	work = Vector<ComplexFloat>(bs);
	CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);

	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of CGELSD is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}

Matrix<ComplexDouble> Lapack::Gelsd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
	if(A.Rows() != B.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	int info;
	int rank;
	Matrix<ComplexDouble> F = A.Clone();
	Matrix<ComplexDouble> X;
	if(B.Rows() >= A.Columns())
	{
		X = B.Clone();
	}
	else
	{
		X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
		X.ReadFromMatrix(B);
	}
	int m = F.Rows();
	int nrhs = X.Columns();
	int n = F.Columns();
	int ldb = X.Rows();

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int smlsiz = 25;
	int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
	nlvl = (nlvl>0)?nlvl:0;
	int bs = -1;

	Vector<ComplexDouble> work(1);
	int lrwork = 0;
	if(m>=n)
	{
		lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	else
	{
		lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
	}
	Vector<double> rwork(lrwork);
	int iwork_size = 3*min*nlvl + 11*min;
	Vector<int> iwork(iwork_size);
	Vector<double> s(min);
	double rcond = max*numeric_limits<double>::epsilon();
	
	ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
	bs = (int)real(work[0]);
	work = Vector<ComplexDouble>(bs);
	ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);

	if(rank < min && info == 0)
	{
			char s[500];
			sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
			Utility::Warning(s);
	}
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of ZGELSD is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	if(B.Rows() > A.Columns())
	{
		X = X(0, A.Columns()-1, 0, nrhs-1);
	}
	return X;
}






//=========================================================
// Generalized Linear Least Squares problems
//=========================================================


// Gglse
Vector<float> Lapack::Gglse(Matrix<float>& A, Matrix<float>& B, Vector<float>& c, Vector<float>& d)
{
	if(A.Columns() != B.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	if(A.Rows() != c.Length() || B.Rows() != d.Length())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions does not match!");
	}

	if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
	}

	int info;
	Matrix<float> A2 = A.Clone();
	Matrix<float> B2 = B.Clone();
	Vector<float> c2 = c.Clone();
	Vector<float> d2 = d.Clone();

	int m = A2.Rows();
	int n = A2.Columns();
	int p = B2.Rows();
	Vector<float> x(n);

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int nb = 64;
	int bs = p + min + max*nb;
	Vector<float> work(bs);

	SGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);

	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SGGLSE is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "Problem experienced when running SGGLSE!");
			Utility::RunTimeError(s);
		}
	}
	return x;
}



Vector<double> Lapack::Gglse(Matrix<double>& A, Matrix<double>& B, Vector<double>& c, Vector<double>& d)
{
	if(A.Columns() != B.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}
	if(A.Rows() != c.Length() || B.Rows() != d.Length())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions does not match!");
	}

	if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
	}

	int info;
	Matrix<double> A2 = A.Clone();
	Matrix<double> B2 = B.Clone();
	Vector<double> c2 = c.Clone();
	Vector<double> d2 = d.Clone();

	int m = A2.Rows();
	int n = A2.Columns();
	int p = B2.Rows();
	Vector<double> x(n);

	int min = (m>n)?n:m;
	int max = (m>n)?m:n;
	int nb = 64;
	int bs = p + min + max*nb;
	Vector<double> work(bs);

	DGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);

	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of DGGLSE is an illegal value!", -info);
			Utility::RunTimeError(s);
		}
		else if(info > 0)  
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;