lapackdriver.cpp

图像分割算法

			char s[500];
			sprintf(s, "Problem experienced when running DGGLSE!");
			Utility::RunTimeError(s);
		}
	}
	return x;
}


Vector<ComplexFloat> Lapack::Gglse(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& c, Vector<ComplexFloat>& 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<ComplexFloat> A2 = A.Clone();
	Matrix<ComplexFloat> B2 = B.Clone();
	Vector<ComplexFloat> c2 = c.Clone();
	Vector<ComplexFloat> d2 = d.Clone();

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

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

	CGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &m, reinterpret_cast<MKL_Complex8*>(B2.Data()), &p, reinterpret_cast<MKL_Complex8*>(c2.Data()), reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(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 CGGLSE 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 CGGLSE!");
			Utility::RunTimeError(s);
		}
	}
	return x;
}


Vector<ComplexDouble> Lapack::Gglse(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& c, Vector<ComplexDouble>& 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<ComplexDouble> A2 = A.Clone();
	Matrix<ComplexDouble> B2 = B.Clone();
	Vector<ComplexDouble> c2 = c.Clone();
	Vector<ComplexDouble> d2 = d.Clone();

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

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

	ZGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &m, reinterpret_cast<MKL_Complex16*>(B2.Data()), &p, reinterpret_cast<MKL_Complex16*>(c2.Data()), reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(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 ZGGLSE 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 ZGGLSE!");
			Utility::RunTimeError(s);
		}
	}
	return x;
}





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

	if(d.Length() != A.Rows())
	{
		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> d2 = d.Clone();

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

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

	SGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.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 SGGGLM 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 SGGGLM!");
			Utility::RunTimeError(s);
		}
	}

}


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

	if(d.Length() != A.Rows())
	{
		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> d2 = d.Clone();

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

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

	DGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.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 DGGGLM 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 DGGGLM!");
			Utility::RunTimeError(s);
		}
	}
}


void Lapack::Ggglm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& d, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
{
	if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	if(d.Length() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
	}

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

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

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

	CGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(y.Data()), reinterpret_cast<MKL_Complex8*>(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 CGGGLM 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 CGGGLM!");
			Utility::RunTimeError(s);
		}
	}
}


void Lapack::Ggglm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& d, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
{
	if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix dimensions does not match!");
	}

	if(d.Length() != A.Rows())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
	}

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

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

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

	ZGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(y.Data()), reinterpret_cast<MKL_Complex16*>(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 ZGGGLM 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 ZGGGLM!");
			Utility::RunTimeError(s);
		}
	}
}






//=========================================================
// Symmetric eigenproblem
//=========================================================



// Syevd or Heevd
Vector<float> Lapack::Syevd(Matrix<float>& A)
{
	if(A.Rows() != A.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	
	int info;
	Matrix<float> F = A.Clone();
	int n = F.Rows();
	Vector<float> w(n);

	int bs = 2*n+1;
	Vector<float> work(bs);

	int iwork_size = 1;
	Vector<int> iwork(iwork_size);

	SSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);

	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SSYEVD 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 failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	return w;
}


Vector<float> Lapack::Syevd(Matrix<float>& A, Matrix<float>& EV)
{
	if(A.Rows() != A.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	
	int info;
	EV = A.Clone();
	int n = EV.Rows();
	Vector<float> w(n);

	int k = 1;
	if(n>1)
	{
		k = (int)(log((double)(n-1))/log(2.0))+1;
	}
	int bs = 3*n*n+(5+2*k)*n+1;
	Vector<float> work(bs);

	int iwork_size = 5*n+3;
	Vector<int> iwork(iwork_size);

	SSYEVD("V", "U", &n, EV.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
	
	if(info != 0)
	{
		if(info < 0)
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			char s[500];
			sprintf(s, "%i'th argument of SSYEVD 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 failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
			Utility::RunTimeError(s);
		}
	}
	return w;
}


Vector<double> Lapack::Syevd(Matrix<double>& A)
{
	if(A.Rows() != A.Columns())
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	
	int info;
	Matrix<double> F = A.Clone();
	int n = F.Rows();
	Vector<double> w(n);

	int bs = 2*n+1;
	Vector<double> work(bs);

	int iwork_size = 1;
	Vector<int> iwork(iwork_size);

	DSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);

	if(info != 0)