summaryrefslogtreecommitdiff
path: root/src/c/linearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'src/c/linearAlgebra')
-rw-r--r--src/c/linearAlgebra/includes/svd.h2
-rw-r--r--src/c/linearAlgebra/interfaces/int_sva.h2
-rw-r--r--src/c/linearAlgebra/interfaces/int_svd.h10
-rw-r--r--src/c/linearAlgebra/sva/dsvaa.c18
-rw-r--r--src/c/linearAlgebra/svd/dsvda.c82
5 files changed, 88 insertions, 26 deletions
diff --git a/src/c/linearAlgebra/includes/svd.h b/src/c/linearAlgebra/includes/svd.h
index 260b87f..dea681f 100644
--- a/src/c/linearAlgebra/includes/svd.h
+++ b/src/c/linearAlgebra/includes/svd.h
@@ -18,7 +18,7 @@
extern "C" {
#endif
-void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \
+double dsvda(double tol,double *in1,int row,int col,double in2,double nout,double *out1, \
double *out2,double *out3);
void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *out1,\
doubleComplex *out2,doubleComplex *out3);
diff --git a/src/c/linearAlgebra/interfaces/int_sva.h b/src/c/linearAlgebra/interfaces/int_sva.h
index 4a2ec56..f1f8260 100644
--- a/src/c/linearAlgebra/interfaces/int_sva.h
+++ b/src/c/linearAlgebra/interfaces/int_sva.h
@@ -18,7 +18,7 @@ extern "C" {
#endif
#define d2svad2d2d2(in1,size,out1,out2,out3) dsvaa(1,in1,size[0],size[1],0,out1,out2,out3);
-#define d2d0svad2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvaa(2,in1,size1[0],size1[1],in2,out1,out2,out3);
+#define d2d0svad2d2d2(in1,size1,in2,out1,out2,out3) dsvaa(2,in1,size1[0],size1[1],in2,out1,out2,out3);
#ifdef __cplusplus
} /* extern "C" */
diff --git a/src/c/linearAlgebra/interfaces/int_svd.h b/src/c/linearAlgebra/interfaces/int_svd.h
index 449ee74..8f40bff 100644
--- a/src/c/linearAlgebra/interfaces/int_svd.h
+++ b/src/c/linearAlgebra/interfaces/int_svd.h
@@ -17,9 +17,13 @@
extern "C" {
#endif
-#define d2svdd2(in1,size1,out1) dsvda(in1,size1[0],size1[1],0,1,out1,NULL,NULL)
-#define d2g2svdd2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvda(in1,size1[0],size1[1],1,3,out1,out2,out3);
-#define d2svdd2d2d2(in1,size1,out1,out2,out3) dsvda(in1,size1[0],size1[1],0,3,out1,out2,out3);
+#define d2svdd2(in1,size1,out1) dsvda(0,in1,size1[0],size1[1],0,1,out1,NULL,NULL)
+#define d2g2svdd2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],1,3,out1,out2,out3);
+#define d2svdd2d2d2(in1,size1,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],0,3,out1,out2,out3);
+
+#define d2svdd2d2d2d0(in1,size1,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],0,4,out1,out2,out3);
+#define d2d0svdd2d2d2d0(in1,size1,tol,out1,out2,out3) dsvda(tol,in1,size1[0],size1[1],0,4,out1,out2,out3);
+
#define z2svdz2(in1,size1,out2) zsvda(in1,size1[0],size1[1],0,1,NULL,out2,NULL);
#define z2g2svdz2z2z2(in1,size1,in2,size2,out1,out2,out3) zsvda(in1,size1[0],size1[1],1,3,out1,out2,out3);
#define z2svdz2z2z2(in1,size1,out1,out2,out3) zsvda(in1,size1[0],size1[1],0,3,out1,out2,out3);
diff --git a/src/c/linearAlgebra/sva/dsvaa.c b/src/c/linearAlgebra/sva/dsvaa.c
index ee27eef..b7d07d8 100644
--- a/src/c/linearAlgebra/sva/dsvaa.c
+++ b/src/c/linearAlgebra/sva/dsvaa.c
@@ -33,14 +33,14 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \
/* Calculation of svd of a given matrix */
double *U,*S,*V;
- U = (double *)malloc((double)row*Min(row,col)*sizeof(double));
- S = (double *)malloc((double)Min(row,col)*Min(row,col)*sizeof(double));
- V = (double *)malloc((double)col*Min(row,col)*sizeof(double));
+ U = (double *)malloc((double)row*min(row,col)*sizeof(double));
+ S = (double *)malloc((double)min(row,col)*min(row,col)*sizeof(double));
+ V = (double *)malloc((double)col*min(row,col)*sizeof(double));
- dsvda(in1,M,N,1,3,U,S,V);
+ dsvda(0,in1,M,N,1,3,U,S,V);
if (ninp == 1){ /* [u,s,v] = sva(A) when input is only matrix */
- tol = Max(row,col)*S[0]*eps;
+ tol = max(row,col)*S[0]*eps;
rk = 0;
for(i=0;i<col;i++){
if(S[i+i*row] > tol){
@@ -52,7 +52,7 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \
tol = in2;
if(tol > 1){
rk = tol;
- if(rk > Min(row,col)){
+ if(rk > min(row,col)){
printf("ERROR: Wrong value for input argument !");
out1 = NULL;
out2 = NULL;
@@ -70,20 +70,20 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \
}
}
arow = M;
- acol = Min(M,N);
+ acol = min(M,N);
for(i=0;i<arow;i++){
for(j=0;j<rk;j++){
out1[i+j*row]=U[i+j*arow];
}
}
- arow = Min(M,N);
+ arow = min(M,N);
for(i=0;i<rk;i++){
for(j=0;j<rk;j++){
out2[i+j*(int)rk] = S[i+j*arow];
}
}
arow = N;
- acol = Min(M,N);
+ acol = min(M,N);
for(i=0;i<arow;i++){
for(j=0;j<rk;j++){
out3[i+j*arow] = V[i+j*arow];
diff --git a/src/c/linearAlgebra/svd/dsvda.c b/src/c/linearAlgebra/svd/dsvda.c
index e6af300..c3bcfc2 100644
--- a/src/c/linearAlgebra/svd/dsvda.c
+++ b/src/c/linearAlgebra/svd/dsvda.c
@@ -27,26 +27,40 @@ int max(int a,int b);
extern double dgesvd_(char*,char*,int*,int*,double*,int*,double*,double*,int*,\
double*,int*,double *,int*,int*);
+#define eps 2.22044604925e-16 /* pow(2,-52) */
+
/* DGESVD computes the singular value decomposition (SVD) of a real
M-by-N matrix A, optionally computing the left and/or right singular
vectors. The SVD is written
A = U * SIGMA * transpose(V) */
-void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \
- double *out2,double *out3){
+/*Function support -
+
+s=svd(X)
+[U,S,V]=svd(X)
+[U,S,V]=svd(X,0) (obsolete)
+[U,S,V]=svd(X,"e")
+[U,S,V,rk]=svd(X [,tol])
+*/
+
+double dsvda(double tol,double *in1,int row,int col,double in2,double nout,double *out1, \
+ double *out2,double *out3){
+
char JOBU,JOBVT;
- int j,k;
+ int i,j,k;
int LDU=1; /*Leading Dimension of U */
int LDVT=1; /*Leading Dimension of VT */
int M = row;
int N = col;
double *buf;
double *S,*U,*VT;
- double *WORK;
+ double *WORK;
- if((nout > 1 && in2 == 1) && (M != N)){ /* [U,S,VT] = svd(x,'e') */
+ int rk; /*Fourth output if needed */
+
+ /*if((nout > 1 && in2 == 1) && (M != N)){ // [U,S,VT] = svd(x,'e')
if(M > N){
JOBU = 'S';
JOBVT = 'A';
@@ -61,7 +75,7 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \
U = (double*) malloc((double) (LDU)*min(M,N)*sizeof(double));
VT = (double*) malloc((double) (LDVT)*N*sizeof(double));
}
- else if(nout > 1){ /* [U,S,VT = svd(x)] */
+ else */if(nout > 1){ /* [U,S,VT = svd(x)] */
JOBU = 'A'; /*If JOBU = 'A', U contains the M-by-M orthogonal matrix U */
JOBVT = 'A'; /*JOBVT = 'A': all N rows of V**T are returned in the array VT;*/
LDU = M;
@@ -74,7 +88,7 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \
JOBVT = 'N';
}
int LDA = max(1,M);
-
+
/* Making a copy of input matrix */
buf = (double*) malloc((double)M*N*sizeof(double));
memcpy(buf,in1,M*N*sizeof(double));
@@ -100,19 +114,63 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \
if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j);
else *((out2+j*(min(M,N)))+k) = 0;
}
- }
- dtransposea(VT,LDVT,N,out3);
+ }
+
+ //dtransposea(VT,LDVT,N,out3);
+ /*As there is some patch of error in SVD, these lines are added */
+
+ for(j=1;j<=N;j++){
+ for(i=j;i<=N;i++){
+ *(out3+i+(j-1)*N-1) = VT[j+(i-1)*N-1];
+ *(out3+j+(i-1)*N-1) = VT[i+(j-1)*N-1];
+ }
+ }
+ /*for(i=0;i<N;i++){
+ for(j=0;j<N;j++){
+ printf("%lf ",VT[i*row+j]);
+ }
+ printf("\n");
+ }*/
}
else{
- memcpy(out1,U,LDU*min(row,col)*sizeof(double));
+ memcpy(out1,U,M*min(M,N)*sizeof(double));
for(j=0;j<min(M,N);j++){
for(k=0;k<min(M,N);k++){
if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j);
else *((out2+j*(min(M,N)))+k) = 0;
}
}
- dtransposea(VT,LDVT,N,out3);
- }
+ //dtransposea(VT,LDVT,N,out3);
+ /*As there is some patch of error in DGESVD, these lines are added */
+ /* out3 first taken in some array then will be copied from it. */
+ double *outV;
+ outV = (double *)malloc(N*N*sizeof(double));
+ for(j=1;j<=N;j++){
+ for(i=j;i<=N;i++){
+ *(outV+i+(j-1)*N-1) = VT[j+(i-1)*N-1];
+ *(outV+j+(i-1)*N-1) = VT[i+(j-1)*N-1];
+ }
+ }
+
+ for(j=0;j<min(M,N)*N;j++){
+ *(out3+j) = *(outV+j);
+ }
+ }
+
+ /* From the fortran file of scilab code - if(tol.eq.0.0d0) tol=dble(max(M,N))*eps*stk(lSV) */
+ if(tol == 0){
+ tol = (double)max(M,N)*eps*S[0];
+ }
+ if(nout == 4){ /*[U,S,VT,rk] = svd(X,tol) where tol - tolerance*/
+ rk = 0;
+ for(i=0;i<min(M,N);i++){
+ if(S[i] > tol){
+ rk = i+1;
+ }
+ }
+ return rk;
+ }
+ return 0;
}
int min(int a,int b){
generated by cgit (git 2.34.1) at 2024-12-20 19:39:32 +0000