From bbd59785a0607149aca4b809d6e43a4b64e2d797 Mon Sep 17 00:00:00 2001
From: Sandeep Gupta
Date: Tue, 20 Jun 2017 17:18:04 +0530
Subject: 20-06-17

---
 .../ToolInitialization/INIT_FillSCI2LibCDirs.bin   | Bin 1628308 -> 1640092 bytes
 .../ToolInitialization/INIT_FillSCI2LibCDirs.sci   |  49 +++++++
 2.3-1/macros/findDeps/getAllSources.bin            | Bin 248996 -> 249160 bytes
 2.3-1/macros/findDeps/getAllSources.sci            |   1 +
 2.3-1/src/c/linearAlgebra/includes/spec.h          |   3 +
 2.3-1/src/c/linearAlgebra/interfaces/int_spec.h    |   7 +
 2.3-1/src/c/linearAlgebra/spec/dspec1a.c           | 151 +++++++++++++++++++--
 2.3-1/src/c/string/disp/zdispa.c                   |   2 +-
 8 files changed, 204 insertions(+), 9 deletions(-)

diff --git a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin
index 55269e36..1c1a52c5 100644
Binary files a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin and b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin differ
diff --git a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci
index 008d7b3c..c44ef3b9 100644
--- a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci
+++ b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci
@@ -4264,6 +4264,49 @@ PrintStringInfo('OUT(2).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y'
 PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
 PrintStringInfo('OUT(2).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
 
+// Edited by Sandeep Gupta, IITB FOSSEE Project.
+PrintStringInfo('NIN=          2',ClassFileName,'file','y');
+PrintStringInfo('NOUT=         1',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+
+PrintStringInfo('NIN=          2',ClassFileName,'file','y');
+PrintStringInfo('NOUT=         2',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP=    IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= ''1''',ClassFileName,'file','y');
+
+PrintStringInfo('NIN=          2',ClassFileName,'file','y');
+PrintStringInfo('NOUT=         3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP=    IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+
+PrintStringInfo('NIN=          2',ClassFileName,'file','y');
+PrintStringInfo('NOUT=         4',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP=    IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).TP=    FA_TP_COMPLEX(IN(1).TP)',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+
 // --- Function List Class. ---
 //NUT: no mixed data types
 ClassFileName = fullfile(SCI2CLibCFLClsDir,ClassName+ExtensionCFuncListCls);
@@ -4292,6 +4335,12 @@ PrintStringInfo('d2'+ArgSeparator+'z2z2',ClassFileName,'file','y');
 PrintStringInfo('c2'+ArgSeparator+'c2c2',ClassFileName,'file','y');
 PrintStringInfo('z2'+ArgSeparator+'z2z2',ClassFileName,'file','y');
 
+PrintStringInfo('d2d2'+ArgSeparator+'z2',ClassFileName,'file','y');
+PrintStringInfo('d2d2'+ArgSeparator+'z2d2',ClassFileName,'file','y');
+PrintStringInfo('d2d2'+ArgSeparator+'z2d2z2',ClassFileName,'file','y');
+PrintStringInfo('d2d2'+ArgSeparator+'z2d2z2z2',ClassFileName,'file','y');
+
+
 // --- Annotation Function And Function List Function. ---
 FunctionName = 'spec'; // AS : Done  AS : Float_Done
 PrintStringInfo('      Adding Function: '+FunctionName+'.',GeneralReport,'file','y');
diff --git a/2.3-1/macros/findDeps/getAllSources.bin b/2.3-1/macros/findDeps/getAllSources.bin
index 401fb991..7bc8fd0d 100644
Binary files a/2.3-1/macros/findDeps/getAllSources.bin and b/2.3-1/macros/findDeps/getAllSources.bin differ
diff --git a/2.3-1/macros/findDeps/getAllSources.sci b/2.3-1/macros/findDeps/getAllSources.sci
index 5dced783..f0f883be 100644
--- a/2.3-1/macros/findDeps/getAllSources.sci
+++ b/2.3-1/macros/findDeps/getAllSources.sci
@@ -1194,6 +1194,7 @@ function allSources = getAllSources(SharedInfo)
 	  "src/c/linearAlgebra/svd/zsvda.c"
 	  "src/c/linearAlgebra/hess/dhessa.c"
 	  "src/c/linearAlgebra/sva/dsvaa.c"
+	  "src/c/linearAlgebra/spec/dspec1a.c"
       "src/c/linearAlgebra/rcond/drconda.c"];
 
   //Files to be inserted only if output format selected is 'Arduino'.
diff --git a/2.3-1/src/c/linearAlgebra/includes/spec.h b/2.3-1/src/c/linearAlgebra/includes/spec.h
index b7b7d35c..3cd0b186 100644
--- a/2.3-1/src/c/linearAlgebra/includes/spec.h
+++ b/2.3-1/src/c/linearAlgebra/includes/spec.h
@@ -49,6 +49,9 @@ EXTERN_MATOPS void dspec2a(double* in, int rows, double* eigenvalues,double* eig
 EXTERN_MATOPS void cspec2a(floatComplex* in, int rows, floatComplex* eigenvalues,floatComplex* eigenvectors);
 EXTERN_MATOPS void zspec2a(doubleComplex* in, int rows,doubleComplex* eigenvalues,doubleComplex* eigenvectors);
 
+/*Edited by - Sandeep Gupta, IITB FOSSEE*/
+void dspec1a(double *in1,double *in2,int size,int nout,doubleComplex *out1,double *out2,doubleComplex *out3,doubleComplex *out4);
+
 #ifdef  __cplusplus
 } /* extern "C" */
 #endif
diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_spec.h b/2.3-1/src/c/linearAlgebra/interfaces/int_spec.h
index c9463292..8772a6a6 100644
--- a/2.3-1/src/c/linearAlgebra/interfaces/int_spec.h
+++ b/2.3-1/src/c/linearAlgebra/interfaces/int_spec.h
@@ -69,4 +69,11 @@
 
 #define z2specz2z2(in,size,out1,out2)		zspec2a(in, size[0], out2, out1)
 
+/*Edited by Sandeep Gupta, IITB, FOSSEE Project.spec(A,B) */
+
+#define d2d2specz2(in1,size,in2,size2,out1) dspec1a(in1,in2,size[0],1,out1,NULL,NULL,NULL)
+#define d2d2specz2d2(in1,size,in2,size2,out1,out2) dspec1a(in1,in2,size[0],2,out1,out2,NULL,NULL)
+#define d2d2specz2d2z2(in1,size,in2,size2,out1,out2,out3) dspec1a(in1,in2,size[0],3,out1,out2,out3,NULL)
+#define d2d2specz2d2z2z2(in1,size,in2,size2,out1,out2,out3,out4) dspec1a(in1,in2,size[0],4,out1,out2,out3,out4)
+
 #endif /* !__INT_SPEC_H__ */
diff --git a/2.3-1/src/c/linearAlgebra/spec/dspec1a.c b/2.3-1/src/c/linearAlgebra/spec/dspec1a.c
index 069454da..28440be6 100644
--- a/2.3-1/src/c/linearAlgebra/spec/dspec1a.c
+++ b/2.3-1/src/c/linearAlgebra/spec/dspec1a.c
@@ -12,30 +12,165 @@
 
 /*This function finds the hessenberg form of a matrix A.*/
 
-#include "hess.h"
+#include "spec.h"
 #include <stdio.h>
 #include "string.h"
 #include "stdlib.h"
 #include "lapack.h"
 #include "matrixTranspose.h"
 #include "matrixMultiplication.h"
+#include "doubleComplex.h"
 
-void dspec1a(double *in1,double *in2,int size,int nout,double *out1,double *out2,double *out3,double *out4){
+extern int dggev_(char *,char *,int *,double *,int *,double *,int *,double *,double *,double *,double *,int *,double *,int *,double *,int *,int *);
+
+void assembleEigenvectorsInPlace(int N,double *ALPHAI,double *EVreal,double *EVimg){
+	int j,i;
+	j=0;
+	while(j<N){
+		if(ALPHAI[j] == 0){
+			//printf(" * ");
+			j+=1;	
+		}	
+		else{	
+			int ij;
+			int ij1;
+			for(i=0;i<N;i++){
+				ij = i+j*N;
+				ij1 = i+(j+1)*N;	
+				EVimg[ij] = EVreal[ij1];
+				EVimg[ij1] = -EVreal[ij1];	
+				EVreal[ij1] = EVreal[ij];
+			}
+			j+=2;						
+		}
+	}
+}
+
+void dspec1a(double *in1,double *in2,int size,int nout,doubleComplex *out1,double *out2,doubleComplex *out3,doubleComplex *out4){
+	
+	int i,j;
 	char JOBVL;
 	char JOBVR;
-	int N;
+	int N=size;
+	
 	double *A;
-	int LDA;
+	int LDA=N;
+	A = (double *)malloc(N*N*sizeof(double));	
+	memcpy(A,in1,N*N*sizeof(double));
+	
 	double *B;
-	int LDA;
+	int LDB=N;
+	B = (double *)malloc(N*N*sizeof(double));
+	memcpy(B,in2,N*N*sizeof(double));	
+	
 	double *ALPHAR;
+	ALPHAR = (double *)malloc(N*sizeof(double));
+	
 	double *ALPHAI;
+	ALPHAI = (double *)malloc(N*sizeof(double));
+		
 	double *BETA;
+	BETA = (double *)malloc(N*sizeof(double));
+	
 	double *VL;
-	int LDVL;
+	VL = (double *)malloc(N*N*sizeof(double));
+	
+	int LDVL=N;
+
 	double *VR;
-	int LDVR;
+	VR = (double *)malloc(N*N*sizeof(double));
+	int LDVR=N;
+	
+	int LWORK=8*N;
+
 	double *WORK;
-	int LWORK;
+	WORK = (double *)malloc(LWORK*sizeof(double));
+	
 	int INFO;
+	if(nout == 1){			/*out1 = spec(A,B)*/
+		JOBVL = 'N';
+		JOBVR = 'N';
+		dggev_(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&INFO);
+		for(i=0;i<N;i++){
+			out1[i] = DoubleComplex(ALPHAR[i]/BETA[i],ALPHAI[i]/BETA[i]);
+		}			
+	}
+	else if(nout == 2){		/*[out1,out2] = spec(A,B)*/
+		JOBVL = 'N';
+		JOBVR = 'N';
+		dggev_(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&INFO);
+		for(i=0;i<N;i++){
+			out1[i] = DoubleComplex(ALPHAR[i],ALPHAI[i]);		
+		}
+		memcpy(out2,BETA,N*sizeof(double));
+	}
+	else if(nout == 3){	 /* [out1,out2,out3] = spec(A,B) */
+		JOBVL = 'N';
+		JOBVR = 'V';	
+		dggev_(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&INFO);
+		for(i=0;i<N;i++){
+			out1[i] = DoubleComplex(ALPHAR[i],ALPHAI[i]);		
+		}	
+		memcpy(out2,BETA,N*sizeof(double));
+		
+		/*Because lapack routine doesn't give result in actual format, \
+			so we have to change the VR little-bit and then return the function */
+		
+		/*See the Scilab code || see the lapack subroutine libary - DGGEV where \
+			it is very explantory and explains all this.			
+		*/
+		double *EVimg;
+		EVimg = (double *)malloc(N*N*sizeof(double));
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				EVimg[i+j*N] = 0;					
+			}		
+		}		
+		assembleEigenvectorsInPlace(N,ALPHAI,VR,EVimg);		
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				out3[i+j*N] = DoubleComplex(VR[i+j*N],EVimg[i+j*N]);			
+			}		
+		}
+	}
+	else if(nout == 4){
+		JOBVL = 'V';
+		JOBVR = 'V';	
+
+		dggev_(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&INFO);
+
+		for(i=0;i<N;i++){
+			out1[i] = DoubleComplex(ALPHAR[i],ALPHAI[i]);		
+		}	
+
+		memcpy(out2,BETA,N*sizeof(double));
+		
+		double *EVimg;
+		EVimg = (double *)malloc(N*N*sizeof(double));
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				EVimg[i+j*N] = 0;					
+			}		
+		}		
+		assembleEigenvectorsInPlace(N,ALPHAI,VR,EVimg);		
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				out4[i+j*N] = DoubleComplex(VR[i+j*N],EVimg[i+j*N]);			
+			}		
+		}		
+		
+		double *EVimg1;
+		EVimg1 = (double *)malloc(N*N*sizeof(double));
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				EVimg1[i+j*N] = 0;			
+			}		
+		}		
+		assembleEigenvectorsInPlace(N,ALPHAI,VL,EVimg1);
+		for(i=0;i<N;i++){
+			for(j=0;j<N;j++){
+				out3[i+j*N] = DoubleComplex(VL[i+j*N],EVimg1[i+j*N]);			
+			}		
+		}
+	}
 }
diff --git a/2.3-1/src/c/string/disp/zdispa.c b/2.3-1/src/c/string/disp/zdispa.c
index 106775ef..bc71de4b 100644
--- a/2.3-1/src/c/string/disp/zdispa.c
+++ b/2.3-1/src/c/string/disp/zdispa.c
@@ -16,7 +16,7 @@ double zdispa (doubleComplex* in, int rows, int columns){
 	int i = 0,j = 0;
 	
 	for (i = 0; i < rows; ++i) {
-		for (j=0;j<columns;j++) printf("  %1.6lf + %1.6lfi  " ,zreals(in[i+j*rows]) ,zimags(in[i+j*rows])); 
+		for (j=0;j<columns;j++) printf("  %1.20lf + %1.20lfi  " ,zreals(in[i+j*rows]) ,zimags(in[i+j*rows])); 
 		printf("\n");
 	}
 	return 0;
-- 
cgit