summaryrefslogtreecommitdiff
path: root/2.3-1/src/c
diff options
context:
space:
mode:
Diffstat (limited to '2.3-1/src/c')
-rw-r--r--2.3-1/src/c/linearAlgebra/includes/spec.h3
-rw-r--r--2.3-1/src/c/linearAlgebra/interfaces/int_spec.h7
-rw-r--r--2.3-1/src/c/linearAlgebra/spec/dspec1a.c151
-rw-r--r--2.3-1/src/c/string/disp/zdispa.c2
4 files changed, 154 insertions, 9 deletions
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;
generated by cgit (git 2.34.1) at 2024-12-20 15:54:31 +0000