diff options
Diffstat (limited to '2.3-1/src/c')
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/includes/sva.h | 29 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/includes/svd.h | 30 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/interfaces/int_sva.h | 28 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/interfaces/int_svd.h | 35 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/sva/dsvaa.c | 92 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/svd/dsvda.c | 184 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c | 2 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c | 2 | ||||
-rw-r--r-- | 2.3-1/src/c/string/disp/ddispa.c | 2 |
9 files changed, 401 insertions, 3 deletions
diff --git a/2.3-1/src/c/linearAlgebra/includes/sva.h b/2.3-1/src/c/linearAlgebra/includes/sva.h new file mode 100644 index 00000000..ea628a32 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/sva.h @@ -0,0 +1,29 @@ + /* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __SVA_H__ +#define __SVA_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ + double *out2,double *out3); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__SVA_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/includes/svd.h b/2.3-1/src/c/linearAlgebra/includes/svd.h new file mode 100644 index 00000000..dea681fc --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/svd.h @@ -0,0 +1,30 @@ + /* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ +#ifndef __SVD_H__ +#define __SVD_H__ +#include "types.h" +#include "doubleComplex.h" + +#ifdef __cplusplus +extern "C" { +#endif + +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); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__SVD_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h b/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h new file mode 100644 index 00000000..f1f8260a --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h @@ -0,0 +1,28 @@ + /* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_SVA_H__ +#define __INT_SVA_H__ + +#ifdef __cplusplus +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,out1,out2,out3) dsvaa(2,in1,size1[0],size1[1],in2,out1,out2,out3); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_SVA_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h b/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h new file mode 100644 index 00000000..8f40bffe --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h @@ -0,0 +1,35 @@ + /* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_SVD_H__ +#define __INT_SVD_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#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); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_SVD_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/sva/dsvaa.c b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c new file mode 100644 index 00000000..b7d07d8c --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c @@ -0,0 +1,92 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + + */ +#include "sva.h" +#include "svd.h" +#include "lapack.h" +#include <stdio.h> +#include <stdlib.h> +#include "string.h" +#include "matrixTranspose.h" + +#define eps 2.22044604925e-16 + +void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ + double *out2,double *out3){ + + double tol; + double rk=0; + int N=col,M=row; + if(row == 0 && col == 0) return; + int i,j; + int arow; /*Actual row of given matrix*/ + int acol; /*Actual col of given matrix*/ + + /* 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)); + + 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; + rk = 0; + for(i=0;i<col;i++){ + if(S[i+i*row] > tol){ + rk+=1; + } + } + } + else{ /*[u,s,v] = sva(A,tol) when two input's are there */ + tol = in2; + if(tol > 1){ + rk = tol; + if(rk > min(row,col)){ + printf("ERROR: Wrong value for input argument !"); + out1 = NULL; + out2 = NULL; + out3 = NULL; + return; + } + } + else{ + rk = 0; + for(i=0;i<col;i++){ + if(S[i+i*row] > tol){ + rk+=1; + } + } + } + } + arow = M; + 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); + 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); + for(i=0;i<arow;i++){ + for(j=0;j<rk;j++){ + out3[i+j*arow] = V[i+j*arow]; + } + } +} diff --git a/2.3-1/src/c/linearAlgebra/svd/dsvda.c b/2.3-1/src/c/linearAlgebra/svd/dsvda.c new file mode 100644 index 00000000..c3bcfc29 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/svd/dsvda.c @@ -0,0 +1,184 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Sandeep Gupta + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + + */ + +/*Funtion to find singular value decomposition of given matrix */ + +#include "lapack.h" +#include <stdio.h> +#include <stdlib.h> +#include "string.h" +#include <math.h> +#include "svd.h" +#include "matrixTranspose.h" + +int min(int a,int b); +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) */ + +/*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 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; + + 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'; + LDVT = N; + } + else{ + JOBU = 'A'; + JOBVT = 'S'; + LDVT = min(M,N); + } + LDU = M; + 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)] */ + 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; + LDVT = N; + U = (double*) malloc((double) M*M*sizeof(double)); + VT = (double*) malloc((double) N*N*sizeof(double)); + } + else{ /* ans = svd(x) */ + JOBU = 'N'; + 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)); + + S = (double*)malloc((double)min(col,row)*sizeof(double)); + + int LWORK = 5*min(M,N); + WORK = (double*)malloc((double)LWORK*sizeof(double)); + int INFO = 0; /*For successful exit */ + + dgesvd_(&JOBU,&JOBVT,&M,&N,buf,&LDA,S,U,&LDU,VT,&LDVT,WORK,&LWORK,&INFO); + /*Subroutine DGESVD from Lapack lib. */ + + if (nout == 1){ /* ans = svd(x)*/ + memcpy(out1,S,min(row,col)*sizeof(double)); + //printf("%lf %lf %lf",*(S),*(S+1),*(S+2)); + } /* [U,S,VT] = svd(x) */ + else if(in2 == 0 && nout > 1){ + memcpy(out1,U,LDU*M*sizeof(double)); + //memcpy(out3,VT,LDVT*min(row,col)*sizeof(double)); + for(j=0;j<M;j++){ + for(k=0;k<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); + /*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,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); + /*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){ + if(a > b) return b; + return a; +} + +int max(int a,int b){ + if(a > b) return a; + return b; +} diff --git a/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c b/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c index 5bc27929..bd15e5b9 100644 --- a/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c @@ -26,7 +26,7 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) int x,y; float res,ress; float re,z,q,m; - floatComplex accu = DoubleComplex(0, 0); + floatComplex accu = FloatComplex(0, 0); floatComplex temp,mm; if(sign==-1) { diff --git a/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c b/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c index ec0df0c7..e6c746ce 100644 --- a/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c @@ -26,7 +26,7 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) int x,y; float res,ress; float re,z,q,m; - floatComplex accu = DoubleComplex(0, 0); + floatComplex accu = FloatComplex(0, 0); floatComplex temp,mm; if(row==1) { diff --git a/2.3-1/src/c/string/disp/ddispa.c b/2.3-1/src/c/string/disp/ddispa.c index 5e6bb84e..29aaceeb 100644 --- a/2.3-1/src/c/string/disp/ddispa.c +++ b/2.3-1/src/c/string/disp/ddispa.c @@ -16,7 +16,7 @@ double ddispa (double* in, int rows, int columns){ int i = 0,j = 0; for (i = 0; i < rows; ++i) { - for (j=0;j<columns;j++) printf (" %1.20f ", in[i+j*rows]); + for (j=0;j<columns;j++) printf (" %e ", in[i+j*rows]); printf("\n"); } return 0; |