diff options
Diffstat (limited to '2.3-1/src/c/linearAlgebra')
23 files changed, 1290 insertions, 20 deletions
diff --git a/2.3-1/src/c/linearAlgebra/fullrf/dfullrfa.c b/2.3-1/src/c/linearAlgebra/fullrf/dfullrfa.c new file mode 100644 index 00000000..a409ae35 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/fullrf/dfullrfa.c @@ -0,0 +1,112 @@ +/* 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 + */ + +/* FULL Rank factorization function in scilab */ + +/* + //[Q,M,rk]=fullrf(A) + //Full rank factorization : A=Q.M + //with range(Q)=range(A) and ker(M)=ker(A), + //Q full column rank , M full row rank + // rk = rank(A) = #columns(Q) = #rows(M) + //F.D. +*/ + +#include "fullrf.h" +#include <stdio.h> +#include <stdlib.h> +#include "svd.h" +#include <math.h> +#include "norm.h" +#include "matrixTranspose.h" +#include "matrixMultiplication.h" + +double dfullrfa(int ninp,double *inp1,int row,int col,double tol,double *out1,double *out2){ + + int i,j; + /* norm inp1 - norm(inp1,1)*/ + double na1; + na1 = dnorma(inp1,row,col,1); + + if(ninp == 1){ + tol = sqrt(pow(2,-52)); + } + + if(na1 < pow(1,-10)){ + out1 = NULL; + out2 = NULL; + return 0; + } + + double tol1; + tol1 = tol*na1; + + double *U,*S,*V; + U = (double *)malloc(row*row*sizeof(double)); + S = (double *)malloc(row*col*sizeof(double)); + V = (double *)malloc(col*col*sizeof(double)); + + double rk; + rk = dsvda(tol1,inp1,row,col,0,4,U,S,V); + + /* sq = sqrt(s) */ + for(i=0;i<row;i++){ + for(j=0;j<col;j++){ + if(i == j){ + S[i*row+j] = pow(S[i*row+j],0.5); + } + else{ + S[i*row+j] = 0; + } + } + } + + double *Q; + Q = (double *)malloc(row*col*sizeof(double)); + dmulma(U,row,row,S,row,col,Q); + + double *VT; + VT = (double *)malloc(col*col*sizeof(double)); + dtransposea(V,col,col,VT); + + /* multiplication of sq*V' or S*VT */ + double *M; + M = (double *)malloc(row*col*sizeof(double)); + dmulma(S,row,col,VT,col,col,M); + + /* This Program is not yet completed properly, as it outputs the whole matrix, instead of the exact output. + + if anyone finds, how to fix the size in INITFILLscidir.sci +please change there and change below few lines of codes accordingly. + +*/ + for(i=0;i<row;i++){ + for(j=0;j<col;j++){ + //if(j < rk) + out1[i*col+j] = Q[i*col+j]; + //else + // out1[i*col+j] = 0; + } + //printf("\n"); + } + + for(i=0;i<row;i++){ + for(j=0;j<col;j++){ + //if(i < rk) + out2[i*col+j] = M[i*col+j]; + //else + // out2[i*col+j] = 0; + } + } + + return rk; +} diff --git a/2.3-1/src/c/linearAlgebra/givens/dgivensa.c b/2.3-1/src/c/linearAlgebra/givens/dgivensa.c new file mode 100644 index 00000000..9bf0637b --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/givens/dgivensa.c @@ -0,0 +1,76 @@ +/* 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 + */ + +/* GIVENS scilab function + Syntax : u=givens(xy) + u=givens(x,y) + xy = [x;y], u=givens(xy) + returns a 2*2 matrix u such that u*xy=[r;0]. + c is equal to u*xy + givens(x,y)=givens([x;y]) +*/ + +#include "givens.h" +#include <stdio.h> +#include <math.h> +#include "norm.h" + +/* All variable names are according to scilab code */ + +void dgivensa(int ninp,double *inp1,int row,int col,double *inp2,int row1,int col1,int nout,double *out1,double *out2){ + double *x; + double r; + x = (double *)malloc((double)2*sizeof(double)); + if(ninp == 2){ + if(row != 1 || col != 1 || row1 != 1 || col1 != 1){ + printf("Wrong size for input argument: A column vector expected.\n"); + return; + } + *(x) = *(inp1); + *(x+1) = *(inp2); + } + else{ + //printf("(%d %d)",row,col); + if(row != 2 || col != 1){ + printf("Wrong size for input argument: A column vector expected.\n"); + return; + } + //printf("(%lf %lf)",inp1[0],inp1[1]); + x[0] = inp1[0]; + x[1] = inp1[1]; + } + if(*(x+1) != 0){ + //printf("(%lf %lf)",x[0],x[1]); + /*Norm of type 2 - find the maximum singular value*/ + r = dnorma(x,2,1,2); + //printf("%lf \n",r); + *(out1) = (*(x))/r; + *(out1+1) = -(*(x+1))/r; + *(out1+2) = (*(x+1))/r; + *(out1+3) = (*(x))/r; + if(nout == 2){ + *(out2) = r; + *(out2+1) = 0; + } + + } + else{ + *(out1) = 1; + *(out1+1) = 0; + *(out1+2) = 1; + *(out1+3) = 0; + if(nout == 2){ + *(out2) = *(x); + *(out2+1) = *(x+1); + } + } +} diff --git a/2.3-1/src/c/linearAlgebra/householder/dhouseholdera.c b/2.3-1/src/c/linearAlgebra/householder/dhouseholdera.c new file mode 100644 index 00000000..5a98bfae --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/householder/dhouseholdera.c @@ -0,0 +1,90 @@ +/* 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 + */ + +/* Householder orthogonal reflexion matrix */ + +/* +Syntax :- + //u=householder(v [,w]) + //Description + //given 2 column vectors v w of same size householder(v,w) returns a unitary + //column vector u, such that (eye-2*u*u')*v is proportional to w. + //(eye-2*u*u') is the orthogonal Householder reflexion matrix + // + // w default value is eye(v). In this case vector (eye-2*u*u')*v is the + // vector eye(v)*(+-norm(v)) +*/ +#include <stdio.h> +#include <stdlib.h> +#include "householder.h" +#include "eye.h" +#include "matrixTranspose.h" +#include "matrixMultiplication.h" +#include <math.h> +#include "norm.h" + +void dhouseholdera(int ninp,double *inp1,int row,double *inp2,double *out1){ + + int i,j; + double *x; + x = (double *)malloc(row*sizeof(double)); + + for(i=0;i<row;i++) + x[i] = inp1[i]; + + if(ninp < 2){ + deyea(x,row,1); + } + else{ + for(i=0;i<row;i++){ + x[i] = inp2[i]; + } + } + + /* vt transpose of inp1 */ + + double *vt; + vt = (double *)malloc(row*sizeof(double)); + dtransposea(inp1,row,1,vt); + + /*wt transpose of inp2 */ + double *wt; + wt = (double *)malloc(row*sizeof(double)); + dtransposea(x,row,1,wt); + + /* vvt = inp1*vt */ + double *vvt; + vvt = (double *)malloc(1*1*sizeof(double)); + dmulma(vt,1,row,inp1,row,1,vvt); + + /* wwt = inp2*wt */ + double *wwt; + wwt = (double *)malloc(1*1*sizeof(double)); + dmulma(wt,1,row,x,row,1,wwt); + + /* a=-sqrt((v'*v)/(w'*w)) */ + double a; + //a = (double *)malloc(1*1*sizeof(double)); + a = -sqrt(vvt[0]/wwt[0]); + + for(i=0;i<row;i++){ + out1[i] = x[i]*a+inp1[i]; + } + + /* norm of out1 */ + double r; + r = dnorma(out1,row,1,2); + + for(i=0;i<row;i++){ + out1[i]=out1[i]/r; + } +} diff --git a/2.3-1/src/c/linearAlgebra/includes/fullrf.h b/2.3-1/src/c/linearAlgebra/includes/fullrf.h new file mode 100644 index 00000000..cc0a33d0 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/fullrf.h @@ -0,0 +1,26 @@ + /* 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 __FULLRF_H__ +#define __FULLRF_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +double dfullrfa(int ninp,double *in1,int row,int col,double tol,double *out1,double *out2); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__FULLRF_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/includes/givens.h b/2.3-1/src/c/linearAlgebra/includes/givens.h new file mode 100644 index 00000000..4aac91b2 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/givens.h @@ -0,0 +1,25 @@ + /* 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 __GIVENS_H__ +#define __GIVENS_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +void dgivensa(int ninp,double *inp1,int row,int col,double *inp2,int row1,int col1,int nout,double *out1,double *out2); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__givens_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/includes/householder.h b/2.3-1/src/c/linearAlgebra/includes/householder.h new file mode 100644 index 00000000..64350a15 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/householder.h @@ -0,0 +1,26 @@ + /* 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 __HOUSEHOLDER_H__ +#define __HOUSEHOLDER_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void dhouseholdera(int ninp,double *inp1,int row,double *inp2,double *out1); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__HOUSEHOLDER_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/includes/qr.h b/2.3-1/src/c/linearAlgebra/includes/qr.h new file mode 100644 index 00000000..2ed12e3a --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/qr.h @@ -0,0 +1,26 @@ + /* 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 __QR_H__ +#define __QR_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +double dqra(int ninp,int nout,double *inp1,int M,int N,double tol,double *out1,double *out2,double *out3); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__QR_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/includes/rowcomp.h b/2.3-1/src/c/linearAlgebra/includes/rowcomp.h new file mode 100644 index 00000000..faf5a2a7 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/rowcomp.h @@ -0,0 +1,26 @@ + /* 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 __ROWCOMP_H__ +#define __ROWCOMP_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +double drowcompa(int ninp,double *A,int row,int col,char *flag,double tol,double *w); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__ROWCOMP_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/includes/sqroot.h b/2.3-1/src/c/linearAlgebra/includes/sqroot.h new file mode 100644 index 00000000..9c1d9652 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/includes/sqroot.h @@ -0,0 +1,26 @@ + /* 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 __SQROOT_H__ +#define __SQROOT_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void dsqroota(double *inp,int row,int col,double *out); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__SQROOT_H__*/ diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_fullrf.h b/2.3-1/src/c/linearAlgebra/interfaces/int_fullrf.h new file mode 100644 index 00000000..1b8a067b --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_fullrf.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_FULLRF_H__ +#define __INT_FULLRF_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2fullrfd2d2d0(in1,size,out1,out2) dfullrfa(1,in1,size[0],size[1],0,out1,out2); +#define d2d0fullrfd2d2d0(in1,size,in2,out1,out2) dfullrfa(2,in1,size[0],size[1],in2,out1,out2); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_FULLRF_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_givens.h b/2.3-1/src/c/linearAlgebra/interfaces/int_givens.h new file mode 100644 index 00000000..ba30dbce --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_givens.h @@ -0,0 +1,32 @@ + /* 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_GIVENS_H__ +#define __INT_GIVENS_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2givensd2(in1,size,out1) dgivensa(1,in1,size[0],size[1],NULL,0,0,1,out1,NULL); +#define d2givensd2d2(in1,size,out1,out2) dgivensa(1,in1,size[0],size[1],NULL,0,0,2,out1,out2); + +#define d0d0givensd2d2(in1,in2,out1,out2) dgivensa(2,&in1,1,1,&in2,1,1,2,out1,out2); +#define d0d0givensd2(in1,out1) dgivensa(2,&in1,1,1,&in2,1,1,1,out1,NULL); + + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_GIVENS_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_householder.h b/2.3-1/src/c/linearAlgebra/interfaces/int_householder.h new file mode 100644 index 00000000..f8637197 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_householder.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_HOUSEHOLDER_H__ +#define __INT_HOUSEHOLDER_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2householderd2(in1,size,out1) dhouseholdera(1,in1,size[0],NULL,out2); +#define d2d2householderd2(in1,size1,in2,size2,out1) dhouseholdera(2,in1,size1[0],in2,out1); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_HOUSEHOLDER_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_qr.h b/2.3-1/src/c/linearAlgebra/interfaces/int_qr.h new file mode 100644 index 00000000..d34d8f41 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_qr.h @@ -0,0 +1,34 @@ + /* 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_QR_H__ +#define __INT_QR_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2qrd2d2(in1,size,out1,out2) dqra(1,2,in1,size[0],size[1],0,out1,out2,NULL); +#define d2qrd2d2d2(in1,size,out1,out2,out3) dqra(1,3,in1,size[0],size[1],0,out1,out2,out3); +#define d2g2qrd2d2(in1,size,in2,size1,out1,out2) dqra(2,2,in1,size[0],size[1],0,out1,out2,NULL); +#define d2g2qrd2d2d2(in1,size,in2,size1,out1,out2,out3) dqra(2,3,in1,size[0],size[1],0,out1,out2,out3); +#define d2d0qrd2d2d0d2(in1,size,in2,out1,out2,out3) dqra(2,4,in1,size[0],size[1],in2,out1,out2,out3); +#define d2qrd2d2d0d2(in1,size,out1,out2,out3) dqra(1,4,in1,size[0],size[1],0,out1,out2,out3); + + + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_HESS_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_rowcomp.h b/2.3-1/src/c/linearAlgebra/interfaces/int_rowcomp.h new file mode 100644 index 00000000..b72687d5 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_rowcomp.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 __INT_ROWCOMP_H__ +#define __INT_ROWCOMP_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2rowcompd2d0(in1,size,out1) drowcompa(1,in1,size[0],size[1],NULL,0,out1); +#define d2g2rowcompd2d0(in1,size,flag,size1,out1) drowcompa(2,in1,size[0],size[1],flag,0,out1); +#define d2g2d0rowcompd2d0(in1,size,flag,size1,tol,out1) drowcompa(2,in1,size[0],size[1],flag,tol,out1); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_ROWCOMP_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_sqroot.h b/2.3-1/src/c/linearAlgebra/interfaces/int_sqroot.h new file mode 100644 index 00000000..57af2c08 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_sqroot.h @@ -0,0 +1,27 @@ + /* 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_SQROOT_H__ +#define __INT_SQROOT_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +#define d2sqrootd2(inp,size,out) dsqroota(inp,size[0],size[1],out); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__INT_HESS_H__*/ + diff --git a/2.3-1/src/c/linearAlgebra/proj/dproja.c b/2.3-1/src/c/linearAlgebra/proj/dproja.c new file mode 100644 index 00000000..e27cd6f2 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/proj/dproja.c @@ -0,0 +1,73 @@ +/* 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 + */ + +/* Scilab function proj code in C */ + +#include <stdio.h> +#include <stdlib.h> +#include "proj.h" +#include "eye.h" +#include "matrixTranspose.h" +#include "matrixMultiplication.h" +#include <math.h> +#include "norm.h" + + +double dproja(double *x1,int l,int k,double *x2,int m2,int n2,double *y){ + int i,j; + double *w; + w = (double *)malloc(l*l*sizeof(double)); + double rk; + rk = drowcompa(1,x1,l,k,NULL,0,w); + + double *w1; + w1 = (double *)malloc(rk*l*sizeof(double)); + + for(i=0;i<l*rk;i++){ + w1[i]=w[i]; + } + + double *x1t; + x1t = (double *)malloc((l-n)*l*sizeof(double); + + for(i=n;i<l;i++){ + for(j=0;j<l;j++){ + x1t[i-n+j*l] = w1[i+j*l]; + } + } + + double x1x2; + x1x2 = (double *)malloc((l-n+1)*n2*sizeof(double)); + dmulma(x1t,l-n+1,l,x2,m2,n2,x1x2); + + double *inx1x2; + inx1x2 = (double *)malloc(); + + + + + + + + + + + + + + + + + + + +} diff --git a/2.3-1/src/c/linearAlgebra/projspec/dprojspeca.c b/2.3-1/src/c/linearAlgebra/projspec/dprojspeca.c new file mode 100644 index 00000000..6bb6fa52 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/projspec/dprojspeca.c @@ -0,0 +1,70 @@ +/* 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 + */ + +/* PROJSPEC function in scilab */ + +/* + //[S,P,D,index]=projspec(A) + //Spectral characteristics of A at 0 + //S = reduced resolvent at 0 (S=-Drazin_inverse(A)) + //P = spectral projection at 0 + //D = Nilpotent operator at 0 + //index = index of the 0 eigenvalue + //! + +*/ + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include "norm.h" +#include "eye.h" + +#define eps pow(2,-52) + +double dprojspeca(double *inp1,int row,int col,double *out1,double *out2,double *out3){ + double tol = pow(10,-6); + int i,j,index; + /*norm(A,1)*/ + double nor; + nor = dnorma(inp1,row,col,1); + + /* P=eye(A),D=A,S=0*P;index=1; */ + if(nor < eps*row*row){ + memcpy(out2,inp1,row*col*sizeof(double)); + deyea(out2,row,col); + memcpy(out3,inp1,row*col*sizeof(double)); + for(i=0;i<row;i++){ + for(j=0;j<row;j++){ + out1[i*row+j] = 0; + } + } + index = 1; + } + + /* rcond(A) */ + double *rcon; + rcon = rcond(inp1,row); + if(rcon > tol){ + dinverma(inp1,out1,row); + for(i=0;i<row*col;i++){ + out2[i]=0; + out3[i]=0; + } + index = 0; + return index; + } + + index = 1; + + +} diff --git a/2.3-1/src/c/linearAlgebra/qr/dqra.c b/2.3-1/src/c/linearAlgebra/qr/dqra.c new file mode 100644 index 00000000..bae4bc27 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/qr/dqra.c @@ -0,0 +1,298 @@ +/* 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 C code is used to generate function for QR decomposition */ + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include "qr.h" +#include "lapack.h" +#include "string.h" +#include "matrixTranspose.h" + +/*For reference check Scilab source code & lapack library websites +Names of variable are almost same for convience. +*/ +/* +int min(int M,int N){ + if(M > N) return N; + return M; +} + +int max(int M,int N){ + if(M > N) return N; + return N; +}*/ + +/* External Function used of lapack library */ +extern double dgeqrf_(int *,int *,double *,int *,double *,double *,int *,int *); +extern double dlacpy_(char *,int *,int *,double *,int *,double *,int *); +extern double dgeqpf_(int *,int *,double *,int *,int *,double *,double *,int *); +extern double dorgqr_(int *,int *,int *,double *,int *,double *,double *,int *,int *); +extern void dlaset_(char *,int *,int *,double *,double *,double *,int *); + +/* function for finding qr */ +double dqra(int ninp,int nout,double *inp1,int M,int N,double tol,double *out1,double *out2,double *out3){ + int i,j; + char choice; + double alpha=0.0,beta=0.0; + int minMN = min(M,N); + + double *A; + A = (double *)malloc(M*N*sizeof(double)); + memcpy(A,inp1,M*N*sizeof(double)); + + if(M <= 0 || N <= 0){ + out1 = NULL; + out2 = NULL; + return 0; + } + /* doldqr */ + if(nout == 4){ /* [Q,R,rk,E]=qr(X [,tol]) */ + if(ninp == 1){ + tol = -1; + } + + int INFO,rk; + + int *JPVT; + JPVT = ( int *)malloc(N*sizeof(int)); + + double *TAU; + TAU = (double *)malloc(minMN*sizeof(double)); + + int LWORK = 3*N; + + double *WORK; + WORK = (double *)malloc(LWORK*sizeof(double)); + + for(i=1;i<=M;i++){ + JPVT[i-1]=0.0; + } + dgeqpf_(&M,&N,A,&M,JPVT,TAU,WORK,&INFO); + + choice = 'U'; + dlacpy_(&choice,&M,&N,A,&M,out2,&M); + + if(M > N){ + for(j=1;j<=N;j++){ + for(i=j+1;i<=M;i++){ + out2[i-1+(j-1)*M] = 0.0; + } + } + } + else{ + for(j=1;j<=M-1;j++){ + for(i=j+1;i<=M;i++){ + out2[i-1+(j-1)*M] = 0.0; + } + } + } + + if(M > N){ + choice = 'F'; + dlacpy_(&choice,&M,&N,A,&M,out1,&M); + for(j=N+1;j<=M;j++){ + for(i=1;i<=M;i++){ + out1[i-1+(j-1)*M] = 0.0; + } + } + } + else{ + choice = 'F'; + dlacpy_(&choice,&M,&M,A,&M,out1,&M); + } + + dorgqr_(&M,&M,&minMN,out1,&M,TAU,WORK,&LWORK,&INFO); + + choice = 'F'; + dlaset_(&choice,&N,&N,&alpha,&beta,out3,&N); + + for(j=1;j<=N;j++){ + i = JPVT[j-1]; + out3[i-1+(j-1)*N] = 1.0; + } + + double tt = abs(out2[0]); + + if(tol == -1){ + tol = (double)max(M,N)*pow(2,-52)*tt; + } + double ch; + //printf("%lf ",tol); + for(j=1;j<=minMN;j++){ + //printf("%e ",out2[(j-1)+(j-1)*M]); + if(out2[(j-1)+(j-1)*M] < 0){ + ch = -out2[(j-1)+(j-1)*M]; + } + else{ + ch = out2[(j-1)+(j-1)*M]; + } + if(ch > tol){ + rk = j; + } + else{ + break; + } + } + //printf("\n"); + return rk; + } + else if(ninp == 1){ + /* (intdgeqpf3) + [Q,R]=qr(A) + [Q,R,E]=qr(A) + */ + int LDA = M; + + double *TAU; + TAU = (double *)malloc(min(M,N)*sizeof(double)); + + int LWORK; + if(nout <= 2) + LWORK = N; + else + LWORK = 3*N; + + double *WORK; + WORK = (double *)malloc((LWORK+1)*sizeof(double)); + + int INFO; + + int *JPVT; + JPVT = (int *)malloc(N*sizeof(int)); + + if(nout <= 2){ + dgeqrf_(&M,&N,A,&M,TAU,WORK,&LWORK,&INFO); + } + else{ + for(i=0;i<N;i++){ + *(JPVT+i)=0.0; + } + dgeqpf_(&M,&N,A,&M,JPVT,TAU,WORK,&INFO); + } + + /* Copying code from A to R */ + choice = 'U'; + dlacpy_(&choice,&M,&N,A,&M,out2,&M); + /*for(i=0;i<M;i++){ + for(j=0;j<N-i;j++){ + out2[i+j*M] = A[i+j*M]; + } + }*/ + + if(M > N){ + for(j=1;j<=N;j++){ + for(i=j+1;i<=M;i++){ + out2[i-1+(j-1)*M] = 0.0; + } + } + } + else{ + for(j=1;j<=M-1;j++){ + for(i=j+1;i<=M;i++){ + out2[i-1+(j-1)*M] = 0.0; + } + } + } + /* lQ - out1 */ + if(M > N){ + choice = 'F'; + dlacpy_(&choice,&M,&N,A,&M,out1,&M); + /*for(i=0;i<M*N;i++){ + out1[i]=A[i]; + }*/ + for(j=N+1;j<=M;j++){ + for(i=1;i<=M;i++){ + out1[i+(j-1)*M] = 0.0; + } + } + } + else{ + choice = 'F'; + dlacpy_(&choice,&M,&M,A,&M,out1,&M); + /*for(i=0;i<M;i++){ + for(j=0;j<M;j++){ + out1 + } + }*/ + } + dorgqr_(&M,&M,&minMN,out1,&M,TAU,WORK,&LWORK,&INFO); + + if(nout > 2){ + choice = 'F'; + dlaset_(&choice,&N,&N,&alpha,&beta,out3,&N); + for(j=1;j<=N;j++){ + i = *(JPVT+j-1); + //printf("%d ",i-1+(j-1)*N); + *(out3+i-1+(j-1)*N) = 1.0; + } + //printf("\n"); + } + } + else{/* + [[Q,R]=qr(A,'e') + [Q,R,E]=qr(A,'e') ] */ + int *JPVT; + JPVT = (int *)malloc(N*sizeof(int)); + + double *TAU; + TAU = (double *)malloc(min(M,N)*sizeof(double)); + int LWORK; + if(nout <= 2){ + LWORK = N; + } + else{ + LWORK = 3*N; + } + + double *WORK; + WORK = (double *)malloc(LWORK*sizeof(double)); + + int INFO; + if(nout <= 2) + dgeqrf_(&M,&N,A,&M,TAU,WORK,&LWORK,&INFO); + else{ + for(i=1;i<=N;i++){ + JPVT[i-1]=0.0; + } + dgeqpf_(&M,&N,A,&M,JPVT,TAU,WORK,&INFO); + } + choice = 'U'; + + dlacpy_(&choice,&minMN,&N,A,&M,out2,&minMN); + + if(N >= 2){ + for(j=1;j<=N-1;j++){ + if(j+1 <= minMN){ + for(i=j+1;i<=minMN;i++){ + out2[i-1+(j-1)*minMN] = 0.0; + } + } + } + } + choice = 'F'; + dlacpy_(&choice,&M,&minMN,A,&M,out1,&M); + dorgqr_(&M,&minMN,&minMN,out1,&M,TAU,WORK,&LWORK,&INFO); + + if(nout > 2){ + choice = 'F'; + dlaset_(&choice,&N,&N,&alpha,&beta,out3,&N); + for(j=1;j<=N;j++){ + i = JPVT[j-1]; + out3[i-1+(j-1)*N] = 1.0; + } + } + } + return 0; +} diff --git a/2.3-1/src/c/linearAlgebra/rowcomp/drowcompa.c b/2.3-1/src/c/linearAlgebra/rowcomp/drowcompa.c new file mode 100644 index 00000000..3161a2d6 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/rowcomp/drowcompa.c @@ -0,0 +1,79 @@ +/* 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 is used to find row compression, range */ + +#include "rowcomp.h" +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include "svd.h" +#include "norm.h" +#include "eye.h" +#include "matrixTranspose.h" +#include "qr.h" + +/* All variable names, are in consideration of scilab documentation. for reference please check the scilab code.*/ + +double drowcompa(int ninp,double *A,int row,int col,char *flag,double tol,double *w){ + double rk; + double *U; + double *S; + double *V; + double *q,*r,*e; + if(row == 0 || col == 0){ + w = NULL; + return 0; + } + + double nA1 = dnorma(A,row,col,1); + if(nA1 < sqrt(pow(2,-52))/10){ + deyea(w,row,col); + return 0; + } + if(ninp == 1){ + flag = "svd"; + tol = sqrt(pow(2,-52))*nA1; + } + else if(ninp == 2){ + tol = sqrt(pow(2,-52))*nA1; + } + else{ + if(tol < 0){ /* if tolerance is negative */ + printf(" Wrong values for input argument #: Non-negative scalar expected"); + } + } + int M = row,N=col; + int minMN = min(M,N); + char check[3]="qr"; + //printf(" %s ",flag); + if(strcmp(check,flag) == 0){ + /* calling qr function*/ + //printf(" * "); + q = (double *)malloc(M*min(M,N)*sizeof(double)); + r = (double *)malloc(minMN*N*sizeof(double)); + e = (double *)malloc(N*N*sizeof(double)); + rk = dqra(2,4,A,M,N,tol,q,r,e); + memcpy(w,q,row*col*sizeof(double)); + dtransposea(q,row,row,w); + return rk; + } + else{ + /* svd function type */ + U = (double *)malloc(row*row*sizeof(double)); + S = (double *)malloc(row*col*sizeof(double)); + V = (double *)malloc(col*col*sizeof(double)); + rk = dsvda(tol,A,row,col,0,4,U,S,V); + dtransposea(U,row,row,w); + return rk; + } +} diff --git a/2.3-1/src/c/linearAlgebra/sqroot/dsqroota.c b/2.3-1/src/c/linearAlgebra/sqroot/dsqroota.c new file mode 100644 index 00000000..038afc98 --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/sqroot/dsqroota.c @@ -0,0 +1,130 @@ +/* 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 + */ + +/* Function - sqroot of scilab, W*W' hermitian factorization */ + +#include <stdio.h> +#include "stdlib.h" +#include "string.h" +#include <math.h> +#include "matrixTranspose.h" +#include "svd.h" +#include "matrixMultiplication.h" + +#define eps 2.22044604925e-16 + +/*It would be good, if you are doing conversoins for only valid inputs before + running the program, check all the conditions before hand. +*/ + +void dsqroota(double *inp,int row,int col, double *out){ + if(row != col){ + printf("Enter valid inputs only - matrix should be symetric\n"); + return; + } + int i,j; + double *U,*S,*V; + double *A,*B; + int rk; + + U = (double *)malloc((double)row*row*sizeof(double)); + S = (double *)malloc((double)Min(row,col)*Min(row,col)*sizeof(double)); + V = (double *)malloc((double)col*col*sizeof(double)); + A = (double *)malloc(rk*rk*sizeof(double)); + B = (double *)malloc(rk*row*sizeof(double)); + + double *Q1; + Q1 = (double *)malloc(row*col*sizeof(double)); + + dtransposea(inp,row,col,Q1); + + double *Q2; + Q2 = (double *)malloc(col*row*sizeof(double)); + + /* Q2 = (inp+inp1')/2; */ + for(i=0;i<row;i++){ + for(j=0;j<row;j++){ + Q2[i+j*row] = ((inp[i+j*row]+Q1[i+j*row])/2); + } + } + + /* norm(Q2-Q1,1) - finding the max value from sum of columns */ + double sum = 0; + double maxi=0; + for(i=0;i<row;i++){ + sum = 0; + for(j=0;j<col;j++) sum+=(Q2[i*row+j]-inp[i*row+j]); + if(maxi < sum){ + maxi = sum; + } + } + + /* if norm(Q1-Q,1) > 100*%eps then */ + if(maxi > 100*eps){ + printf("Warning: Wrong size for input argument and Symmetric expected\n"); + } + maxi = 0; + for(i=0;i<col;i++){ + sum=0; + for(j=0;j<row;j++) sum+=Q1[i*row+j]; + if(maxi < sum){ + maxi = sum; + } + } + + /*if norm(Q,1) < sqrt(%eps) then S=[];return;end*/ + + if(maxi < sqrt(eps)){ + out = NULL; + } + else{ + rk = dsvda(0,inp,row,col,0,4,U,S,V); + + /*Will be used in complex numbers*/ + //C = (double *)malloc(rk*row*sizeof(double)); + + /*for(i=0;i<row;i++){ + for(j=0;j<row;j++){ + printf("%lf ",S[i*row+j]); + } + printf("\n"); + }*/ + /*sqrt of S*/ + printf("%d ",rk); + for(i=0;i<rk;i++){ + for(j=0;j<rk;j++){ + A[i*rk+j] = sqrt(S[i*row+j]); + //printf("%lf ",A[i*rk+j]); + } + } + for(i=0;i<col*rk;i++){ + B[i] = V[i]; + printf("%lf ",B[i]); + } + printf("\n"); + /*for(i=0;i<rk;i++){ + for(j=0;j<rk;j++){ + //A[i*rk+j] = sqrt(S[i*row+j]); + printf("%lf ",A[i*rk+j]); + } + printf("\n"); + }*/ + /*for(i=0;i<col;i++){ + for(j=0;j<rk;j++){ + //B[i*col+j] = VT[i*col+j]; + printf("%lf ",B[i*col+j]); + } + printf("\n"); + }*/ + dmulma(B,col,rk,A,rk,rk,out); + } +} diff --git a/2.3-1/src/c/linearAlgebra/sva/dsvaa.c b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c index b7d07d8c..0d34d8d1 100644 --- a/2.3-1/src/c/linearAlgebra/sva/dsvaa.c +++ b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c @@ -33,14 +33,14 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ /* 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)); + 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; + tol = Max(row,col)*S[0]*eps; rk = 0; for(i=0;i<col;i++){ if(S[i+i*row] > tol){ @@ -52,7 +52,7 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ tol = in2; if(tol > 1){ rk = tol; - if(rk > min(row,col)){ + if(rk > Min(row,col)){ printf("ERROR: Wrong value for input argument !"); out1 = NULL; out2 = NULL; @@ -70,20 +70,20 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ } } arow = M; - acol = min(M,N); + 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); + 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); + 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/.1.c.swp b/2.3-1/src/c/linearAlgebra/svd/.1.c.swp Binary files differnew file mode 100644 index 00000000..81d9e9cf --- /dev/null +++ b/2.3-1/src/c/linearAlgebra/svd/.1.c.swp diff --git a/2.3-1/src/c/linearAlgebra/svd/zsvda.c b/2.3-1/src/c/linearAlgebra/svd/zsvda.c index 0d360222..12a07aaf 100644 --- a/2.3-1/src/c/linearAlgebra/svd/zsvda.c +++ b/2.3-1/src/c/linearAlgebra/svd/zsvda.c @@ -113,7 +113,16 @@ void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *o out3[i+j*N] = zconjs(VT[j+i*N]); out3[j+i*N] = zconjs(VT[i+j*N]); } - } + } + /* output from zgesvd is copied to out2 variables in required format*/ + for(j=0;j<M;j++){ + for(k=0;k<N;k++){ + if(j == k) + out2[j*(Min(M,N))+k] = DoubleComplex(S[j],0); + else + out2[j*(Min(M,N))+k] = DoubleComplex(0,0); + } + } //ztransposea(VT,LDVT,Min(M,N),out3); /*for(i=0;i<N;i++){ for(j=0;j<N;j++){ @@ -124,7 +133,7 @@ void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *o //free(U); //free(VT); } - else{ + else{ /*svd(x,'e')*/ LDA = M; LDU = M; if(M > N){ @@ -156,18 +165,18 @@ void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *o for(j=0;j<N;j++){ out3[j+i*N] = zconjs(VT[i+j*Min(M,N)]); } - } + } + /* output from zgesvd is copied to out2 variables in required format*/ + 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] = DoubleComplex(S[j],0); + else + out2[j*(Min(M,N))+k] = DoubleComplex(0,0); + } + } //free(U); //free(VT); } - /* output from zgesvd is copied to out2 variables in required format*/ - 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] = DoubleComplex(S[j],0); - else - out2[j*(Min(M,N))+k] = DoubleComplex(0,0); - } - } } } |