From b43eccd4cffed5bd1017c5821524fb6e49202f78 Mon Sep 17 00:00:00 2001 From: Sandeep Gupta Date: Sun, 18 Jun 2017 23:55:40 +0530 Subject: First commit --- 2.3-1/src/c/linearAlgebra/balanc/dbalanca.c | 75 ++++++++++ 2.3-1/src/c/linearAlgebra/bdiag/dbdiaga.c | 23 +++ 2.3-1/src/c/linearAlgebra/hess/dhessa.c | 81 ++++++++++ 2.3-1/src/c/linearAlgebra/includes/balanc.h | 27 ++++ 2.3-1/src/c/linearAlgebra/includes/hess.h | 26 ++++ 2.3-1/src/c/linearAlgebra/includes/rcond.h | 26 ++++ 2.3-1/src/c/linearAlgebra/includes/schur.h | 30 ++++ 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_balanc.h | 29 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_hess.h | 28 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_rcond.h | 25 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_schur.h | 60 ++++++++ 2.3-1/src/c/linearAlgebra/interfaces/int_sva.h | 28 ++++ 2.3-1/src/c/linearAlgebra/interfaces/int_svd.h | 31 ++++ 2.3-1/src/c/linearAlgebra/rcond/drconda.c | 46 ++++++ 2.3-1/src/c/linearAlgebra/schur/dgschura.c | 161 ++++++++++++++++++++ 2.3-1/src/c/linearAlgebra/schur/dschura.c | 126 ++++++++++++++++ 2.3-1/src/c/linearAlgebra/schur/zgschura.c | 164 ++++++++++++++++++++ 2.3-1/src/c/linearAlgebra/sva/dsvaa.c | 92 ++++++++++++ 2.3-1/src/c/linearAlgebra/svd/.1.c.swp | Bin 0 -> 12288 bytes 2.3-1/src/c/linearAlgebra/svd/dsvda.c | 126 ++++++++++++++++ 2.3-1/src/c/linearAlgebra/svd/zsvda.c | 173 ++++++++++++++++++++++ 23 files changed, 1436 insertions(+) create mode 100644 2.3-1/src/c/linearAlgebra/balanc/dbalanca.c create mode 100644 2.3-1/src/c/linearAlgebra/bdiag/dbdiaga.c create mode 100644 2.3-1/src/c/linearAlgebra/hess/dhessa.c create mode 100644 2.3-1/src/c/linearAlgebra/includes/balanc.h create mode 100644 2.3-1/src/c/linearAlgebra/includes/hess.h create mode 100644 2.3-1/src/c/linearAlgebra/includes/rcond.h create mode 100644 2.3-1/src/c/linearAlgebra/includes/schur.h 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_balanc.h create mode 100644 2.3-1/src/c/linearAlgebra/interfaces/int_hess.h create mode 100644 2.3-1/src/c/linearAlgebra/interfaces/int_rcond.h create mode 100644 2.3-1/src/c/linearAlgebra/interfaces/int_schur.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/rcond/drconda.c create mode 100644 2.3-1/src/c/linearAlgebra/schur/dgschura.c create mode 100644 2.3-1/src/c/linearAlgebra/schur/dschura.c create mode 100644 2.3-1/src/c/linearAlgebra/schur/zgschura.c create mode 100644 2.3-1/src/c/linearAlgebra/sva/dsvaa.c create mode 100644 2.3-1/src/c/linearAlgebra/svd/.1.c.swp create mode 100644 2.3-1/src/c/linearAlgebra/svd/dsvda.c create mode 100644 2.3-1/src/c/linearAlgebra/svd/zsvda.c (limited to '2.3-1/src/c/linearAlgebra') diff --git a/2.3-1/src/c/linearAlgebra/balanc/dbalanca.c b/2.3-1/src/c/linearAlgebra/balanc/dbalanca.c new file mode 100644 index 00000000..a86a1967 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/balanc/dbalanca.c @@ -0,0 +1,75 @@ +/* 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: Siddhesh Wani + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +/*Funtion to balance given square matrix to improve its condition number*/ +#include "balanc.h" +#include "lapack.h" +#include "stdlib.h" +#include "string.h" +#include "eye.h" + +void dbalanca(double* in1, int rows, double* in2, double* out1, \ + double* out2, double* out3, double* out4) +{ + char JOB = 'B'; + char RSIDE = 'R'; + char LSIDE = 'L'; + + int ILO, IHI, INFO; + double *buf1, *buf2, *LSCALE, *RSCALE; + double *LWORK; + + + if(in2 == NULL) + { + /*Single input matrix*/ + buf1 = (double*) malloc((double) rows*rows*sizeof(double)); + LSCALE = (double*) malloc((double) rows*sizeof(double)); + + /*copy input to temp buf as lapack function modifies input*/ + memcpy(buf1,in1,rows*rows*sizeof(double)); + + dgebal_(&JOB,&rows,buf1,&rows,&ILO,&IHI,LSCALE,&INFO); + deyea(out2,rows,rows); + dgebak_(&JOB,&RSIDE,&rows,&ILO,&IHI,LSCALE,&rows,out2,&rows,&INFO); + + /*copy computed matrix to output*/ + memcpy(out1,buf1,rows*rows*sizeof(double)); + } + else + { + /*two matrices as inputs*/ + buf1 = (double*) malloc((double) rows*rows*sizeof(double)); + buf2 = (double*) malloc((double) rows*rows*sizeof(double)); + LSCALE = (double*) malloc((double) rows*sizeof(double)); + RSCALE = (double*) malloc((double) rows*sizeof(double)); + LWORK = (double*) malloc((double) 6*rows*sizeof(double)); + + /*copy input to temp buf as lapack function modifies input*/ + memcpy(buf1,in1,rows*rows*sizeof(double)); + memcpy(buf2,in2,rows*rows*sizeof(double)); + + dggbal_(&JOB,&rows,buf1,&rows,buf2,&rows,&ILO,&IHI,LSCALE,RSCALE, \ + LWORK,&INFO); + deyea(out3,rows,rows); + deyea(out4,rows,rows); + + dgebak_(&JOB,&LSIDE,&rows,&ILO,&IHI,LSCALE,&rows,out3,&rows,&INFO); + dgebak_(&JOB,&RSIDE,&rows,&ILO,&IHI,LSCALE,&rows,out4,&rows,&INFO); + + /*copy computed matrix to output*/ + memcpy(out1,buf1,rows*rows*sizeof(double)); + memcpy(out2,buf2,rows*rows*sizeof(double)); + + } + +} \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/bdiag/dbdiaga.c b/2.3-1/src/c/linearAlgebra/bdiag/dbdiaga.c new file mode 100644 index 00000000..fd455ba7 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/bdiag/dbdiaga.c @@ -0,0 +1,23 @@ +/* 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 + */ + +/*This function performs the block-diagonalization of matrix A.*/ + +/*--------------------------------------------------------------------------*/ +/* [Ab [,X [,bs]]]=bdiag(A [,rMax]) */ +/*--------------------------------------------------------------------------*/ + +#include +#include "string.h" +#include "stdlib.h" +#include "lapack.h" + diff --git a/2.3-1/src/c/linearAlgebra/hess/dhessa.c b/2.3-1/src/c/linearAlgebra/hess/dhessa.c new file mode 100644 index 00000000..57f81b35 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/hess/dhessa.c @@ -0,0 +1,81 @@ +/* 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 + */ + +/*This function finds the hessenberg form of a matrix A.*/ + +#include "hess.h" +#include +#include "string.h" +#include "stdlib.h" +#include "lapack.h" +#include "matrixTranspose.h" +#include "matrixMultiplication.h" + +extern int dgehrd_(int *, int *,int *,double *,int *,double *,double *,int *,int *); +extern int dorghr_(int *, int *,int *,double *,int *,double *,double *,int *,int *); + + +void dhessa(double *in1,int size,int nout,double *out1, double *out2){ + int i,j,k; + int N = size; + int ILO=1; + int IHI=N; + double *A; + int LDA=N; + double *TAU; + double *WORK; + int LWORK = N; + int INFO; + A = (double *)malloc((double)size*size*sizeof(double)); + memcpy(A,in1,size*size*sizeof(double)); + TAU = (double *)malloc((double)size*sizeof(double)); + WORK = (double *)malloc((double)LWORK*sizeof(double)); + dgehrd_(&N,&ILO,&IHI,A,&LDA,TAU,WORK,&N,&INFO); + + for(i=0;i 1){ + dorghr_(&N,&ILO,&IHI,A,&LDA,TAU,WORK,&LWORK,&INFO); + for(i=0;i=ILO-1;i--){ + tau = TAU[i]; + double V[size],v[size],v1[size*size]; + for(j=0;j + +/*flag --> 0: nothing + --> 1: continuous + --> 2: discrete +*/ + +lapack_logical selctg21( double* in1, double* in2, double* in3); +lapack_logical selctg22( double* in1, double* in2, double* in3); + +double dgschura(double* in1, int size, double* in2, int flag, int nout, \ + double* out1, double* out2, double* out3, double* out4) +{ + char JOBVSL = 'N'; + char JOBVSR = 'N'; + char SORT = 'N'; + int SDIM = 0; + int LDVSL = size, LDVSR = size; + int LWORK = 8*size+16, INFO; + double *ALPHAR, *ALPHAI, *BETA, *VSL, *VSR, *WORK; + int *BWORK; + double ret = 0; + double *buf1, *buf2; /*input is copied to buf, since lapack function directly + modifies the input variable*/ + + /*Used incase of flag > 0*/ + LAPACK_D_SELECT3 selctg = &selctg21; + + if(nout >= 2){ + JOBVSL = 'V'; + JOBVSR = 'V'; + } + if(flag > 0) SORT = 'S'; + + buf1 = (double*) malloc((double) size*size*sizeof(double)); + buf2 = (double*) malloc((double) size*size*sizeof(double)); + ALPHAR = (double*) malloc((double) size*sizeof(double)); + ALPHAI = (double*) malloc((double) size*sizeof(double)); + BETA = (double*) malloc((double) size*sizeof(double)); + VSL = (double*) malloc((double) LDVSL*size*sizeof(double)); + VSR = (double*) malloc((double) LDVSR*size*sizeof(double)); + WORK = (double*) malloc((double) LWORK*sizeof(double)); + BWORK = (int*) malloc((double) size*sizeof(double)); + + + memcpy(buf1,in1,size*size*sizeof(double)); + memcpy(buf2,in2,size*size*sizeof(double)); + + dgges_(&JOBVSL,&JOBVSR,&SORT,selctg,&size,buf1,&size,buf2,&size,&SDIM, \ + ALPHAR,ALPHAI,BETA,VSL,&LDVSL,VSR,&LDVSR,WORK,&LWORK,BWORK,&INFO); + + /*if (INFO != 0) + { + out1 = NULL; + return 0; + }*/ + + if(nout == 1) + { + return(SDIM); + } + else if(nout == 2) + { + if(flag == 0) + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(double)); + memcpy(out2,buf2,size*size*sizeof(double)); + } + else + { + /*copy VSR to out1 and return SDIM */ + memcpy(out1,VSR,size*size*sizeof(double)); + ret = SDIM; + } + } + else if(nout == 3) + { + /*copy VSL to out1, VSR to out2, return SDIM*/ + memcpy(out1,VSL,size*size*sizeof(double)); + memcpy(out2,VSR,size*size*sizeof(double)); + ret = SDIM; + } + else if(nout == 4) + { + if(flag == 0) + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(double)); + memcpy(out2,buf2,size*size*sizeof(double)); + /*copy VSL to out3 and VSR to out4*/ + memcpy(out3,VSL,size*size*sizeof(double)); + memcpy(out4,VSR,size*size*sizeof(double)); + } + else + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(double)); + memcpy(out2,buf2,size*size*sizeof(double)); + /*copy VSR to out3 and return SDIM*/ + memcpy(out3,VSR,size*size*sizeof(double)); + ret = SDIM; + } + } + else /*nout = 5*/ + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(double)); + memcpy(out2,buf2,size*size*sizeof(double)); + /*copy VSL to out3 and VSR to out4*/ + memcpy(out3,VSL,size*size*sizeof(double)); + memcpy(out4,VSR,size*size*sizeof(double)); + /*return SDIM*/ + ret = SDIM; + } + + free(buf1); + free(buf2); + free(VSL); + free(VSR); + free(ALPHAR); + free(ALPHAI); + free(BETA); + free(WORK); + free(BWORK); + + return ret; +} + +lapack_logical selctg21(double* in1, double* in2, double* in3) +{ + if((sqrt(*in1**in1+*in2**in2)/ *in3) < 1) + return 1; + else + return 0; +} + +lapack_logical selctg22(double* in1, double* in2, double* in3) +{ + if((sqrt(*in1**in1+*in2**in2)/ *in3) < 1) + return 1; + else + return 0; +} \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/schur/dschura.c b/2.3-1/src/c/linearAlgebra/schur/dschura.c new file mode 100644 index 00000000..802caa81 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/schur/dschura.c @@ -0,0 +1,126 @@ +/* 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: Siddhesh Wani + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +/*Funtion to find schur decomposition of given square matrix */ +#include "schur.h" +#include "lapack.h" +#include "stdlib.h" +#include "string.h" +#include + +/*flag --> 0: nothing + --> 1: continuous + --> 2: discrete +*/ + +lapack_logical selctg11( double* in1, double* in2); +lapack_logical selctg12( double* in1, double* in2); + +double dschura(double* in1, int size, int flag, int nout, double* out1, \ + double* out2) +{ + char JOBVS = 'N'; + char SORT = 'N'; + int SDIM = 0; + int LVDS = size; + int LWORK = 3*size, INFO; + double *WR, *WI, *VS, *WORK; + int *BWORK; + double ret = 0; + double *buf; /*input is copied to buf, since lapack function direclty + modifies the input variable*/ + + + /*Used incase of flag > 0*/ + LAPACK_D_SELECT2 selctg; + + if(flag == 1 || flag == 0) + selctg = &selctg11; + else if(flag == 2) + selctg = &selctg12; + + if(nout >= 2) JOBVS = 'V'; + if(flag > 0) SORT = 'S'; + + buf = (double*) malloc((double) size*size*sizeof(double)); + WR = (double*) malloc((double) size*sizeof(double)); + WI = (double*) malloc((double) size*sizeof(double)); + VS = (double*) malloc((double) LVDS*size*sizeof(double)); + WORK = (double*) malloc((double) LWORK*sizeof(double)); + BWORK = (int*) malloc((double) size*sizeof(double)); + + + memcpy(buf,in1,size*size*sizeof(double)); + + dgees_(&JOBVS,&SORT,selctg,&size,buf,&size,&SDIM,WR,WI,VS,&LVDS, \ + WORK,&LWORK,BWORK,&INFO); + + /*if (INFO != 0) + { + out1 = NULL; + return 0; + } */ + + if(nout == 1) + { + /*Copy result in in1 to out1*/ + memcpy(out1,buf,size*size*sizeof(double)); + } + else if(nout == 2) + { + if(flag == 0) + { + /*copy in1 to out2 and VS to out1*/ + memcpy(out2,buf,size*size*sizeof(double)); + memcpy(out1,VS,size*size*sizeof(double)); + } + else + { + /*copy VS to out1 and SDIM to out2*/ + memcpy(out1,VS,size*size*sizeof(double)); + ret = SDIM; + } + } + else + { + /*copy VS to out1, SDIM to out2, in1 to out3*/ + memcpy(out1,VS,size*size*sizeof(double)); + memcpy(out2,buf,size*size*sizeof(double)); + ret = SDIM; + } + + free(buf); + free(WI); + free(WR); + free(VS); + free(WORK); + free(BWORK); + + return ret; +} + +lapack_logical selctg11(double* in1, double* in2) +{ + if(*in1 <= 0) + return 1; + else + return 0; +} + +lapack_logical selctg12(double* in1, double* in2) +{ + + if(sqrt(*in1**in1+*in2**in2) < 1) + return 1; + else + return 0; +} \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/schur/zgschura.c b/2.3-1/src/c/linearAlgebra/schur/zgschura.c new file mode 100644 index 00000000..a49cbe51 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/schur/zgschura.c @@ -0,0 +1,164 @@ +/* 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 + */ + +/*Fucntion to find generalised schur decomposition of given square complex matrices */ + +#include "schur.h" +#include "lapack.h" +#include "stdlib.h" +#include "string.h" +#include +#include "doubleComplex.h" +#include "matrixTranspose.h" + +/*flag --> 0: nothing + --> 1: continuous + --> 2: discrete +*/ + +lapack_logical selctg21( doubleComplex* in1, doubleComplex* in2, doubleComplex* in3); +lapack_logical selctg22( doubleComplex* in1, doubleComplex* in2, doubleComplex* in3); + +doubleComplex dgschura(doubleComplex* in1, int size, doubleComplex* in2, int flag, int nout, \ + doubleComplex* out1, doubleComplex* out2, doubleComplex* out3, doubleComplex* out4) +{ + char JOBVSL = 'N'; + char JOBVSR = 'N'; + char SORT = 'N'; + int SDIM = 0; + int LDVSL = size, LDVSR = size; + int LWORK = 2*size, INFO; + doubleComplex *ALPHA, *BETA, *VSL, *VSR, *WORK; + double *RWORK; + int *BWORK; + doubleComplex ret = 0; + doubleComplex *buf1, *buf2; /*input is copied to buf, since lapack function directly + modifies the input variable*/ + + /*Used incase of flag > 0*/ + LAPACK_D_SELECT3 selctg = &selctg21; + + if(nout >= 2){ + JOBVSL = 'V'; + JOBVSR = 'V'; + } + if(flag > 0) SORT = 'S'; + + buf1 = (doubleComplex*) malloc((doubleComplex) size*size*sizeof(doubleComplex)); + buf2 = (doubleComplex*) malloc((doubleComplex) size*size*sizeof(doubleComplex)); + ALPHA = (doubleComplex*) malloc((doubleComplex) size*sizeof(doubleComplex)); + RWORK = (double *)malloc((double)8*size*sizeof(double)); + BETA = (doubleComplex*) malloc((doubleComplex) size*sizeof(doubleComplex)); + VSL = (doubleComplex*) malloc((doubleComplex) LDVSL*size*sizeof(doubleComplex)); + VSR = (doubleComplex*) malloc((doubleComplex) LDVSR*size*sizeof(doubleComplex)); + WORK = (doubleComplex*) malloc((doubleComplex) LWORK*sizeof(doubleComplex)); + BWORK = (int*) malloc((doubleComplex) size*sizeof(doubleComplex)); + + + memcpy(buf1,in1,size*size*sizeof(doubleComplex)); + memcpy(buf2,in2,size*size*sizeof(doubleComplex)); + + zgges_(&JOBVSL,&JOBVSR,&SORT,selctg,&size,buf1,&size,buf2,&size,&SDIM, \ + ALPHA,BETA,VSL,&LDVSL,VSR,&LDVSR,WORK,&LWORK,RWORK,BWORK,&INFO); + + /*if (INFO != 0) + { + out1 = NULL; + return 0; + }*/ + + if(nout == 1) + { + return(SDIM); + } + else if(nout == 2) + { + if(flag == 0) + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(doubleComplex)); + memcpy(out2,buf2,size*size*sizeof(doubleComplex)); + } + else + { + /*copy VSR to out1 and return SDIM */ + memcpy(out1,VSR,size*size*sizeof(doubleComplex)); + ret = SDIM; + } + } + else if(nout == 3) + { + /*copy VSL to out1, VSR to out2, return SDIM*/ + memcpy(out1,VSL,size*size*sizeof(doubleComplex)); + memcpy(out2,VSR,size*size*sizeof(doubleComplex)); + ret = SDIM; + } + else if(nout == 4) + { + if(flag == 0) + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(doubleComplex)); + memcpy(out2,buf2,size*size*sizeof(doubleComplex)); + /*copy VSL to out3 and VSR to out4*/ + memcpy(out3,VSL,size*size*sizeof(doubleComplex)); + memcpy(out4,VSR,size*size*sizeof(doubleComplex)); + } + else + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(doubleComplex)); + memcpy(out2,buf2,size*size*sizeof(doubleComplex)); + /*copy VSR to out3 and return SDIM*/ + memcpy(out3,VSR,size*size*sizeof(doubleComplex)); + ret = SDIM; + } + } + else /*nout = 5*/ + { + /*copy in1 to out1 and in2 to out2*/ + memcpy(out1,buf1,size*size*sizeof(doubleComplex)); + memcpy(out2,buf2,size*size*sizeof(doubleComplex)); + /*copy VSL to out3 and VSR to out4*/ + memcpy(out3,VSL,size*size*sizeof(doubleComplex)); + memcpy(out4,VSR,size*size*sizeof(doubleComplex)); + /*return SDIM*/ + ret = SDIM; + } + + free(buf1); + free(buf2); + free(VSL); + free(VSR); + free(ALPHA); + free(BETA); + free(WORK); + free(BWORK); + + return ret; +} + +lapack_logical selctg21(doubleComplex* in1, doubleComplex* in2, doubleComplex* in3) +{ + if((sqrt(*in1**in1+*in2**in2)/ *in3) < 1) + return 1; + else + return 0; +} + +lapack_logical selctg22(doubleComplex* in1, doubleComplex* in2, doubleComplex* in3) +{ + if((sqrt(*in1**in1+*in2**in2)/ *in3) < 1) + return 1; + else + return 0; +} 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..ee27eef7 --- /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(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*); + +/* 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){ + + char JOBU,JOBVT; + int 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; + + 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 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/linearAlgebra/svd/zsvda.c b/2.3-1/src/c/linearAlgebra/svd/zsvda.c new file mode 100644 index 00000000..0d360222 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/svd/zsvda.c @@ -0,0 +1,173 @@ +/* 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 "svd.h" +#include "lapack.h" +#include +#include +#include "string.h" +#include "doubleComplex.h" +#include "matrixTranspose.h" +#include "conj.h" + +extern doubleComplex zgesvd_( char* , char* , int* , int* ,doubleComplex *,\ + int* , double* ,doubleComplex* , int* ,doubleComplex* , int* ,\ + doubleComplex* , int* , double* , int* ); + + +int Min(int a,int b){ + if(a > b) + return b; + return a; +} + +int Max(int a,int b){ + if(a > b) + return a; + else + return b; +} + +void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *out1,\ + doubleComplex *out2,doubleComplex *out3){ + + /* Allocating memory and copying the input in buf*/ + doubleComplex *buf; + buf = (doubleComplex *)malloc(row*col*sizeof(doubleComplex)); + memcpy(buf,in1,row*col*sizeof(doubleComplex)); + + /* Type of variable used */ + int i,j,k; + char JOBU,JOBVT; + int M = row; + int N = col; + int LDA,LDU,LDVT,LWORK,INFO; + + /*double precision array to store Sigma*/ + double *S; + S = (double *)malloc(Min(M,N)*sizeof(double)); + + /* amount of memory needed for work */ + LWORK = Max(1,2*Min(M,N)+Max(M,N)); + doubleComplex *WORK = malloc(Max(1,2*LWORK)*sizeof(doubleComplex)); + + double *RWORK; + RWORK = (double *)malloc(5*Min(M,N)*sizeof(double)); + + INFO = 0; + + if(nout == 1){ + JOBU = 'N'; + JOBVT = 'N'; + LDA = M; + LDU = M; + LDVT = N; + //doubleComplex *U,*VT; + //U = malloc(sizeof(doubleComplex)); + //VT = malloc(sizeof(doubleComplex)); + zgesvd_(&JOBU,&JOBVT,&M,&N,buf,&LDA,S,NULL,&LDU,NULL,&LDVT,WORK,&LWORK,RWORK,&INFO); + + //memcpy(out2,S,Min(M,N)*sizeof(double)); + for(i=0;i N){ + JOBU = 'S'; + JOBVT = 'A'; + LDVT = N; + } + else{ + JOBU = 'A'; + JOBVT = 'S'; + LDVT = Min(M,N); + } + doubleComplex *U; + U = malloc(LDU*Min(M,N)*sizeof(doubleComplex)); + doubleComplex *VT; + VT = malloc(LDVT*N*sizeof(doubleComplex)); + zgesvd_(&JOBU,&JOBVT,&M,&N,buf,&LDA,S,U,&LDU,VT,&LDVT,WORK,&LWORK,RWORK,&INFO); + memcpy(out1,U,M*Min(M,N)*sizeof(doubleComplex)); + //ztransposea(VT,LDVT,Min(row,col),out3); + + /* These lines are added to patch an error of ZGESVD */ + /* + ij = i+(j-1)*N + ji = j+(i-1)*N + zstk(lV+ij-1) = conjg(zstk(lVT+ji-1)) + zstk(lV+ji-1) = conjg(zstk(lVT+ij-1)) + */ + for(i=0;i +#include "spec.h" + + + +void cspeca(floatComplex* in, int rows, floatComplex* out){ + /* As we use Lapack to find the eigenvalues, we must cast the floatComplex input into doubleComplex + and the doubleComplex output of dspeca into floatComplex*/ + + int i; + doubleComplex* dblin; + doubleComplex* dblout; + + dblin=(doubleComplex*)malloc((unsigned int)(rows*rows)*sizeof(doubleComplex)); + dblout=(doubleComplex*)malloc((unsigned int)rows*sizeof(doubleComplex)); + + for (i=0;i +#include "string.h" +#include "stdlib.h" +#include "lapack.h" +#include "matrixTranspose.h" +#include "matrixMultiplication.h" + +void dspec1a(double *in1,double *in2,int size,int nout,double *out1,double *out2,double *out3,double *out4){ + char JOBVL; + char JOBVR; + int N; + double *A; + int LDA; + double *B; + int LDA; + double *ALPHAR; + double *ALPHAI; + double *BETA; + double *VL; + int LDVL; + double *VR; + int LDVR; + double *WORK; + int LWORK; + int INFO; +} diff --git a/2.3-1/src/c/linearAlgebra/spec/dspeca.c b/2.3-1/src/c/linearAlgebra/spec/dspeca.c new file mode 100644 index 00000000..1131ecf2 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/dspeca.c @@ -0,0 +1,71 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * 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 + * + */ + +#include +#include "spec.h" +#include "lapack.h" +#include "zeros.h" +#include "max.h" + +void dspeca(double* in, int rows,double* out){ + int i=0, j=0; + int symmetric=0; + int INFO=0; + int iWorkSize = 0; + double* pdblWork; + double* outReal; + double* outImag; + double* pdblLeftvectors; + double* pdblRightvectors; + double* inCopy; + + inCopy = (double*)malloc((unsigned int)(rows*rows) * sizeof(double)); + outReal = (double*)malloc((unsigned int)rows * sizeof(double)); + outImag = (double*)malloc((unsigned int)rows * sizeof(double)); + pdblLeftvectors=NULL; + pdblRightvectors=NULL; + + + iWorkSize = 3*rows; + pdblWork = (double*)malloc((unsigned int)iWorkSize * sizeof(double)); + + for(i=0;i +#include "spec.h" + + + +void sspeca(float* in, int rows, float* out){ + /* As we use Lapack to find the eigenvalues, we must cast the float input into double + and the doubleComplex output of dspeca into floatComplex*/ + + int i; + double* dblin; + double* dblout; + + dblin=(double*)malloc((unsigned int)(rows*rows)*sizeof(double)); + dblout=(double*)malloc((unsigned int)rows*sizeof(double)); + + for (i=0;i +#include "spec.h" +#include "stdio.h" +#include "assert.h" +#include "math.h" + +static void dspecaTest(void){ + double in[4]={1,1,1,3}; + double resultR[2]={0.5857864376269050765700,3.4142135623730949234300}; + + + double *out; + + int i; + + out=(double*)malloc((unsigned int)2*sizeof(double)); + + dspeca(in,2,out); + for(i=0;i<2;i++){ + if (out[i]>1e-16) assert( fabs(out[i]-resultR[i]) / fabs(out[i]) <3e-16); + else assert(1); + } + +} + + + + +static void zspecaTest(void){ + double inR[4]={1,1,1,3}; + double inI[4]={0,0,0,0}; + double resultR[2]={0.5857864376269050765700,3.4142135623730949234300}; + double resultI[2]={0,0}; + + double in2R[4]={1,1,-2,3}; + double in2I[4]={0,0,0,0}; + double result2R[2]={1.9999999999999997779554,1.9999999999999997779554}; + double result2I[2]={0.9999999999999997779554,-0.9999999999999997779554}; + + double in3R[16]={0.0603054538369178771973,0.631347181741148233414 ,0.0241030259057879447937,0.1527438252232968807221, + 0.9481177683919668197632,0.2744265575893223285675,0.4794727200642228126526,0.485576837789267301559, + 0.676400367170572280884,0.8321249918080866336823,0.0125875836238265037537,0.5453779897652566432953, + 0.8426716942340135574341,0.9029165101237595081329,0.4409482078626751899719,0.8332359003834426403046 }; + double in3I[16]={0.7233976423740386962891,0.4377150186337530612946,0.3080607382580637931824,0.8749813153408467769623, + 0.5355882328003644943237 ,0.3085998897440731525421,0.3354632416740059852600,0.2342486302368342876434, + 0.2589411698281764984131,0.8521509231068193912506,0.4821739485487341880798,0.6095217890106141567230, + 0.9872822705656290054321,0.9811426815576851367950,0.3303113700821995735169,0.3589145573787391185761}; + double result3R[4]={1.9486046375762748894545,0.2315060459861970343365,- 0.7694743345806510648188,- 0.2300808535483104266817}; + double result3I[4]={2.1242015857341254303492,- 0.4115101262891645017561,- 0.1980716835789462781925,0.3584662621795701720195}; + + doubleComplex *in,*in2,*in3,out[2],out2[4]; + + int i; + + in=DoubleComplexMatrix(inR,inI,4); + in2=DoubleComplexMatrix(in2R,in2I,4); + in3=DoubleComplexMatrix(in3R,in3I,16); + + zspeca(in,2,out); + for(i=0;i<2;i++){ + if (zreals(out[i])>1e-16) assert( fabs(zreals(out[i])-resultR[i]) / fabs(zreals(out[i])) <3e-16); + else assert(1); + if (zimags(out[i])>1e-16) assert( fabs(zimags(out[i])-resultI[i]) / fabs(zimags(out[i])) <1e-16); + else assert(1); + } + + + zspeca(in2,2,out); + for(i=0;i<2;i++){ + if (zreals(out[i])>1e-16) assert( fabs(zreals(out[i])-result2R[i]) / fabs(zreals(out[i])) <3e-16); + else assert(1); + if (zimags(out[i])>1e-16) assert( fabs(zimags(out[i])-result2I[i]) / fabs(zimags(out[i])) <3e-15); + else assert(1); + } + + + zspeca(in3,4,out2); + for(i=0;i<4;i++){ + if (zreals(out2[i])>1e-16) assert( fabs(zreals(out2[i])-result3R[i]) / fabs(zreals(out2[i])) <3e-15); + else assert(1); + if (zimags(out2[i])>1e-16) assert( fabs(zimags(out2[i])-result3I[i]) / fabs(zimags(out2[i])) <3e-15); + else assert(1); + } + + +} + + +static int testSpec(void){ + printf(">>> Double Spec Tests <<<\n"); + printf(">>> Double <<<\n"); + dspecaTest(); + printf(">>> DoubleComplex <<<\n"); + zspecaTest(); + return 0; +} + + +int main(void){ + assert (testSpec()==0); + + + return 0; +} diff --git a/2.3-1/src/c/linearAlgebra/spec/testFloatSpec.c b/2.3-1/src/c/linearAlgebra/spec/testFloatSpec.c new file mode 100644 index 00000000..9bd46adc --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/testFloatSpec.c @@ -0,0 +1,110 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * 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 + * + */ + +#include +#include "spec.h" +#include "stdio.h" +#include "assert.h" +#include "math.h" + +static void sspecaTest(void){ + float in[4]={1.0f,1.0f,1.0f,3.0f}; + float resultR[2]={0.5857864376269050765700f,3.4142135623730949234300f}; + + float in1[4]={0.0f,4.0f,1.0f,0.0f}; + float result[2]={2.0f,-2.0f}; + + float in2[4]={1.0f,1.0f,-2.0f,3.0f}; + + float *out; + + int i; + + out=(float *)malloc((unsigned int)2*sizeof(float)); + + sspeca(in,2,out); + for (i=0;i<2;i++) printf("%f\n",out[i]); + for(i=0;i<2;i++){ + if (out[i]>1e-16) assert( fabs(out[i]-resultR[i]) / fabs(out[i]) <1e-16); + else assert(1); + } + + + sspeca(in1,2,out); + for (i=0;i<2;i++) printf("%f\n",out[i]); + for(i=0;i<2;i++){ + if (out[i]>1e-16) assert( fabs(out[i]-result[i]) / fabs(out[i]) <1e-16); + else assert(1); + } + + sspeca(in2,2,out); + for (i=0;i<2;i++) printf("%f\n",out[i]); + +} + + +static void cspecaTest(void){ + float inR[4]={1.0f,1.0f,1.0f,3.0f}; + float inI[4]={0.0f,0.0f,0.0f,0.0f}; + float resultR[2]={0.5857864376269050765700f,3.4142135623730949234300f}; + float resultI[2]={0,0}; + + float in2R[4]={1.0f,1.0f,-2.0f,3.0f}; + float in2I[4]={0.0f,0.0f,0.0f,0.0f}; + float result2R[2]={1.9999999999999997779554f,1.9999999999999997779554f}; + float result2I[2]={0.9999999999999997779554f,-0.9999999999999997779554f}; + + floatComplex *in,*in2,out[4]; + + int i; + + in=FloatComplexMatrix(inR,inI,4); + in2=FloatComplexMatrix(in2R,in2I,4); + + cspeca(in,2,out); + for(i=0;i<2;i++){ + if (creals(out[i])>1e-16) assert( fabs(creals(out[i])-resultR[i]) / fabs(creals(out[i])) <1e-15); + else assert(1); + if (cimags(out[i])>1e-16) assert( fabs(cimags(out[i])-resultI[i]) / fabs(cimags(out[i])) <1e-16); + else assert(1); + } + + + cspeca(in2,2,out); + for(i=0;i<2;i++){ + if (creals(out[i])>1e-16) assert( fabs(creals(out[i])-result2R[i]) / fabs(creals(out[i])) <1e-15); + else assert(1); + if (cimags(out[i])>1e-16) assert( fabs(cimags(out[i])-result2I[i]) / fabs(cimags(out[i])) <1e-15); + else assert(1); + } + +} + + + + +static int testSpec(void){ + printf(">>> Float Spec Tests <<<\n"); + printf(">>> Float <<<\n"); + sspecaTest(); + printf(">>> FloatComplex <<<\n"); + cspecaTest(); + return 0; +} + + +int main(void){ + assert (testSpec()==0); + + + return 0; +} diff --git a/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj b/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj new file mode 100644 index 00000000..f696bff4 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj @@ -0,0 +1,178 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {651B0E61-1047-4575-BE31-D9CB28062CBA} + testDoubleSpec + Win32Proj + + + + Application + Unicode + true + + + Application + Unicode + true + + + Application + Unicode + + + Application + Unicode + + + + + + + + + + + + + + + + + + + <_ProjectFileVersion>10.0.40219.1 + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + true + true + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + false + false + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + EditAndContinue + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + MachineX86 + + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + MachineX86 + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + + + + + CompileAsCpp + CompileAsCpp + CompileAsCpp + CompileAsCpp + + + + + {9b1bd750-1fef-4d6b-9422-782d16181cee} + false + + + {fd335544-52bf-4736-a34e-77f591d158d5} + false + + + + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj.filters b/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj.filters new file mode 100644 index 00000000..dcdaea5c --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/test_DoubleSpec/testDoubleSpec.vcxproj.filters @@ -0,0 +1,22 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj b/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj new file mode 100644 index 00000000..0631bb40 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj @@ -0,0 +1,178 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {B7BEF8A2-3E42-4FA3-94A7-22618594057D} + testFloatSpec + Win32Proj + + + + Application + Unicode + true + + + Application + Unicode + true + + + Application + Unicode + + + Application + Unicode + + + + + + + + + + + + + + + + + + + <_ProjectFileVersion>10.0.40219.1 + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + true + true + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + false + false + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + EditAndContinue + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + MachineX86 + + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + MachineX86 + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + + + + + CompileAsCpp + CompileAsCpp + CompileAsCpp + CompileAsCpp + + + + + {9b1bd750-1fef-4d6b-9422-782d16181cee} + false + + + {fd335544-52bf-4736-a34e-77f591d158d5} + false + + + + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj.filters b/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj.filters new file mode 100644 index 00000000..2ba2e59f --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/test_FloatSpec/testFloatSpec.vcxproj.filters @@ -0,0 +1,22 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec/zspeca.c b/2.3-1/src/c/linearAlgebra/spec/zspeca.c new file mode 100644 index 00000000..a2b1c25a --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec/zspeca.c @@ -0,0 +1,80 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * 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 + * + */ + +#include +#include "spec.h" +#include "lapack.h" +#include "zeros.h" +#include "conj.h" + +void zspeca(doubleComplex* in, int rows,doubleComplex* out){ + int i = 0, j = 0; + int hermitian = 0; + int INFO = 0; + int iWorkSize = 0; + doubleComplex* pdblWork; + doubleComplex* pdblLeftvectors; + doubleComplex* pdblRightvectors; + doubleComplex* pdblRWork; + double* outReal; + double* outImag; + doubleComplex* inCopy; + + inCopy = (doubleComplex*)malloc((unsigned int)(rows*rows) * sizeof(doubleComplex)); + outReal = (double*)malloc((unsigned int)rows * sizeof(double)); + outImag = (double*)malloc((unsigned int)rows * sizeof(double)); + pdblLeftvectors=NULL; + pdblRightvectors=NULL; + + iWorkSize = 2*rows; + pdblWork = (doubleComplex*)malloc((unsigned int)iWorkSize * sizeof(doubleComplex)); + pdblRWork = (doubleComplex*)malloc((unsigned int)(3*rows) * sizeof(doubleComplex)); + + for(i=0;i +#include "spec.h" + + + +void cspec2a(floatComplex* in, int rows, floatComplex* eigenvalues,floatComplex* eigenvectors){ + /* As we use Lapack to find the eigenvalues, we must cast the floatComplex input into doubleComplex + and the doubleComplex outputs of dspeca into floatComplex*/ + + int i; + doubleComplex* dblin; + doubleComplex* dbleigenvalues,*dbleigenvectors; + + dblin=(doubleComplex*)malloc((unsigned int)(rows*rows)*sizeof(doubleComplex)); + dbleigenvalues=(doubleComplex*)malloc((unsigned int)(rows*rows)*sizeof(doubleComplex)); + dbleigenvectors=(doubleComplex*)malloc((unsigned int)(rows*rows)*sizeof(doubleComplex)); + + for (i=0;i +#include "spec.h" +#include "lapack.h" +#include "zeros.h" +#include "max.h" + +void dspec2a(double* in, int rows,double* eigenvalues,double* eigenvectors){ + int i=0, j=0, ij=0, ij1=0; + int symmetric=0; + int INFO=0; + int iWorkSize = 0; + double* pdblWork; + double* outReal; + double* outImag; + double* pdblLeftvectors; + double* pdblRightvectors; + double* inCopy; + + /* FIXME : malloc here */ + inCopy = (double*)malloc((unsigned int)(rows*rows) * sizeof(double)); + outReal = (double*)malloc((unsigned int)rows * sizeof(double)); + outImag = NULL; + pdblLeftvectors=NULL; + pdblRightvectors=NULL; + + + iWorkSize = 4*rows; + pdblWork = (double*)malloc((unsigned int)iWorkSize * sizeof(double)); + + for(i=0;i +#include "spec.h" + + + +void sspec2a(float* in, int rows, float* eigenvalues,float* eigenvectors){ + /* As we use Lapack to find the eigenvalues, we must cast the float input into double + and the doubleComplex outputs of dspec2a into floatComplex*/ + + int i; + double* dblin; + double *dbleigenvalues,*dbleigenvectors; + + dblin=(double*)malloc((unsigned int)(rows*rows)*sizeof(double)); + dbleigenvalues = (double*)malloc((unsigned int)(rows*rows)*sizeof(double)); + dbleigenvectors = (double*)malloc((unsigned int)(rows*rows)*sizeof(double)); + + for (i=0;i1e-16) assert( fabs(out3[i]-resultValues3R[i]) / fabs(out3[i]) <3e-16); + else assert(1); + } + for(i=0;i<9;i++){ + if (out4[i]>1e-16) assert( fabs(out4[i]-resultVectors3R[i]) / fabs(out4[i]) <3e-16); + else assert(1); + } + + + dspec2a(in,2,out1,out2); + for(i=0;i<4;i++) printf("%f\n",out1[i]); + for(i=0;i<4;i++){ + if (out1[i]>1e-16) assert( fabs(out1[i]-resultValuesR[i]) / fabs(out1[i]) <3e-16); + else assert(1); + } + for(i=0;i<4;i++){ + if (out2[i]>1e-16) assert( fabs(out2[i]-resultVectorsR[i]) / fabs(out2[i]) <3e-16); + else assert(1); + } + + + dspec2a(in2,2,out1,out2); + + for(i=0;i<4;i++){ + if (out1[i]>1e-16) assert( fabs(out1[i]-resultValues2R[i]) / fabs(out1[i]) <1e-16); + else assert(1); + } + for(i=0;i<4;i++){ + if (out2[i]>1e-16) assert( fabs(out2[i]-resultVectors2R[i]) / fabs(out2[i]) <3e-16); + else assert(1); + } + +} + + + + +static void zspec2aTest(void){ + double inR[4]={1,1,1,3}; + double inI[4]={0,0,0,0}; + double resultValuesR[4]={0.5857864376269050765700,0,0,3.4142135623730949234300}; + double resultValuesI[4]={0,0,0,0}; + double resultVectorsR[4]={- 0.9238795325112867384831,0.3826834323650897817792, + 0.3826834323650897817792,0.9238795325112867384831}; + double resultVectorsI[4]={0,0,0,0}; + + double in2R[4]={1,1,-2,3}; + double in2I[4]={0,0,0,0}; + double resultValues2R[4]={1.9999999999999997779554,0,0,1.9999999999999997779554}; + double resultValues2I[4]={0.9999999999999997779554,0,0,-0.9999999999999997779554}; + double resultVectors2R[4]={0.8164965809277261454824,- 0.4082482904638631282523, + 0.8164965809277261454824,- 0.4082482904638631282523}; + double resultVectors2I[4]={0,- 0.4082482904638629062077,0,0.4082482904638629062077}; + + + double in3R[16]={ 0.4685934986919164657593, 0.4262020816095173358917, 0.4217656338587403297424, 0.1034853602759540081024 , + 0.4279759414494037628174, 0.7860729382373392581940, 0.8568082069978117942810, 0.1993437460623681545258 , + 0.7131301630288362503052, 0.5208952468819916248322, 0.9311723159626126289368, 0.4143836158327758312225, + 0.5980196148157119750977, 0.5549105503596365451813, 0.8552952585741877555847, 0.3097750707529485225678}; + double in3I[16]={0.9446128141134977340698, 0.2442680452950298786163, 0.8760446915403008460999, 0.4874215493910014629364, + 0.3844018988311290740967, 0.0922345430590212345123, 0.0705418726429343223572, 0.7338807261548936367035, + 0.7792180571705102920532, 0.7827638057060539722443, 0.2161567779257893562317, 0.3420197847299277782440, + 0.430093832314014434814 , 0.6261752569116652011871, 0.9554250938817858695984, 0.45415506651625037193 }; + double resultValues3R[16]={2.0556597713281337114211 , 0. , 0. , 0. , + 0. , - 0.3025728286579054682193 , 0. , 0. , + 0. , 0. , 0.1432776981289703988054 , 0. , + 0. , 0. , 0. , 0.59924918284561956571821}; + double resultValues3I[16]={2.0772175543447914947137 , 0. , 0. , 0. , + 0. , - 0.5223086999196728807959 , 0. , 0. , + 0. , 0. , 0.3871735516160564882782 , 0. , + 0. , 0. , 0. , - 0.2349232044266161556489 }; + double resultVectors3R[16]={ 0.5762578366606958546470 , 0.4529873647084416976583 , 0.5438116534400605495137 , 0.3599855200589946457868 , + - 0.3285074902800487195442, - 0.0923654569732548225147, 0.6473744750124794666846, - 0.1812302216765642182139 , + 0.7395301236826927748780, - 0.4588526195219964631811, 0.0997916953568371567496, - 0.3603263850824164338249 , + 0.0717862182160008133192, 0.6967198304315765922112, - 0.1153149977500917411355 , - 0.3232658081022470875077 }; + double resultVectors3I[16]={ 0. , - 0.1082260232744242933745, - 0.1403164362264210929254 , 0.0775054168635924967123 , + 0.0347047096401541868560, - 0.4780157518240965019984, 0. , 0.4493438228998393735303, + 0. , - 0.0484367977162360741072 , 0.0664282497469102783949 , - 0.3098340648435350952461 , + 0.1768623397230159322024, 0. , - 0.4621680558325051979551 , 0.3830799998945915163517}; + + + double in4R[49]={0.0292230211198329925537, 0.1069206790998578071594, 9.7463708464056253433228, 7.7042609406635165214539 , 6.3966313377022743225098, 9.5358861843124032020569, 4.1484833415597677230835, + 4.6615396952256560325623, 2.2384551353752613067627, 3.4569733263924717903137, 2.3406236339360475540161, 6.6193414805456995964050,2.6911795139312744140625,9.8302489006891846656799 , + 7.3933512251824140548706, 2.9021358629688620567322 , 3.9959496073424816131592 , 9.792278115637600421906 , 0.0755135808140039443970, 7.0210226578637957572937 , 8.9075061306357383728027 , + 6.5269742021337151527405 , 4.136228552088141441345 , 8.5816909139975905418396 , 2.3960896767675876617432 , 4.8200417729094624519348, 9.848585547879338264465 , 4.4697216479107737541199, + 2.3298137634992599487305, 2.2345139319077134132385, 8.0025654565542936325073, 7.7507343282923102378845, 4.8250066302716732025146, 7.8448831336572766304016 , 1.5346793178468942642212 , + 7.7228771103546023368835 , 8.9520217850804328918457, 9.4673257926478981971741, 6.9818257447332143783569, 9.5591608202084898948669, 0.4840173013508319854736, 5.3235206427052617073059, + 2.9468670953065156936646, 6.9771366892382502555847 , 9.3762038648128509521484 , 1.0755608463659882545471 , 7.7634243946522474288940, 6.3734178384765982627869 , 0.6630615703761577606201 }; + double in4I[49]={8.1694598542526364326477, 3.5473910067230463027954, 5.860544512979686260223, 3.4610254690051078796387, 4.4247510144487023353577, 8.8226650562137365341187, 7.6436930662021040916443 , + 7.764139864593744277954, 2.8076809318736195564270, 6.9086006004363298416138, 7.1747286943718791007996, 3.7201813608407974243164 , 9.3254965776577591896057, 2.2566775511950254440308, + 7.5211019208654761314392 , 3.0747506581246852874756 , 5.9793261485174298286438, 1.9238903466612100601196, 2.3870888305827975273132 , 7.4708331003785133361816, 9.7131536761298775672913, + 4.3703553732484579086304, 6.2248750543221831321716, 3.2913279719650745391846, 3.2963873865082859992981, 2.7686371374875307083130 , 4.8729835031554102897644, 7.7320465445518493652344 , + 9.0775218093767762184143, 5.0122931879013776779175, 5.1595458528026938438416, 5.79267887398600578308 , 1.5463936375454068183899 , 7.6111377868801355361938, 6.9079177780076861381531 , + 2.8732293471693992614746, 9.64253133628517389297 , 5.6607243325561285018921, 5.7821379369124770164490, 2.6634209789335727691650 , 4.7470985027030110359192, 6.0735465306788682937622, + 0.4097307054325938224793, 5.0125684589147567749023 , 5.2959309751167893409729, 6.2594583164900541305542, 2.2193526616320013999939 , 3.467419948428869247437, 3.9511676924303174018860 }; + double resultValues4R[49]={39.018457752476393807228 ,0,0,0,0,0,0, + 0,2.2455352415739073812517,0,0,0,0,0, + 0,0,- 8.4848591657394489828903,0,0,0,0, + 0,0,0,- 9.617033196016185669919,0,0,0, + 0,0,0,0,- 3.7815852014367528077798,0,0, + 0,0,0,0,0,- 4.44858427601393291440,0, + 0,0,0,0,0,0,- 0.3001282122401837670900}; + double resultValues4I[49]={36.525725390582877594170,0,0,0,0,0,0, + 0,- 5.9058762755733944516123,0,0,0,0,0, + 0,0,- 3.5959527108686115681735,0,0,0,0, + 0,0,0,1.8928051444081777088257,0,0,0, + 0,0,0,0,- 3.1889864758923365251064,0,0, + 0,0,0,0,0,3.2418188738592150777151,0, + 0,0,0,0,0,0,1.5279802073147508156836}; + + double resultVectors4R[49]={ 0.3499995148567048852684, 0.3131141942665225941234, 0.4518966551544936205431, 0.3429358568207974133912, 0.3073919828504245721490, 0.4236868694783915145763, 0.3928739146959387973368 , + 0.2654822689469810148566 , - 0.3607918093136298631762, 0.0473515544359270068586 , 0.5097969509656486986060 , -0.3771134435310480315096 , 0.2011534205484152293408, - 0.4415200724898113993078, + 0.5454886553461798515130, - 0.2173248800239625522224, - 0.2171708216301463378883 , - 0.2460117034233289534662 , -0.3694743444265859433351 , 0.3601018742104445391483, 0.1530612374027510713681, + 0.6385835932752577104310, - 0.0013906755423099548263, - 0.3368505708673739662551 , 0.1394203608831885433 , -0.2292783363046769218308 , -0.2740379414191142504187, 0.1926279946047058377889 , + - 0.6106238336849327819067, 0.0128604884735820379493, 0.2248434453925002574071, - 0.0955950290268089419854 , 0.6322383898844411431739 , - 0.1362122761803175874373, - 0.0030545591332723984190, + 0.5431360047322275619308, 0.1304324282476526930541, - 0.3713242513274904177401, - 0.2693820083406764376299 , - 0.1431429883314006001882 , 0.0533021313694299267438, 0.0071247423818797811501 , + - 0.156705245229635081738, 0.0034676864415526253982 , - 0.3351413975178426096768 , 0.2927797449896782921996 , 0.5714932440406852443005 , 0.0738366305178744797288, - 0.4509415569247051669422 }; + double resultVectors4I[49]={ 0.0720991445669864616796 , 0.0844739028302603361942 , 0. , 0.0173431366502564965337, -0.0428710162141596462515, 0.0695226101295536302871, 0.0999385617869267273150, + 0.2510960793336660668018 , - 0.1123901435582568414384, 0.0091069807153507792430, 0. , - 0.2545493692188907641771, 0.1149714203529902251111, 0.0762796210842419941667, + 0. , 0.0611827470444509316505, 0.2295643426501567385678, - 0.2950878328821502361024, 0.198556748820033573955, 0.0646801132828138114483, - 0.2660572521172601678785, + 0. , 0.2688905463703794573860, 0.2230635356255223633593 , - 0.2013432448241845862391, - 0.0386922946770884290668, 0.0358963301038670995480, - 0.3593711366897305303780 , + - 0.1386717901043315326337, 0.0600994116374372472356, - 0.1167558596038882207102 , 0.317965599565130263571, 0. , -0.0642522754144960878131, - 0.0853665514099945371695, + 0. , 0.2715431731544318272320, - 0.1184277560110737481658 , 0.3655618467608299226868, - 0.4598766797780974302512 , 0.0188698345957667955319, - 0.1446429340244405892246, + 0.346095581161616094867 , 0.0750882555805542639682, - 0.2494867985551498246188, - 0.2328934592128872882455, 0. , -0.0047307828167728396829, - 0.0078254955461786041004}; + + + + + doubleComplex *in,*in2,*in3,*in4,out1[4],out2[4],out3[16],out4[16],out5[49],out6[49]; + + int i; + + in=DoubleComplexMatrix(inR,inI,4); + in2=DoubleComplexMatrix(in2R,in2I,4); + in3=DoubleComplexMatrix(in3R,in3I,16); + in4=DoubleComplexMatrix(in4R,in4I,49); + zspec2a(in,2,out1,out2); + + + for(i=0;i<4;i++){ + if (zreals(out1[i])>1e-16) assert( fabs(zreals(out1[i])-resultValuesR[i]) / fabs(zreals(out1[i])) <3e-16); + else assert(1); + if (zimags(out1[i])>1e-16) assert( fabs(zimags(out1[i])-resultValuesI[i]) / fabs(zimags(out1[i])) <1e-16); + else assert(1); + } + + for(i=0;i<4;i++){ + if (zreals(out2[i])>1e-16) assert( fabs(zreals(out2[i])-resultVectorsR[i]) / fabs(zreals(out2[i])) <3e-16); + else assert(1); + if (zimags(out2[i])>1e-16) assert( fabs(zimags(out2[i])-resultVectorsI[i]) / fabs(zimags(out2[i])) <1e-16); + else assert(1); + } + + + zspec2a(in2,2,out1,out2); + for(i=0;i<4;i++){ + if (zreals(out1[i])>1e-16) assert( fabs(zreals(out1[i])-resultValues2R[i]) / fabs(zreals(out1[i])) <3e-16); + else assert(1); + if (zimags(out1[i])>1e-16) assert( fabs(zimags(out1[i])-resultValues2I[i]) / fabs(zimags(out1[i])) <3e-15); + else assert(1); + } + + for(i=0;i<4;i++){ + if (zreals(out2[i])>1e-16) assert( fabs(zreals(out2[i])-resultVectors2R[i]) / fabs(zreals(out2[i])) <3e-16); + else assert(1); + if (zimags(out2[i])>1e-16) assert( fabs(zimags(out2[i])-resultVectors2I[i]) / fabs(zimags(out2[i])) <3e-15); + else assert(1); + } + + + + zspec2a(in3,4,out3,out4); + for(i=0;i<16;i++){ + if (zreals(out3[i])>1e-16) assert( fabs(zreals(out3[i])-resultValues3R[i]) / fabs(zreals(out3[i])) <3e-15); + else assert(1); + if (zimags(out3[i])>1e-16) assert( fabs(zimags(out3[i])-resultValues3I[i]) / fabs(zimags(out3[i])) <3e-15); + else assert(1); + } + for(i=0;i<16;i++){ + if (zreals(out4[i])>1e-16) assert( fabs(zreals(out4[i])-resultVectors3R[i]) / fabs(zreals(out4[i])) <3e-15); + else assert(1); + if (zimags(out4[i])>1e-16) assert( fabs(zimags(out4[i])-resultVectors3I[i]) / fabs(zimags(out4[i])) <3e-14); + else assert(1); + } + + /* FIXME : assert : 3e-16 maybe du to little values like 0.001... */ + zspec2a(in4,7,out5,out6); + + for(i=0;i<49;i++){ + if (zreals(out5[i])>1e-16) assert( fabs(zreals(out5[i])-resultValues4R[i]) / fabs(zreals(out5[i])) <3e-15); + else assert(1); + if (zimags(out5[i])>1e-16) assert( fabs(zimags(out5[i])-resultValues4I[i]) / fabs(zimags(out5[i])) <3e-14); + else assert(1); + } + for(i=0;i<49;i++){ + if (zreals(out6[i])>1e-16) assert( fabs(zreals(out6[i])-resultVectors4R[i]) / fabs(zreals(out6[i])) <3e-14); + else assert(1); + if (zimags(out6[i])>1e-16) assert( fabs(zimags(out6[i])-resultVectors4I[i]) / fabs(zimags(out6[i])) <3e-13); + else assert(1); + } + + + +} + + +static int testSpec2(void){ + printf(">>> Double Spec2 Tests <<<\n"); + printf(">>> Double <<<\n"); + dspec2aTest(); + printf(">>> DoubleComplex <<<\n"); + zspec2aTest(); + return 0; +} + + +int main(void){ + assert (testSpec2()==0); + + + return 0; +} diff --git a/2.3-1/src/c/linearAlgebra/spec2/testFloatSpec2.c b/2.3-1/src/c/linearAlgebra/spec2/testFloatSpec2.c new file mode 100644 index 00000000..fb2bc255 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/testFloatSpec2.c @@ -0,0 +1,134 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * 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 + * + */ + + +#include "spec.h" +#include "stdio.h" +#include "assert.h" +#include "math.h" + +static void sspec2aTest(void){ + float in[4]={1.0f,1.0f,1.0f,3.0f}; + float resultValuesR[4]={0.5857864376269050765700f,0,0,3.4142135623730949234300f}; + float resultVectorsR[4]={- 0.9238795325112867384831f,0.3826834323650897817792f, + 0.3826834323650897817792f,0.9238795325112867384831f}; + + float in2[4]={1.0f,1.0f,-2.0f,3.0f}; + float resultValues2R[4]={1.9999999999999997779554f,0,0,1.9999999999999997779554f}; + float resultVectors2R[4]={0.8164965809277261454824f,- 0.4082482904638631282523f, + 0.8164965809277261454824f,- 0.4082482904638631282523f}; + float eigenvalues[4],eigenvectors[4]; + + int i; + + sspec2a(in,2,eigenvalues,eigenvectors); + + for(i=0;i<4;i++){ + if (eigenvalues[i]>1e-6) assert( fabs(eigenvalues[i]-resultValuesR[i]) / fabs(eigenvalues[i]) <1e-16); + else assert(1); + } + for(i=0;i<4;i++){ + if (eigenvectors[i]>1e-6) assert( fabs(eigenvectors[i]-resultVectorsR[i]) / fabs(eigenvectors[i]) <1e-16); + else assert(1); + } + + + sspec2a(in2,2,eigenvalues,eigenvectors); + + + + for(i=0;i<4;i++){ + if (eigenvalues[i]>1e-6) assert( fabs(eigenvalues[i]-resultValues2R[i]) / fabs(eigenvalues[i]) <1e-16); + else assert(1); + } + for(i=0;i<4;i++){ + if (eigenvectors[i]>1e-6) assert( fabs(eigenvectors[i]-resultVectors2R[i]) / fabs(eigenvectors[i]) <1e-16); + else assert(1); + } +} + + + + +static void cspec2aTest(void){ + float inR[4]={1.0f,1.0f,1.0f,3.0f}; + float inI[4]={0.0f,0.0f,0.0f,0.0f}; + float resultValuesR[4]={0.5857864376269050765700f,0,0,3.4142135623730949234300f}; + float resultValuesI[4]={0,0,0,0}; + float resultVectorsR[4]={- 0.9238795325112867384831f,0.3826834323650897817792f, + 0.3826834323650897817792f,0.9238795325112867384831f}; + float resultVectorsI[4]={0,0,0,0}; + + float in2R[4]={1.0f,1.0f,-2.0f,3.0f}; + float in2I[4]={0.0f,0.0f,0.0f,0.0f}; + float resultValues2R[4]={1.9999999999999997779554f,0,0,1.9999999999999997779554f}; + float resultValues2I[4]={0.9999999999999997779554f,0,0,-0.9999999999999997779554f}; + float resultVectors2R[4]={0.8164965809277261454824f,- 0.4082482904638631282523f, + 0.8164965809277261454824f,- 0.4082482904638631282523f}; + float resultVectors2I[4]={0,- 0.4082482904638629062077f,0,0.4082482904638629062077f}; + + floatComplex *in,*in2,out1[4],out2[4]; + + int i; + + in=FloatComplexMatrix(inR,inI,4); + in2=FloatComplexMatrix(in2R,in2I,4); + + cspec2a(in,2,out1,out2); + for(i=0;i<4;i++){ + if (creals(out1[i])>1e-16) assert( fabs(creals(out1[i])-resultValuesR[i]) / fabs(creals(out1[i])) <1e-15); + else assert(1); + if (cimags(out1[i])>1e-16) assert( fabs(cimags(out1[i])-resultValuesI[i]) / fabs(cimags(out1[i])) <1e-16); + else assert(1); + } + for(i=0;i<4;i++){ + if (creals(out2[i])>1e-16) assert( fabs(creals(out2[i])-resultVectorsR[i]) / fabs(creals(out2[i])) <1e-15); + else assert(1); + if (cimags(out2[i])>1e-16) assert( fabs(cimags(out2[i])-resultVectorsI[i]) / fabs(cimags(out2[i])) <1e-16); + else assert(1); + } + + + cspec2a(in2,2,out1,out2); + for(i=0;i<4;i++){ + if (creals(out1[i])>1e-16) assert( fabs(creals(out1[i])-resultValues2R[i]) / fabs(creals(out1[i])) <1e-15); + else assert(1); + if (cimags(out1[i])>1e-16) assert( fabs(cimags(out1[i])-resultValues2I[i]) / fabs(cimags(out1[i])) <1e-15); + else assert(1); + } + for(i=0;i<4;i++){ + if (creals(out2[i])>1e-16) assert( fabs(creals(out2[i])-resultVectors2R[i]) / fabs(creals(out2[i])) <1e-15); + else assert(1); + if (cimags(out2[i])>1e-16) assert( fabs(cimags(out2[i])-resultVectors2I[i]) / fabs(cimags(out2[i])) <1e-16); + else assert(1); + } +} + + + + +static int testSpec2(void){ + printf(">>> Float Spec2 Tests <<<\n"); + printf(">>> Float <<<\n"); + sspec2aTest(); + printf(">>> FloatComplex <<<\n"); + cspec2aTest(); + return 0; +} + + +int main(void){ + assert (testSpec2()==0); + + + return 0; +} diff --git a/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj b/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj new file mode 100644 index 00000000..7c36943b --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj @@ -0,0 +1,178 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {FD0A4495-C372-41EF-B8EF-0A0939DA6C1E} + testDoubleSpec2 + Win32Proj + + + + Application + Unicode + true + + + Application + Unicode + true + + + Application + Unicode + + + Application + Unicode + + + + + + + + + + + + + + + + + + + <_ProjectFileVersion>10.0.40219.1 + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + true + true + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + false + false + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + EditAndContinue + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + MachineX86 + + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + MachineX86 + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + + + + + CompileAsCpp + CompileAsCpp + CompileAsCpp + CompileAsCpp + + + + + {9b1bd750-1fef-4d6b-9422-782d16181cee} + false + + + {fd335544-52bf-4736-a34e-77f591d158d5} + false + + + + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj.filters b/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj.filters new file mode 100644 index 00000000..04f409dc --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/test_DoubleSpec2/testDoubleSpec2.vcxproj.filters @@ -0,0 +1,22 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj b/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj new file mode 100644 index 00000000..f4154eb1 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj @@ -0,0 +1,178 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {E20D0BC2-B28A-412C-9715-5B93B319A4A3} + testFloatSpec2 + Win32Proj + + + + Application + Unicode + true + + + Application + Unicode + true + + + Application + Unicode + + + Application + Unicode + + + + + + + + + + + + + + + + + + + <_ProjectFileVersion>10.0.40219.1 + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + true + true + $(SolutionDir)bin\ + $(SolutionDir)bin\ + $(ProjectDir)$(Configuration)\ + $(ProjectDir)$(Configuration)\ + false + false + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + EditAndContinue + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + MachineX86 + + + + + Disabled + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + EnableFastChecks + MultiThreadedDebugDLL + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + MachineX86 + + + + + MaxSpeed + true + ../../../includes;../../../type;../../includes;../../../operations/includes;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + true + + + Level3 + ProgramDatabase + + + $(SolutionDir)bin\$(ProjectName).exe + true + Console + true + true + + + + + CompileAsCpp + CompileAsCpp + CompileAsCpp + CompileAsCpp + + + + + {9b1bd750-1fef-4d6b-9422-782d16181cee} + false + + + {fd335544-52bf-4736-a34e-77f591d158d5} + false + + + + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj.filters b/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj.filters new file mode 100644 index 00000000..0c38d2f2 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/test_FloatSpec2/testFloatSpec2.vcxproj.filters @@ -0,0 +1,22 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + \ No newline at end of file diff --git a/2.3-1/src/c/linearAlgebra/spec2/zspec2a.c b/2.3-1/src/c/linearAlgebra/spec2/zspec2a.c new file mode 100644 index 00000000..e33ed897 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/spec2/zspec2a.c @@ -0,0 +1,96 @@ +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008 - INRIA - Arnaud TORSET + * + * 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 + * + */ + +#include +#include "spec.h" +#include "lapack.h" +#include "zeros.h" +#include "stdio.h" +#include "conj.h" + +void zspec2a(doubleComplex* in, int rows,doubleComplex* eigenvalues, doubleComplex* eigenvectors){ + int i = 0, j = 0; + int hermitian = 0; + int INFO = 0; + int iWorkSize = 0; + doubleComplex* pdblWork; + doubleComplex* pdblLeftvectors; + doubleComplex* pdblRightvectors; + doubleComplex* pdblRWork; + double* pdblRWork2; + double* outReal; + double* outImag; + doubleComplex* inCopy; + + inCopy = (doubleComplex*)malloc((unsigned int)(rows*rows) * sizeof(doubleComplex)); + outReal = (double*)malloc((unsigned int)rows * sizeof(double)); + outImag = (double*)malloc((unsigned int)rows * sizeof(double)); + pdblLeftvectors=NULL; + pdblRightvectors=NULL; + + iWorkSize = 2*rows; + pdblWork = (doubleComplex*)malloc((unsigned int)iWorkSize * sizeof(doubleComplex)); + pdblRWork = NULL; + pdblRWork2 = NULL; + + for(i=0;i #include "string.h" #include "stdlib.h" #include "lapack.h" #include "matrixTranspose.h" #include "matrixMultiplication.h" +#include "doubleComplex.h" -void dspec1a(double *in1,double *in2,int size,int nout,double *out1,double *out2,double *out3,double *out4){ +extern int dggev_(char *,char *,int *,double *,int *,double *,int *,double *,double *,double *,double *,int *,double *,int *,double *,int *,int *); + +void assembleEigenvectorsInPlace(int N,double *ALPHAI,double *EVreal,double *EVimg){ + int j,i; + j=0; + while(j