From 771262669202239258b54394840cc75114f9cb2b Mon Sep 17 00:00:00 2001 From: ukashanoor Date: Tue, 4 Jul 2017 11:55:14 +0530 Subject: after debugging 1 --- 2.3-1/src/c/linearAlgebra/includes/sva.h | 29 ++++ 2.3-1/src/c/linearAlgebra/includes/svd.h | 30 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_sva.h | 28 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_svd.h | 35 ++++ 2.3-1/src/c/linearAlgebra/sva/dsvaa.c | 92 +++++++++++ 2.3-1/src/c/linearAlgebra/svd/dsvda.c | 184 +++++++++++++++++++++ .../src/c/signalProcessing/transforms/dct/cdcta.c | 2 +- .../c/signalProcessing/transforms/idct/cidcta.c | 2 +- 2.3-1/src/c/string/disp/ddispa.c | 2 +- 9 files changed, 401 insertions(+), 3 deletions(-) create mode 100644 2.3-1/src/c/linearAlgebra/includes/sva.h create mode 100644 2.3-1/src/c/linearAlgebra/includes/svd.h create mode 100644 2.3-1/src/c/linearAlgebra/interfaces/int_sva.h create mode 100644 2.3-1/src/c/linearAlgebra/interfaces/int_svd.h create mode 100644 2.3-1/src/c/linearAlgebra/sva/dsvaa.c create mode 100644 2.3-1/src/c/linearAlgebra/svd/dsvda.c (limited to '2.3-1/src') 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 +#include +#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 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 tol){ + rk+=1; + } + } + } + } + arow = M; + acol = min(M,N); + for(i=0;i +#include +#include "string.h" +#include +#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 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*/ #include @@ -24,10 +25,10 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) int i,j,k,u,v; int n; int x,y; - float res,ress; + float res,ress,vv,ff; float re,z,q,m; floatComplex accu = FloatComplex(0, 0); - floatComplex temp,mm; + floatComplex temp,mm,aa,bb,cc; if(sign==-1) { if(row==1) @@ -43,17 +44,25 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) { for(j=0;j*/ #include @@ -24,10 +25,10 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) int i,j,k,u,v; int n; int x,y; - double res,ress; - double re,z,q,m; + double res,ress,vv,ff; + double re,z,q,m; doubleComplex accu = DoubleComplex(0, 0); - doubleComplex temp,mm; + doubleComplex temp,mm,aa,bb,cc; if(sign==-1) { if(row==1) @@ -43,22 +44,30 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) { for(j=0;j*/ #include @@ -24,10 +25,10 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) int i,j,k,u,v; int n=col; int x,y; - float res,ress; + float res,ress,vv,ff; float re,z,q,m; floatComplex accu = FloatComplex(0, 0); - floatComplex temp,mm; + floatComplex temp,mm,aa,bb; if(row==1) { res=1./sqrt(n); @@ -46,12 +47,14 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa=FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa=FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } } } @@ -77,19 +80,37 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) y=row*j+i; mm=in[j*row+i]; z=(float)(((float)v+1.0/2.0)*(float)j); - q=(float)(M_PI/(float)col); - mm=mm*(cos(q*z)); + q=(float)(M_PI/(float)col); + vv=cos(q*z); + aa=FloatComplex(vv,0); + mm=cmuls(mm,aa); if(j==0) - temp=cadds(temp,mm*(1./sqrt((float)col))); + { + vv=1./sqrt((float)col); + aa=FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } else - temp=cadds(temp,mm*sqrt(2./col)); + { + vv=sqrt(2./col); + aa=FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } } z=(float)(((float)u+1.0/2.0)*(float)i); q=(float)(M_PI/(float)row); if(i==0) - out[x]=cadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv=(cos(z*q))*(1./sqrt(row)); + aa=FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } else - out[x]=cadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv=(cos(z*q))*sqrt(2./row); + aa=FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } } } } diff --git a/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c b/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c index 2177b18c..cc01c966 100644 --- a/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "doubleComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include */ #include @@ -24,10 +25,10 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) int i,j,k,u,v; int n=col; int x,y; - double res,ress; + double res,ress,vv,ff; double re,z,q,m; doubleComplex accu = DoubleComplex(0, 0); - doubleComplex temp,mm; + doubleComplex temp,mm,aa,bb; if(row==1) { res=1./sqrt(n); @@ -46,17 +47,19 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa=DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa=DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } } } } - + } } else @@ -77,19 +80,37 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) y=row*j+i; mm=in[j*row+i]; z=(double)(((double)v+1.0/2.0)*(double)j); - q=(double)(M_PI/(double)col); - mm=mm*(cos(q*z)); + q=(double)(M_PI/(double)col); + vv=cos(q*z); + aa=DoubleComplex(vv,0); + mm=zmuls(mm,aa); if(j==0) - temp=zadds(temp,mm*(1./sqrt((double)col))); + { + vv=1./sqrt((double)col); + aa=DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } else - temp=zadds(temp,mm*sqrt(2./col)); + { + vv=sqrt(2./col); + aa=DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } } z=(double)(((double)u+1.0/2.0)*(double)i); q=(double)(M_PI/(double)row); if(i==0) - out[x]=zadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv=(cos(z*q))*(1./sqrt(row)); + aa=DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } else - out[x]=zadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv=(cos(z*q))*sqrt(2./row); + aa=DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } } } } -- cgit