Non-linear least squares added

Add wrapper for (just unbounded so far) non-linear least squares functions.
13 years ago · 8d9546f568
2 changed files with 228 additions and 0 deletions
--- a/src/NativeWrappers/MKL/optimization.cpp
+++ b/src/NativeWrappers/MKL/optimization.cpp
@ -0,0 +1,227 @@
+#include <algorithm>
+#include <complex>
+#define MKL_Complex8 std::complex<float>
+#define MKL_Complex16 std::complex<double>
+
+#include "mkl_lapack.h"
+#include "mkl_cblas.h"
+#include "lapack_common.h"
+#include "wrapper_common.h"
+#include "mkl_lapacke.h"
+#include "mkl.h"
+#include "mkl_trans.h"
+
+template<typename T>
+inline MKL_INT lu_factor_2(MKL_INT m, T a[], MKL_INT ipiv[],
+                         void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*))
+{
+    std::complex<double> x = 5;
+	MKL_INT info = 0;
+    getrf(&m, &m, a, &m, ipiv, &info);
+    shift_ipiv_down(m, ipiv);
+    return info;
+};
+
+extern "C" {
+
+//_Mkl_Api(MKL_INT,dtrnlsp_init,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, MKL_INT*, MKL_INT*, double*))
+//_Mkl_Api(MKL_INT,dtrnlsp_check,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlsp_solve,(_TRNSP_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlsp_get,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*))
+//_Mkl_Api(MKL_INT,dtrnlsp_delete,(_TRNSP_HANDLE_t*))
+//
+//_Mkl_Api(MKL_INT,dtrnlspbc_init,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, double*, MKL_INT*, MKL_INT*, double*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_check,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_solve,(_TRNSPBC_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_get,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_delete,(_TRNSPBC_HANDLE_t*))
+//
+//_Mkl_Api(MKL_INT,djacobi_init,(_JACOBIMATRIX_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*))
+//_Mkl_Api(MKL_INT,djacobi_solve,(_JACOBIMATRIX_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,djacobi_delete,(_JACOBIMATRIX_HANDLE_t*))
+//_Mkl_Api(MKL_INT,djacobi,(USRFCND fcn, MKL_INT*, MKL_INT*, double*, double*, double*))
+//_Mkl_Api(MKL_INT,djacobix,(USRFCNXD fcn, MKL_INT*, MKL_INT*, double*, double*, double*,void*))
+
+	typedef double (__stdcall *function)();
+
+	DLLEXPORT double Test(function func, double input[])
+	{
+		input[0] = 5;
+		return func();
+	}
+
+    DLLEXPORT MKL_INT unbound_nonlinearleastsq_init(_TRNSPBC_HANDLE_t* handle, MKL_INT n, MKL_INT m, double x[], double eps[], MKL_INT iter1, MKL_INT iter2, double rs)
+    {
+		return dtrnlsp_init(handle, &n, &m, x, eps, &iter1, &iter2, &rs);
+	}
+
+    DLLEXPORT MKL_INT unbound_nonlinearleastsq_check(_TRNSPBC_HANDLE_t* handle, MKL_INT n, MKL_INT m, double fjac[], double fvec[], double eps[], int info[])
+    {
+		return dtrnlsp_check(handle, &n, &m, fjac, fvec, eps, info);
+	}
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_solve(_TRNSPBC_HANDLE_t* handle, double fvec[], double fjac[], int* RCI_Request)
+    {
+		return dtrnlsp_solve(handle, fvec, fjac, RCI_Request);
+    }
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_get(_TRNSPBC_HANDLE_t* handle, MKL_INT* iter, MKL_INT* st_cr, double* r1, double* r2)
+    {
+		return dtrnlsp_get(handle, iter, st_cr, r1, r2);
+    }
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_delete(_TRNSPBC_HANDLE_t* handle)
+    {
+		return 	dtrnlsp_delete(handle);
+    }
+
+	DLLEXPORT MKL_INT jacobi_init(_JACOBIMATRIX_HANDLE_t* handle, MKL_INT n, MKL_INT m, double x[], double fjac[], double eps)
+	{
+		return djacobi_init(handle, &n, &m, x, fjac, &eps);
+	}
+
+	DLLEXPORT MKL_INT jacobi_solve(_JACOBIMATRIX_HANDLE_t* handle, double f1[], double f2[], MKL_INT* RCI_Request)
+	{
+		return djacobi_solve(handle, f1, f2, RCI_Request);
+	}
+
+	DLLEXPORT MKL_INT jacobi_delete(_TRNSPBC_HANDLE_t* handle)
+    {
+		return djacobi_delete(handle);
+    }
+
+	DLLEXPORT void FreeBuffers()
+	{
+		MKL_FreeBuffers();
+	}
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq(double parameters[], double parametersInitialGuess[], MKL_INT parametersLength, 
+		double residuals[], MKL_INT residualsLength, 
+		double jacobian[], // size parametersLength * residualsLength
+		double residualsMinus[], function updateResidualsMinus,
+		double residualsPlus[], function updateResidualsPlus,
+		function updateResiduals)
+	{
+		double eps[6]; // stop criteria
+		MKL_INT i;
+		for (i = 0; i < 6; i++)
+			eps[i] = 0.00001;
+
+		MKL_INT RCI_Request;     
+		MKL_INT successful;
+
+		MKL_INT maxIterations = 1000, maxTrialStepIterations = 100;
+
+		double rs = 0.0;
+		_TRNSP_HANDLE_t solverHandle;  
+		_JACOBIMATRIX_HANDLE_t jacobianHandle;
+
+		MKL_INT info[6]; // for parameter checking
+
+		double initialStepBound = 0.0;
+
+		double jacobianPrecision = 0.000001;
+
+		MKL_INT iteration;
+
+		double initialResidual, finalResidual;
+
+		// zero initial values:
+		for (i = 0; i < residualsLength; i++)
+			residuals[i] = 0.0;
+		for (i = 0; i < residualsLength * parametersLength; i++)
+			jacobian[i] = 0.0;
+
+		if (dtrnlsp_init(&solverHandle, &parametersLength, &residualsLength, parameters, eps, &maxIterations, &maxTrialStepIterations, &initialStepBound) !=
+			TR_SUCCESS)
+		{                             
+			MKL_FreeBuffers();
+			return 1;
+		}
+	 
+		if (dtrnlsp_check(&solverHandle, &parametersLength, &residualsLength, jacobian, residuals, eps, info) != TR_SUCCESS)
+		{                               
+			MKL_FreeBuffers();
+			return 1;
+		}
+		else
+		{
+			if (info[0] != 0 || // Handle invalid
+				info[1] != 0 || // Jacobian array not valid
+				info[2] != 0 || // Parameters array not valid
+				info[3] != 0)   // Eps array not valid
+			{                                  
+				MKL_FreeBuffers();
+				return 1;
+			}
+		}
+
+		if (djacobi_init(&jacobianHandle, &parametersLength, &residualsLength, parameters, jacobian, &jacobianPrecision) != TR_SUCCESS)
+		{
+			MKL_FreeBuffers ();
+			return 1;
+		}
+
+		RCI_Request = 0;
+		successful = 0;
+		while (successful == 0)
+		{
+			if (dtrnlsp_solve(&solverHandle, residuals, jacobian, &RCI_Request) != TR_SUCCESS)
+			{                             
+				MKL_FreeBuffers();
+				return 1;
+			}
+			if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 ||
+				RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
+				successful = 1;
+			if (RCI_Request == 1) // recalculate function to update parameters
+			{
+				updateResiduals();
+			}
+			if (RCI_Request == 2)
+			{
+				MKL_INT rci_request = 0;
+				MKL_INT jacobianSuccessful = 0;
+				
+				// update Jacobian matrix:
+				while (jacobianSuccessful == 0)
+				{
+					if (djacobi_solve (&jacobianHandle, residualsPlus, residualsMinus, &rci_request) != TR_SUCCESS)
+					{
+						MKL_FreeBuffers ();
+						return 1;
+					}
+					if (rci_request == 1)
+						updateResidualsPlus();
+					else if (rci_request == 2)
+						updateResidualsMinus();
+					else if (rci_request == 0)
+						jacobianSuccessful = 1;
+				}
+			}
+		}
+
+		MKL_INT stopCriterionNumber, iterations;
+		if (dtrnlsp_get (&solverHandle, &iterations, &stopCriterionNumber, &initialResidual, &finalResidual) != TR_SUCCESS)
+		{                                  
+			MKL_FreeBuffers();
+			return 1;
+		}
+
+		if (dtrnlsp_delete (&solverHandle) != TR_SUCCESS)
+		{                                  
+			MKL_FreeBuffers();
+			return 1;
+		}
+
+		if (djacobi_delete (&jacobianHandle) != TR_SUCCESS)
+		{
+			MKL_FreeBuffers ();
+			return 1;
+		}
+                          
+		MKL_FreeBuffers ();
+		return 0;
+	}
+  
+}
--- a/src/NativeWrappers/Windows/MKL/MKLWrapper.vcxproj
+++ b/src/NativeWrappers/Windows/MKL/MKLWrapper.vcxproj
@ -293,6 +293,7 @@
    <ClCompile Include="..\..\MKL\blas.c" />
    <ClCompile Include="..\..\MKL\lapack.cpp" />
    <ClCompile Include="..\..\MKL\vector_functions.c" />
+	<ClCompile Include="..\..\MKL\optimisation.cpp" />
  </ItemGroup>
  <ItemGroup>
    <ResourceCompile Include="..\..\Common\resource.rc" />