summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--macros/ToolInitialization/INIT_FillSCI2LibCDirs.binbin1628308 -> 1640092 bytes
-rw-r--r--macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci49
-rw-r--r--macros/findDeps/getAllSources.binbin248996 -> 249160 bytes
-rw-r--r--macros/findDeps/getAllSources.sci1
-rw-r--r--src/c/linearAlgebra/includes/spec.h3
-rw-r--r--src/c/linearAlgebra/interfaces/int_spec.h7
-rw-r--r--src/c/linearAlgebra/spec/dspec1a.c151
-rw-r--r--src/c/string/disp/zdispa.c2
8 files changed, 204 insertions, 9 deletions
diff --git a/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin b/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin
index 55269e3..1c1a52c 100644
--- a/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin
+++ b/macros/ToolInitialization/INIT_FillSCI2LibCDirs.bin
Binary files differ
diff --git a/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci b/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci
index 008d7b3..c44ef3b 100644
--- a/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci
+++ b/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/macros/findDeps/getAllSources.bin b/macros/findDeps/getAllSources.bin
index 401fb99..7bc8fd0 100644
--- a/macros/findDeps/getAllSources.bin
+++ b/macros/findDeps/getAllSources.bin
Binary files differ
diff --git a/macros/findDeps/getAllSources.sci b/macros/findDeps/getAllSources.sci
index 5dced78..f0f883b 100644
--- a/macros/findDeps/getAllSources.sci
+++ b/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/src/c/linearAlgebra/includes/spec.h b/src/c/linearAlgebra/includes/spec.h
index b7b7d35..3cd0b18 100644
--- a/src/c/linearAlgebra/includes/spec.h
+++ b/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/src/c/linearAlgebra/interfaces/int_spec.h b/src/c/linearAlgebra/interfaces/int_spec.h
index c946329..8772a6a 100644
--- a/src/c/linearAlgebra/interfaces/int_spec.h
+++ b/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/src/c/linearAlgebra/spec/dspec1a.c b/src/c/linearAlgebra/spec/dspec1a.c
index 069454d..28440be 100644
--- a/src/c/linearAlgebra/spec/dspec1a.c
+++ b/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/src/c/string/disp/zdispa.c b/src/c/string/disp/zdispa.c
index 106775e..bc71de4 100644
--- a/src/c/string/disp/zdispa.c
+++ b/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;
generated by cgit (git 2.34.1) at 2024-12-20 07:31:49 +0000