Signed-off-by: Yimin <[email protected]>

GaZ3ll3 · Jan 7, 2015 · f267172 · f267172
1 parent 390024f
commit f267172
Show file tree

Hide file tree

Showing 2 changed files with 86 additions and 2 deletions.
diff --git a/src/Solver/private/Solver.cc b/src/Solver/private/Solver.cc
@@ -19,6 +19,76 @@ Solver::~Solver() {
 #endif
 }
 
+void Solver::Reference(MatlabPtr& DX, MatlabPtr& DY, MatlabPtr Points){
+
+	auto _numberofpoints = mxGetN(Points);
+
+	auto Vander = mxCreateNumericMatrix(_numberofpoints,_numberofpoints,mxDOUBLE_CLASS, mxREAL);
+	auto VanderX = mxCreateNumericMatrix(_numberofpoints, _numberofpoints, mxDOUBLE_CLASS, mxREAL);
+	auto VanderY = mxCreateNumericMatrix(_numberofpoints, _numberofpoints, mxDOUBLE_CLASS, mxREAL);
+
+
+	auto Vander_ptr  = mxGetPr(Vander);
+	auto VanderX_ptr = mxGetPr(VanderX);
+	auto VanderY_ptr = mxGetPr(VanderY);
+
+	auto nodes_ptr   = mxGetPr(Points);
+
+	int deg = round((sqrt(8*_numberofpoints + 1) - 3)/2);
+
+	if ((deg + 1) * (deg + 2)/2 != _numberofpoints){
+		mexErrMsgTxt("Invalid length of input nodes\n");
+	}
+
+	for (size_t col = 0; col < _numberofpoints; col++){
+		for (size_t i = 0; i < deg + 1; i++){
+			for (size_t j = 0; j < i + 1; j++){
+				*Vander_ptr++ = pow(nodes_ptr[2*col], i - j)*pow(nodes_ptr[2*col + 1], j);
+			}
+		}
+	}
+
+	for (size_t col = 0; col < _numberofpoints; col++){
+		for (size_t i = 0; i < deg + 1; i++){
+			for (size_t j = 0; j < i + 1; j++) {
+				if (j == i){
+					*VanderX_ptr++ = 0.;
+				}
+				else{
+					*VanderX_ptr++ = (i - j) * pow(nodes_ptr[2*col], i - j - 1)*pow(nodes_ptr[2*col + 1], j);
+				}
+			}
+		}
+	}
+
+
+
+	for (size_t col = 0; col < _numberofpoints; col++){
+		for (size_t i = 0; i < deg + 1; i++){
+			for (size_t j = 0; j < i + 1; j++) {
+				if (j == 0){
+					*VanderY_ptr++ = 0.;
+				}
+				else{
+					*VanderY_ptr++ = j* pow(nodes_ptr[2*col], i - j)*pow(nodes_ptr[2*col + 1], j - 1);
+				}
+			}
+		}
+	}
+
+
+	mxArray* RHS_x[] = {Vander, VanderX};
+	mexCallMATLAB(1, &DX, 2, RHS_x, "mldivide");
+
+	mxArray* RHS_y[] = {Vander, VanderY};
+	mexCallMATLAB(1, &DY, 2, RHS_y, "mldivide");
+
+	mxDestroyArray(Vander);
+	mxDestroyArray(VanderX);
+	mxDestroyArray(VanderY);
+
+}
+
 } /* namespace MEX */
 
 using namespace MEX;
@@ -40,6 +110,20 @@ MEX_DEFINE(delete) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
 	Session<Solver>::destroy(input.get(0));
 }
 
+MEX_DEFINE(reference) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
+	InputArguments input(nrhs,  prhs, 2);
+	OutputArguments output(nlhs, plhs, 2);
+
+	Solver* solver = Session<Solver>::get(input.get(0));
+
+	size_t numberofpoints = mxGetN(prhs[1]);
+
+	plhs[0] = mxCreateNumericMatrix(numberofpoints, numberofpoints, mxDOUBLE_CLASS, mxREAL);
+	plhs[1] = mxCreateNumericMatrix(numberofpoints, numberofpoints, mxDOUBLE_CLASS, mxREAL);
+
+	solver->Reference(plhs[0], plhs[1], CAST(prhs[1]));
+}
+
 }
 
 MEX_DISPATCH
diff --git a/src/Solver/private/Solver.h b/src/Solver/private/Solver.h
@@ -34,8 +34,8 @@ class Solver {
 	/*
 	 * public methods
 	 */
-
-
+	void Reference(MatlabPtr &DX, MatlabPtr &DY,
+			MatlabPtr Points);
 };
 
 } /* namespace MEX */