Support for function 'schur' added

q

11 files changed, 575 insertions, 4 deletions

diff --git a/src/c/CACSD/includes/syslin.h b/src/c/CACSD/includes/syslin.h
new file mode 100644
index 0000000..7d361c9
--- /dev/null
+++ b/src/c/CACSD/includes/syslin.h
@@ -0,0 +1,31 @@
+/* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+*/
+ 
+#ifndef __SYSLIN_H__
+#define __SYSLIN_H__
+
+#include <stdlib.h>
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+void dsyslina(double* A, int no_of_states, double* B, int no_of_inputs, \
+	double* C, int no_of_outputs, double* D, double* X0, double* out);
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__SYSLIN_H__*/
\ No newline at end of file
diff --git a/src/c/CACSD/interfaces/int_syslin.h b/src/c/CACSD/interfaces/int_syslin.h
new file mode 100644
index 0000000..2327c38
--- /dev/null
+++ b/src/c/CACSD/interfaces/int_syslin.h
@@ -0,0 +1,78 @@
+/* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+*/
+
+#ifndef __INT_SYSLIN_H__
+#define __INT_SYSLIN_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+#define g2d2d2d2syslinss2(in1,size1,in2,size2,in3,size3,in4,size4,out) \
+	dsyslina(in2,size2[0],in3,size3[1],in4,size4[0],NULL,NULL,out)
+
+#define g2d2d2d2d2syslinss2(in1,size1,in2,size2,in3,size3,in4,size4,in5, \
+	size5,out)	dsyslina(in2,size2[0],in3,size3[1],in4,size4[0],in5,NULL,out)
+
+#define g2d2d2d2d2d2syslinss2(in1,size1,in2,size2,in3,size3,in4,size4, \
+ 	in5,size5,in6,size6,out)	dsyslina(in2,size2[0],in3,size3[1],in4, \
+ 	size4[0],in5,in6,out)
+//
+#define g2d0d0d0syslinss2(in1,size1,in2,in3,in4,out) \
+	dsyslina(in2,1,in3,1,in4,1,NULL,NULL,out)
+
+#define g2d0d2d0syslinss2(in1,size1,in2,in3,size3,in4,out) \
+	dsyslina(in2,1,in3,size3[1],in4,1,NULL,NULL,out)
+
+#define g2d0d0d2syslinss2(in1,size1,in2,in3,in4,size4,out) \
+	dsyslina(in2,1,in3,1,in4,size4[0],NULL,NULL,out)
+
+#define g2d0d2d2syslinss2(in1,size1,in2,in3,size3,in4,size4,out) \
+	dsyslina(in2,1,in3,size3[1],in4,size4[0],NULL,NULL,out)
+//
+#define g2d0d0d0d0syslinss2(in1,size1,in2,in3,in4,in5,out) \
+	dsyslina(in2,1,in3,1,in4,1,in5,NULL,out)
+
+#define g2d0d2d0d2syslinss2(in1,size1,in2,in3,size3,in4,in5,size5, \
+	out) dsyslina(in2,1,in3,size3[1],in4,1,in5,NULL,out)
+
+#define g2d0d0d2d2syslinss2(in1,size1,in2,in3,in4,size4,in5,size5, \
+	out) dsyslina(in2,1,in3,1,in4,size4[0],in5,NULL,out)
+
+#define g2d0d2d2d2syslinss2(in1,size1,in2,in3,size3,in4,size4,in5, \
+	size5, out) dsyslina(in2,1,in3,size3[1],in4,size4[0],in5,NULL,out)
+
+#define g2d2d2d2d0syslinss2(in1,size1,in2,size2,in3,size3,in4,size4, \
+	in5,out) dsyslina(in2,size2[0],in3,size3[1],in4,size4[0],in5,NULL,out)
+//
+#define g2d0d0d0d0d0syslinss2(in1,size1,in2,in3,in4,in5,in6,out) \
+	dsyslina(in2,1,in3,1,in4,1,in5,in6,out)
+
+#define g2d0d2d0d2d0syslinss2(in1,size1,in2,in3,size3,in4,in5,size5, \
+	in6,out) dsyslina(in2,1,in3,size3[1],in4,1,in5,in6,out)
+
+#define g2d0d0d2d2d0syslinss2(in1,size1,in2,in3,in4,size4,in5,size5, \
+	in6,out) dsyslina(in2,1,in3,1,in4,size4[0],in5,in6,out)
+
+#define g2d0d2d2d2d0syslinss2(in1,size1,in2,in3,size3,in4,size4,in5, \
+	size5,in6,out) dsyslina(in2,1,in3,size3[1],in4,size4[0],in5,in6,out)
+
+#define g2d2d2d2d0d2syslinss2(in1,size1,in2,size2,in3,size3,in4,size4, \
+	in5,in6,size6,out) dsyslina(in2,size2[0],in3,size3[1],in4,size4[0], \
+	in5,in6,out)
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+
+#endif /*__INT_SYSLIN_H__*/
\ No newline at end of file
diff --git a/src/c/CACSD/syslin/dsyslina.c b/src/c/CACSD/syslin/dsyslina.c
new file mode 100644
index 0000000..cd7fa6b
--- /dev/null
+++ b/src/c/CACSD/syslin/dsyslina.c
@@ -0,0 +1,94 @@
+/* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+*/
+
+/*Function for declaring state space system*/
+
+#include <stdlib.h>
+
+void dsyslina(double* A, int no_of_states, double* B, int no_of_inputs, \
+	double* C, int no_of_outputs, double* D, double* X0, double* out)
+{
+	int row = 0,col = 0;
+    int no_of_cols = no_of_states + no_of_inputs + 1;
+
+    /*Copy matrix A into out matrix*/
+    for(row = 0; row < no_of_states; row++)
+    {
+    	for(col = 0; col < no_of_states; col++)
+    	{
+    		out[row*no_of_cols + col] = A[row*no_of_states + col];
+    	}
+    	
+    }
+
+    /*Copy matrix B in out matrix*/
+    for(row = 0; row < no_of_states; row++)
+    {
+    	for(col=0; col < no_of_inputs; col++)
+    	{
+    		out[row * no_of_cols + no_of_states + col] = B[row * no_of_inputs + col];
+    	}
+    }
+
+    /*Copy matrix C in out matrix*/
+    for(row = 0; row < no_of_outputs; row++)
+    {
+    	for(col=0; col < no_of_states; col++)
+    	{
+    		out[(no_of_states + row)*no_of_cols + col] = C[row * no_of_states + col];
+    	}
+    }
+
+	/*Copy matrix D in out matrix*/
+	if( D != NULL)
+	{
+	    for(row = 0; row < no_of_outputs; row++)
+		{
+	    	for(col=0; col < no_of_inputs; col++)
+	    	{
+	    		out[(no_of_states+row)*no_of_cols + no_of_states+col] = \
+	    			D[row * no_of_inputs + col];
+	    	}
+	    }
+	}
+	else
+	{
+	    for(row = 0; row < no_of_outputs; row++)
+		{
+	    	for(col=0; col < no_of_inputs; col++)
+	    	{
+	    		out[(no_of_states+row)*no_of_cols + no_of_states+col] = 0;
+	    	}
+	    }
+	}	
+
+	/*Copy matrix X0 in out matrix*/
+	if( X0 != NULL)
+	{
+	    for(row = 0; row < no_of_states; row++)
+		{
+	    	out[(row+1)*(no_of_cols) - 1] = X0[row];
+	    	
+	    }
+	}
+	else
+	{
+	    for(row = 0; row < no_of_states; row++)
+		{
+	    	out[(row+1)*(no_of_cols) - 1] = 0;
+	    	
+	    }
+	}	
+
+
+
+}
\ No newline at end of file
diff --git a/src/c/differential_calculus/includes/ode.h b/src/c/differential_calculus/includes/ode.h
index 7996243..7e11003 100644
--- a/src/c/differential_calculus/includes/ode.h
+++ b/src/c/differential_calculus/includes/ode.h
@@ -14,11 +14,11 @@
 #define __ODE_H__
 
 double dodes(double initial_value, double start_time, double end_time, \
-			int (*ode_function), char *solver_type, double nequs, double eps_abs, double eps_rel, \
+			int (*ode_function)(double, double*, double*, int*), char *solver_type, double nequs, double eps_abs, double eps_rel, \
 			double step_size, int *params);
 
 void dodea(double *initial_value, double start_time, double end_time, \
-			int (*ode_function), char *solver_type, double nequs, double eps_abs, double eps_rel, \
+			int (*ode_function)(double, double*, double*, int*), char *solver_type, double nequs, double eps_abs, double eps_rel, \
 			double step_size, int *params, double *out);
 
 #endif /*__ODE_H__*/
