diff options
-rw-r--r-- | 2.3-1/includes/sci2clib.h | 9 | ||||
-rw-r--r-- | 2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci | 156 | ||||
-rw-r--r-- | 2.3-1/macros/findDeps/getAllHeaders.sci | 2 | ||||
-rw-r--r-- | 2.3-1/macros/findDeps/getAllInterfaces.sci | 2 | ||||
-rw-r--r-- | 2.3-1/macros/findDeps/getAllSources.sci | 2 | ||||
-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 | 99 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c | 121 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c | 45 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c | 43 | ||||
-rw-r--r-- | 2.3-1/src/c/string/disp/ddispa.c | 2 |
16 files changed, 798 insertions, 81 deletions
diff --git a/2.3-1/includes/sci2clib.h b/2.3-1/includes/sci2clib.h index f662325c..c27a16db 100644 --- a/2.3-1/includes/sci2clib.h +++ b/2.3-1/includes/sci2clib.h @@ -18,6 +18,15 @@ extern "C" { /* interfacing lapack */ #include "lapack.h" + +#include "sva.h" +#include "int_sva.h" + + +#include "svd.h" +#include "int_svd.h" + + /* AUXILIARY FUNCTIONS */ /* interfacing abs */ diff --git a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci index 65859740..d310ad89 100644 --- a/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci +++ b/2.3-1/macros/ToolInitialization/INIT_FillSCI2LibCDirs.sci @@ -6933,6 +6933,162 @@ PrintStringInfo(' Adding Function: '+FunctionName+'.',GeneralReport,'file', INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCAnnFunDir,ClassName,GeneralReport,ExtensionCAnnFun);
INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCFLFunDir,ClassName,GeneralReport,ExtensionCFuncListFun);
+//------------------------------------
+//---- Class SVD -------------------
+//------------------------------------
+ClassName = 'SVD';
+
+// --- Class Annotation. ---
+PrintStringInfo(' Adding Class: '+ClassName+'.',GeneralReport,'file','y');
+ClassFileName = fullfile(SCI2CLibCAnnClsDir,ClassName+ExtensionCAnnCls);
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 1',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 2',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 1',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= ''1''',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 2',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 4',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(1)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(2)= ''1''',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 2',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 4',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).TP= ''d''',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(1)= ''1''',ClassFileName,'file','y');
+PrintStringInfo('OUT(4).SZ(2)= ''1''',ClassFileName,'file','y');
+
+ClassFileName = fullfile(SCI2CLibCFLClsDir,ClassName+ExtensionCFuncListCls);
+PrintStringInfo('d2'+ArgSeparator+'d2',ClassFileName,'file','y');
+PrintStringInfo('d2g2'+ArgSeparator+'d2d2d2',ClassFileName,'file','y');
+PrintStringInfo('d2'+ArgSeparator+'d2d2d2',ClassFileName,'file','y');
+PrintStringInfo('d2'+ArgSeparator+'d2d2d2d0',ClassFileName,'file','y');
+PrintStringInfo('d2d0'+ArgSeparator+'d2d2d2d0',ClassFileName,'file','y');
+
+// --- Annotation Function And Function List Function. ---
+FunctionName = 'svd';
+PrintStringInfo(' Adding Function: '+FunctionName+'.',GeneralReport,'file','y');
+INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCAnnFunDir,ClassName,GeneralReport,ExtensionCAnnFun);
+INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCFLFunDir,ClassName,GeneralReport,ExtensionCFuncListFun);
+
+//------------------------------------
+//---- Class SVA ---------------------
+//------------------------------------
+ClassName = 'SVA';
+
+// --- Class Annotation. ---
+PrintStringInfo(' Adding Class: '+ClassName+'.',GeneralReport,'file','y');
+ClassFileName = fullfile(SCI2CLibCAnnClsDir,ClassName+ExtensionCAnnCls);
+
+PrintStringInfo('NIN= 1',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+
+PrintStringInfo('NIN= 2',ClassFileName,'file','y');
+PrintStringInfo('NOUT= 3',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(1)= IN(1).SZ(1)',ClassFileName,'file','y');
+PrintStringInfo('OUT(1).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(1)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(2).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).TP= IN(1).TP',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(1)= IN(1).SZ(2)',ClassFileName,'file','y');
+PrintStringInfo('OUT(3).SZ(2)= FA_MIN(IN(1).SZ(1),IN(1).SZ(2))',ClassFileName,'file','y');
+
+ClassFileName = fullfile(SCI2CLibCFLClsDir,ClassName+ExtensionCFuncListCls);
+PrintStringInfo('d2'+ArgSeparator+'d2d2d2',ClassFileName,'file','y');
+PrintStringInfo('d2d0'+ArgSeparator+'d2d2d2',ClassFileName,'file','y');
+
+// --- Annotation Function And Function List Function. ---
+FunctionName = 'sva';
+PrintStringInfo(' Adding Function: '+FunctionName+'.',GeneralReport,'file','y');
+INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCAnnFunDir,ClassName,GeneralReport,ExtensionCAnnFun);
+INIT_GenAnnFLFunctions(FunctionName,SCI2CLibCFLFunDir,ClassName,GeneralReport,ExtensionCFuncListFun);
+
+
//------------------------------------
//---- Class lqr ---------------------
diff --git a/2.3-1/macros/findDeps/getAllHeaders.sci b/2.3-1/macros/findDeps/getAllHeaders.sci index 7788f2c4..81d83e08 100644 --- a/2.3-1/macros/findDeps/getAllHeaders.sci +++ b/2.3-1/macros/findDeps/getAllHeaders.sci @@ -196,6 +196,8 @@ function allHeaders = getAllHeaders(SharedInfo) "src/c/CACSD/includes/lqr.h" "src/c/CACSD/includes/lqe.h" "src/c/CACSD/includes/obscont.h" + "src/c/linearAlgebra/includes/sva.h" + "src/c/linearAlgebra/includes/svd.h" "src/c/linearAlgebra/includes/schur.h" "src/c/linearAlgebra/includes/balanc.h" "src/c/linearAlgebra/includes/rcond.h"]; diff --git a/2.3-1/macros/findDeps/getAllInterfaces.sci b/2.3-1/macros/findDeps/getAllInterfaces.sci index 9086ba10..fc601321 100644 --- a/2.3-1/macros/findDeps/getAllInterfaces.sci +++ b/2.3-1/macros/findDeps/getAllInterfaces.sci @@ -191,6 +191,8 @@ function allInterfaces = getAllInterfaces(SharedInfo) "src/c/CACSD/interfaces/int_lqr.h" "src/c/CACSD/interfaces/int_lqe.h" "src/c/CACSD/interfaces/int_obscont.h" + "src/c/linearAlgebra/interfaces/int_sva.h" + "src/c/linearAlgebra/interfaces/int_svd.h" "src/c/linearAlgebra/interfaces/int_schur.h" "src/c/linearAlgebra/interfaces/int_balanc.h" "src/c/linearAlgebra/interfaces/int_rcond.h"]; diff --git a/2.3-1/macros/findDeps/getAllSources.sci b/2.3-1/macros/findDeps/getAllSources.sci index 5f0fbff5..bb2718d6 100644 --- a/2.3-1/macros/findDeps/getAllSources.sci +++ b/2.3-1/macros/findDeps/getAllSources.sci @@ -1201,6 +1201,8 @@ function allSources = getAllSources(SharedInfo) "src/c/CACSD/lqr/dlqra.c" "src/c/CACSD/lqe/dlqea.c" "src/c/CACSD/obscont/dobsconta.c" + "src/c/linearAlgebra/sva/dsvaa.c" + "src/c/linearAlgebra/svd/dsvda.c" "src/c/linearAlgebra/schur/dschura.c" "src/c/linearAlgebra/schur/dgschura.c" "src/c/linearAlgebra/balanc/dbalanca.c" 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..7ff83645 100644 --- a/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "floatComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -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 = DoubleComplex(0, 0); - floatComplex temp,mm; + floatComplex accu = FloatComplex(0, 0); + 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<col;j++) { - y=row*j+i; - temp=in[y]*(cos(((M_PI)*(y+1-1./2.)*(x))/n)); + y=row*j+i; + vv = cos(((M_PI)*(y+1-1./2.)*(x))/n); + aa = FloatComplex(vv,0); + temp=cmuls(in[y],aa); out[x]=cadds(out[x],temp); } } if(x==0) - out[x]*=1./(sqrt(n)); + { + vv = 1./(sqrt(n)); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else { float res=2./n; - out[x]*=sqrt(res); + res = sqrt(res); + aa = FloatComplex(res,0); + out[x]=cmuls(out[x],aa); } } } @@ -76,29 +85,53 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) y=j*row+i; z=(float)(((float)j+1.0/2.0)*(float)v); q=(float)(M_PI/(float)col); - mm=in[y]*(cos(q*z)); + vv = cos(q*z); + aa = FloatComplex(vv,0); + mm=cmuls(in[y],aa); temp=cadds(temp,mm); } z=(float)(((float)i+1.0/2.0)*(float)u); q=(float)(M_PI/(float)row); - temp*=cos(q*z); + ff = cos(q*z); + bb = FloatComplex(ff,0); + temp=cmuls(temp,bb); out[x]=cadds(out[x],temp); } if(u==0) { - out[x]*=1./sqrt((float)row); + vv = 1./sqrt((float)row); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((float)col); + { + vv = 1./sqrt((float)col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } } else { - out[x]*=sqrt(2./row); + vv = sqrt(2./row); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((float)col); + { + vv = 1./sqrt((float)col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } } } } @@ -125,12 +158,14 @@ void cdcta(floatComplex *in,int row,int col,int sign,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)); } } } @@ -156,19 +191,37 @@ void cdcta(floatComplex *in,int row,int col,int sign,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/dct/zdcta.c b/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c index 0204d682..3ae2e333 100644 --- a/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "doubleComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -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<col;j++) { - y=row*j+i; - temp=in[y]*(cos(((M_PI)*(y+1-1./2.)*(x))/n)); + y=row*j+i; + vv = cos(((M_PI)*(y+1-1./2.)*(x))/n); + aa = DoubleComplex(vv,0); + temp=zmuls(in[y],aa); out[x]=zadds(out[x],temp); } } if(x==0) - out[x]*=1./(sqrt(n)); + { + vv = 1./(sqrt(n)); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else { - double res=2./n; - out[x]*=sqrt(res); + double res=2./n; + res = sqrt(res); + aa = DoubleComplex(res,0); + out[x]=zmuls(out[x],aa); } } } - } - else + } + else { n=col*row; for(u=0;u<row;u++) @@ -74,31 +83,55 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) for(j=0;j<col;j++) { y=j*row+i; - z=(double)(((double)j+1.0/2.0)*(double)v); - q=(double)(M_PI/(double)col); - mm=in[y]*(cos(q*z)); + z=(double )(((double )j+1.0/2.0)*(double )v); + q=(double )(M_PI/(double )col); + vv = cos(q*z); + aa = DoubleComplex(vv,0); + mm=zmuls(in[y],aa); temp=zadds(temp,mm); } - z=(double)(((double)i+1.0/2.0)*(double)u); - q=(double)(M_PI/(double)row); - temp*=cos(q*z); + z=(double )(((double )i+1.0/2.0)*(double )u); + q=(double )(M_PI/(double )row); + ff = cos(q*z); + bb = DoubleComplex(ff,0); + temp=zmuls(temp,bb); out[x]=zadds(out[x],temp); - } + } if(u==0) { - out[x]*=1./sqrt((double)row); + vv = 1./sqrt((double )row); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((double)col); + { + vv = 1./sqrt((double )col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } } else { - out[x]*=sqrt(2./row); + vv = sqrt(2./row); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((double)col); + { + vv = 1./sqrt((double )col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } } } } @@ -125,12 +158,14 @@ void zdcta(doubleComplex *in,int row,int col,int sign,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)); } } } @@ -155,20 +190,38 @@ void zdcta(doubleComplex *in,int row,int col,int sign,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)); + z=(double )(((double )v+1.0/2.0)*(double )j); + 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); + 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)); + } } } } 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..ae98ba19 100644 --- a/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c @@ -14,7 +14,8 @@ #include "idct.h" #include "addition.h" #include "types.h" -#include "floatComplex.h" +#include "floatComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -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 = DoubleComplex(0, 0); - floatComplex temp,mm; + floatComplex accu = FloatComplex(0, 0); + 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 <fftw3.h>*/ #include <math.h> @@ -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)); + } } } } 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; |