/*
Copyright (C) 2004 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

#include <iostream>

#include <vector>
#include <string>
//#include <string.h>

extern "C"{
#include <mathlink.h>
}

//#include <fstream>

using namespace std;

//ofstream out("info.txt");

#ifdef LIMITED
static long maxdim = 20;
char *limited = "Print[\"This is a limited version, the maximum system dimension is 20. You can find the full version at http://www.slicot.de/\"]";
#endif

double scale = 1.2;

extern "C"{

int ab09ad_(char *dico, char *job, char *equil, char* ordsel, long *n, long *m, long *p, long *nr, double *A, long *lda, double *B, long *ldb, double *C, long *ldc, double *hsv, double *tol, long *iwork, double *dwork, long *ldwork, long *iwarn, long *info);

int ab09bd_(char *dico, char *job, char *equil, char* ordsel, long *n, long *m, long *p, long *nr, double *A, long *lda, double *B, long *ldb, double *C, long *ldc, double *D, long *ldd, double *hsv, double *tol1, double *tol2, long *iwork, double *dwork, long *ldwork, long *iwarn, long *info);

int ab09cd_(char *dico, char *equil, char* ordsel, long *n, long *m, long *p, long *nr, double *A, long *lda, double *B, long *ldb, double *C, long *ldc, double *D, long *ldd, double *hsv, double *tol1, double *tol2, long *iwork, double *dwork, long *ldwork, long *iwarn, long *info);

long max3(long a, long b, long c)
{
	return max(max(a, b), c);
}

void runBTA(char *job, char *equil, char *ordsel, long n, long m, long p, long nr, double *A, long nA, double *B, long nB, double *C, long nC, double tol, long ldwork)
{	
	vector<double> hsv(n);
	vector<long> iwork;
	long dim;
	if (job[0] == 'B')
		dim = 1;
	else if (job[0] == 'N')
		dim = n;
	else
	{
		MLPutString(stdlink, "Parameter job is incorrect");
		return;
	}
	iwork.resize(dim);
	dim = n*(2*n+max3(n, m, p)+5)+n*(n+1)/2;
	if (ldwork < dim)
		dim = dim*scale;
	else
		dim = ldwork;
	vector<double> dwork(dim);
	long iwarn = 0;
	long info;
	
#ifdef LIMITED
if (n > maxdim)
{
	MLEvaluateString(stdlink, limited);
	info = 100;
}
else
#endif
	ab09ad_("C", job, equil, ordsel, &n, &m, &p, &nr, A, &n, B, &n, C, &p, &*hsv.begin(), &tol, &*iwork.begin(), &*dwork.begin(), &dim, &iwarn, &info);
	if (info ==0)
	{
		MLPutFunction(stdlink, "List", 6);
		MLPutInteger(stdlink, nr);
		vector<double> tmp(max3(nr, m, p)*nr);
		int k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = A[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*nr);
		k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < m; ++j)
				tmp[k++] = B[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*m);
		k = 0;
		for (int i = 0; i < p; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = C[i + j*p];
		MLPutRealList(stdlink, &*tmp.begin(), nr*p);
		MLPutRealList(stdlink, &*hsv.begin(), n);
		MLPutFunction(stdlink, "List", 2);
		MLPutInteger(stdlink, dim);	
		MLPutInteger(stdlink, int(dwork[0]));	
	}
	else
	{
		MLPutFunction(stdlink, "List", 2);
		MLPutString(stdlink, "Function has failed");
		MLPutInteger(stdlink, info);	
	}	
}

void runSPA(char *job, char *equil, char *ordsel, long n, long m, long p, long nr, double *A, long nA, double *B, long nB, double *C, long nC, double *D, long nD, double tol1, double tol2, long ldwork)
{	
	vector<double> hsv(n);
	vector<long> iwork(2*n);
	long dim = n*(2*n+max3(n, m, p)+5)+n*(n+1)/2;
	if (ldwork < dim)
		dim = dim*scale;
	else
		dim = ldwork;
	vector<double> dwork(dim);
	long iwarn = 0;
	long info;
#ifdef LIMITED
if (n > maxdim)
{
	MLEvaluateString(stdlink, limited);
	info = 100;
}
else
#endif
	ab09bd_("C", job, equil, ordsel, &n, &m, &p, &nr, A, &n, B, &n, C, &p, D, &p, &*hsv.begin(), &tol1, &tol2, &*iwork.begin(), &*dwork.begin(), &dim, &iwarn, &info);
	if (info ==0)
	{
		MLPutFunction(stdlink, "List", 7);
		MLPutInteger(stdlink, nr);
		vector<double> tmp(max3(nr, m, p)*nr);
		int k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = A[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*nr);
		k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < m; ++j)
				tmp[k++] = B[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*m);
		k = 0;
		for (int i = 0; i < p; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = C[i + j*p];
		MLPutRealList(stdlink, &*tmp.begin(), nr*p);
		k = 0;
		for (int i = 0; i < p; ++i)
			for (int j = 0; j < m; ++j)
				tmp[k++] = D[i + j*p];
		MLPutRealList(stdlink, &*tmp.begin(), p*m);
		MLPutRealList(stdlink, &*hsv.begin(), n);
		MLPutFunction(stdlink, "List", 3);
		MLPutInteger(stdlink, dim);	
		MLPutInteger(stdlink, int(dwork[0]));	
		MLPutInteger(stdlink, int(iwork[0]));	
	}
	else
	{
		MLPutFunction(stdlink, "List", 2);
		MLPutString(stdlink, "Function has failed");
		MLPutInteger(stdlink, info);	
	}	
}

void runHNA(char *equil, char *ordsel, long n, long m, long p, long nr, double *A, long nA, double *B, long nB, double *C, long nC, double *D, long nD, double tol1, double tol2, long ldwork)
{	
	vector<double> hsv(n);
	vector<long> iwork(m);
	long dim = n*(2*n+max3(n, m, p)+5)+n*(n+1)/2;
	long dim2 = n*(m+p+2) + 2*m*p + min(m, n) + max(3*m,min(n,m)+p);
	if (dim < dim2) dim = dim2;	
	if (ldwork < dim)
		dim = dim*scale;
	else
		dim = ldwork;
	vector<double> dwork(dim);
	long iwarn = 0;
	long info;
#ifdef LIMITED
if (n > maxdim)
{
	MLEvaluateString(stdlink, limited);
	info = 100;
}
else
#endif
	ab09cd_("C", equil, ordsel, &n, &m, &p, &nr, A, &n, B, &n, C, &p, D, &p, &*hsv.begin(), &tol1, &tol2, &*iwork.begin(), &*dwork.begin(), &dim, &iwarn, &info);
	if (info ==0)
	{
		MLPutFunction(stdlink, "List", 7);
		MLPutInteger(stdlink, nr);
		vector<double> tmp(max3(nr, m, p)*nr);
		int k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = A[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*nr);
		k = 0;
		for (int i = 0; i < nr; ++i)
			for (int j = 0; j < m; ++j)
				tmp[k++] = B[i + j*n];
		MLPutRealList(stdlink, &*tmp.begin(), nr*m);
		k = 0;
		for (int i = 0; i < p; ++i)
			for (int j = 0; j < nr; ++j)
				tmp[k++] = C[i + j*p];
		MLPutRealList(stdlink, &*tmp.begin(), nr*p);
		k = 0;
		for (int i = 0; i < p; ++i)
			for (int j = 0; j < m; ++j)
				tmp[k++] = D[i + j*p];
		MLPutRealList(stdlink, &*tmp.begin(), p*m);
		MLPutRealList(stdlink, &*hsv.begin(), n);
		MLPutFunction(stdlink, "List", 3);
		MLPutInteger(stdlink, dim);	
		MLPutInteger(stdlink, int(dwork[0]));	
		MLPutInteger(stdlink, int(iwork[0]));	
	}
	else
	{
		MLPutFunction(stdlink, "List", 2);
		MLPutString(stdlink, "Function has failed");
		MLPutInteger(stdlink, info);	
	}	
}

}

int main(int argc, char *argv[]) 
{
  return MLMain(argc, argv);
}