\ No newline at end of file
diff --git a/src/c/differential_calculus/interfaces/int_ode.h b/src/c/differential_calculus/interfaces/int_ode.h
index 28f4399..51764ea 100644
--- a/src/c/differential_calculus/interfaces/int_ode.h
+++ b/src/c/differential_calculus/interfaces/int_ode.h
@@ -46,6 +46,20 @@ extern "C" {
 					in2, in3, func_name, solvertype, size1[1], 1.0e-2, 1.0e-2, \
 					1.0e-6, size1, out)
 
+#define d0d0d0d0d0fn0oded0(in1, in2, in3, in4, in5, func_name) \
+					dodes(in1, in2, in3, func_name, "rkf",1, in5, in4, 1.0e-6, NULL)
+
+#define d2d0d0d0d0fn0oded2(in1, size1, in2, in3, in4, in5, func_name, out) \
+					dodea(in1, in2, in3, func_name, "rkf", size1[1], in5, in4, \
+					1.0e-6, size1, out)
+
+#define d0d0d2d0d0fn0oded2(in1, in2, in3, size3, in4, in5, func_name, out) \
+					dodea(in1, in2, in3, func_name, "rkf", 1, in5, in4, 1.0e-6, size3, out)
+
+#define d2d0d2d0d0fn0oded2(in1, size1, in2, in3, size3, in4, in5, func_name, out) \
+					dodea(in1, in2, in3, func_name, "rkf", size1[1], in5, in4, \
+					1.0e-6, size1, out)
+
 #ifdef  __cplusplus
 } /* extern "C" */
 #endif
