diff options
Diffstat (limited to 'src/c/linearAlgebra/svd/dsvda.c')
| -rw-r--r-- | src/c/linearAlgebra/svd/dsvda.c | 82 |
1 files changed, 70 insertions, 12 deletions
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){ |