diff --git a/src/c/differential_calculus/ode/dodea.c b/src/c/differential_calculus/ode/dodea.c
index 1cb07fa..97d56a1 100644
--- a/src/c/differential_calculus/ode/dodea.c
+++ b/src/c/differential_calculus/ode/dodea.c
@@ -20,7 +20,7 @@
 
 
 void dodea(double *initial_value, double start_time, double end_time, \
-			int (*ode_function), char *solver_type, double nequs, double eps_abs, \
+			int (*ode_function)(double, double*, double*, int*), char *solver_type, double nequs, double eps_abs, \
 			double eps_rel, double step_size, int *params, double *out)
 {
 	double t = start_time;
diff --git a/src/c/differential_calculus/ode/dodes.c b/src/c/differential_calculus/ode/dodes.c
index adef1ba..6d1a867 100644
--- a/src/c/differential_calculus/ode/dodes.c
+++ b/src/c/differential_calculus/ode/dodes.c
@@ -20,7 +20,8 @@
 
 
 double dodes(double initial_value, double start_time, double end_time, \
-			int (*ode_function), char *solver_type, double nequs, double eps_abs, \
+			int (*ode_function)(double, double*, double*, int*), \
+			char *solver_type, double nequs, double eps_abs, \
 			double eps_rel, double step_size, int *params)
 {
 	double out = 0, t = 0;
diff --git a/src/c/linearAlgebra/includes/schur.h b/src/c/linearAlgebra/includes/schur.h
new file mode 100644
index 0000000..27b5c2a
--- /dev/null
+++ b/src/c/linearAlgebra/includes/schur.h
@@ -0,0 +1,30 @@
+ /* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __SCHUR_H__
+#define __SCHUR_H__
+
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+double dschura(double* in1, int size, int flag, int nout, double* out1, \
+		double* out2);
+
+double dgschura(double* in1, int size, double* in2, int flag, int nout, \
+			double* out1, double* out2, double* out3, double* out4);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__SCHUR_H__*/
diff --git a/src/c/linearAlgebra/interfaces/int_schur.h b/src/c/linearAlgebra/interfaces/int_schur.h
new file mode 100644
index 0000000..81324e6
--- /dev/null
+++ b/src/c/linearAlgebra/interfaces/int_schur.h
@@ -0,0 +1,60 @@
+ /* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_SCHUR_H__
+#define __INT_SCHUR_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+#define d2schurd2(in1,size1,out1) dschura(in1,size1[0],0,1,out1,NULL)
+ 
+#define d2schurd2d2(in1,size1,out1,out2) dschura(in1,size1[0],0,2,out1,out2)	
+
+#define d2g2schurd2(in1,size1,in2,size2,out1) (in2[0]=='c')? \
+			dschura(in1,size1[0],1,1,out1,NULL): \
+			dschura(in1,size1[0],2,1,out1,NULL)
+
+#define d2g2schurd2d0(in1,size1,in2,size2,out1) (in2[0]=='c')? \
+			dschura(in1,size1[0],1,2,out1,NULL): \
+			dschura(in1,size1[0],2,2,out1,NULL)
+
+#define d2g2schurd2d0d2(in1,size1,in2,size2,out1,out2) (in2[0]=='c')? \
+			dschura(in1,size1[0],1,3,out1,out2): \
+			dschura(in1,size1[0],2,3,out1,out2)
+
+#define d2d2schurd2d2(in1,size1,in2,size2,out1,out2) dgschura(in1,size1[0], \
+			in2,0,2,out1,out2,NULL,NULL)
+
+#define d2d2schurd2d2d2d2(in1,size1,in2,size2,out1,out2,out3,out4) \
+			dgschura(in1,size1[0],in2,0,4,out1,out2,out3,out4)
+
+#define d2d2g2schurd0(in1,size1,in2,size2,in3,size3) dgschura(in1,size1[0], \
+			in2,1,1,NULL,NULL,NULL,NULL)
+
+#define d2d2g2schurd2d0(in1,size1,in2,size2,in3,size3,out1) \
+			dgschura(in1,size1[0],in2,1,2,out1,NULL,NULL,NULL)
+
+#define d2d2g2schurd2d2d0(in1,size1,in2,size2,in3,size3,out1,out2) \
+			dgschura(in1,size1[0],in2,1,3,out1,out2,NULL,NULL)
+
+#define d2d2g2schurd2d2d2d0(in1,size1,in2,size2,in3,size3,out1,out2,out3) \
+			dgschura(in1,size1[0],in2,1,4,out1,out2,out3,NULL)
+
+#define d2d2g2schurd2d2d2d2d0(in1,size1,in2,size2,in3,size3,out1,out2,out3, \
+			out4) dgschura(in1,size1[0],in2,1,5,out1,out2,out3,out4)
+			
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_SCHUR_H__*/
diff --git a/src/c/linearAlgebra/schur/dgschura.c b/src/c/linearAlgebra/schur/dgschura.c
new file mode 100644
index 0000000..56d4845
--- /dev/null
+++ b/src/c/linearAlgebra/schur/dgschura.c
@@ -0,0 +1,152 @@
+/* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+
+/*Fucntion to find generalised schur decomposition of given square matrices */
+#include "schur.h"
+#include "lapack.h"
+#include "stdlib.h"
+#include "string.h"
+
+
+/*flag  --> 0: nothing
+		--> 1: continuous
+		--> 2: discrete
+*/
+
+lapack_logical selctg2( double* in1, double* in2, double* in3);
+
+double dgschura(double* in1, int size, double* in2, int flag, int nout, \
+			double* out1, double* out2, double* out3, double* out4)
+{
+	char JOBVSL = 'N';
+	char JOBVSR = 'N';
+	char SORT = 'N';
+	int SDIM = 0;
+	int LDVSL = size, LDVSR = size;
+    int LWORK = 8*size+16, INFO;
+	double *ALPHAR, *ALPHAI, *BETA, *VSL, *VSR, *WORK;
+	int *BWORK;
+	double ret = 0;
+	double *buf1, *buf2; /*input is copied to buf, since lapack function directly
+	modifies the input variable*/
+
+	/*Used incase of flag > 0*/ 
+	LAPACK_D_SELECT3 selctg = &selctg2;  
+
+	if(nout >= 2){
+		 JOBVSL = 'V';
+		 JOBVSR = 'V';
+		}
+	if(flag > 0) SORT = 'S';
+
+	buf1   = (double*) malloc((double) size*size*sizeof(double));
+	buf2   = (double*) malloc((double) size*size*sizeof(double));
+	ALPHAR = (double*) malloc((double) size*sizeof(double));
+	ALPHAI = (double*) malloc((double) size*sizeof(double));
+	BETA   = (double*) malloc((double) size*sizeof(double));
+	VSL = (double*) malloc((double) LDVSL*size*sizeof(double));
+	VSR = (double*) malloc((double) LDVSR*size*sizeof(double));
+	WORK = (double*) malloc((double) LWORK*sizeof(double));
+	BWORK = (int*) malloc((double) size*sizeof(double));
+
+
+	memcpy(buf1,in1,size*size*sizeof(double));
+	memcpy(buf2,in2,size*size*sizeof(double));
+
+	dgges_(&JOBVSL,&JOBVSR,&SORT,selctg,&size,buf1,&size,buf2,&size,&SDIM, \
+			ALPHAR,ALPHAI,BETA,VSL,&LDVSL,VSR,&LDVSR,WORK,&LWORK,BWORK,&INFO);
+
+	/*if (INFO != 0) 
+	{
+		out1 = NULL;
+		return 0;
+	}*/
+
+	if(nout == 1)
+	{
+		return(SDIM);
+	}
+	else if(nout == 2)
+	{
+		if(flag == 0)
+		{
+			/*copy in1 to out1 and in2 to out2*/
+			memcpy(out1,buf1,size*size*sizeof(double));
+			memcpy(out2,buf2,size*size*sizeof(double));
+		}
+		else
+		{
+			/*copy VSR to out1 and return SDIM */
+			memcpy(out1,VSR,size*size*sizeof(double));
+			ret = SDIM;
+		}
+	}
+	else if(nout == 3) 
+	{
+		/*copy VSL to out1, VSR to out2, return SDIM*/
+		memcpy(out1,VSL,size*size*sizeof(double));
+		memcpy(out2,VSR,size*size*sizeof(double));
+		ret = SDIM;
+	}
+	else if(nout == 4)
+	{
+		if(flag == 0)
+		{
+			/*copy in1 to out1 and in2 to out2*/
+			memcpy(out1,buf1,size*size*sizeof(double));
+			memcpy(out2,buf2,size*size*sizeof(double));
+			/*copy VSL to out3 and VSR to out4*/
+			memcpy(out3,VSL,size*size*sizeof(double));
+			memcpy(out4,VSR,size*size*sizeof(double));
+		}
+		else
+		{
+			/*copy in1 to out1 and in2 to out2*/
+			memcpy(out1,buf1,size*size*sizeof(double));
+			memcpy(out2,buf2,size*size*sizeof(double));
+			/*copy VSR to out3 and return SDIM*/
+			memcpy(out3,VSR,size*size*sizeof(double));
+			ret = SDIM;
+		}
+	}
+	else /*nout = 5*/
+	{
+		/*copy in1 to out1 and in2 to out2*/
+		memcpy(out1,buf1,size*size*sizeof(double));
+		memcpy(out2,buf2,size*size*sizeof(double));
+		/*copy VSL to out3 and VSR to out4*/
+		memcpy(out3,VSL,size*size*sizeof(double));
+		memcpy(out4,VSR,size*size*sizeof(double));
+		/*return SDIM*/
+		ret = SDIM;
+	}
+
+	free(buf1);
+	free(buf2);
+	free(VSL);
+	free(VSR);
+	free(ALPHAR);
+	free(ALPHAI);
+	free(BETA);
+	free(WORK);
+	free(BWORK);
+
+	return ret;
+}
+
+lapack_logical selctg2(double* in1, double* in2,  double* in3)
+{
+	if((sqrt(*in1**in1+*in2**in2)/ *in3) < 1)
+		return 1;
+	else 
+		return 0;
+}
\ No newline at end of file
diff --git a/src/c/linearAlgebra/schur/dschura.c b/src/c/linearAlgebra/schur/dschura.c
new file mode 100644
index 0000000..c31ddca
--- /dev/null
+++ b/src/c/linearAlgebra/schur/dschura.c
@@ -0,0 +1,111 @@
+/* Copyright (C) 2017 - IIT Bombay - FOSSEE
+
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ Author: Siddhesh Wani
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+
+/*Fucntion to find schur decomposition of given square matrix */
+#include "schur.h"
+#include "lapack.h"
+#include "stdlib.h"
+#include "string.h"
+
+
+/*flag  --> 0: nothing
+		--> 1: continuous
+		--> 2: discrete
+*/
+
+lapack_logical selctg1( double* in1, double* in2);
+
+double dschura(double* in1, int size, int flag, int nout, double* out1, \
+		double* out2)
+{
+	char JOBVS = 'N';
+	char SORT = 'N';
+	int SDIM = 0;
+	int LVDS = size;
+    int LWORK = 3*size, INFO;
+	double *WR, *WI, *VS, *WORK;
+	int *BWORK;
+	double ret = 0;
+	double *buf; /*input is copied to buf, since lapack function direclty
+	modifies the input variable*/
+
+	/*Used incase of flag > 0*/ 
+	LAPACK_D_SELECT2 selctg = &selctg1;  
+
+	if(nout >= 2) JOBVS = 'V';
+	if(flag > 0) SORT = 'S';
+
+	buf = (double*) malloc((double) size*size*sizeof(double));
+	WR = (double*) malloc((double) size*sizeof(double));
+	WI = (double*) malloc((double) size*sizeof(double));
+	VS = (double*) malloc((double) LVDS*size*sizeof(double));
+	WORK = (double*) malloc((double) LWORK*sizeof(double));
+	BWORK = (int*) malloc((double) size*sizeof(double));
+
+
+	memcpy(buf,in1,size*size*sizeof(double));
+
+	dgees_(&JOBVS,&SORT,selctg,&size,buf,&size,&SDIM,WR,WI,VS,&LVDS, \
+		WORK,&LWORK,BWORK,&INFO);
+
+	/*if (INFO != 0)
+	{
+		out1 = NULL;
+		return 0;	
+	} */
+
+	if(nout == 1)
+	{
+		/*Copy result in in1 to out1*/
+		memcpy(out1,buf,size*size*sizeof(double));
+	}
+	else if(nout == 2)
+	{
+		if(flag == 0)
+		{
+			/*copy in1 to out2 and VS to out1*/
+			memcpy(out2,buf,size*size*sizeof(double));
+			memcpy(out1,VS,size*size*sizeof(double));
+		}
+		else
+		{
+			/*copy VS to out1 and SDIM to out2*/
+			memcpy(out1,VS,size*size*sizeof(double));
+			ret = SDIM;
+		}
+	}
+	else 
+	{
+		/*copy VS to out1, SDIM to out2, in1 to out3*/
+		memcpy(out1,VS,size*size*sizeof(double));
+		memcpy(out2,buf,size*size*sizeof(double));
+		ret = SDIM;
+	}
+
+	free(buf);
+	free(WI);
+	free(WR);
+	free(VS);
+	free(WORK);
+	free(BWORK);
+
+	return ret;
+}
+
+lapack_logical selctg1(double* in1, double* in2)
+{
+
+	if(sqrt(*in1**in1+*in2**in2) < 1)
+		return 1;
+	else 
+		return 0;
+}
\ No newline at end of file