summaryrefslogtreecommitdiff
path: root/2.3-1/src/c/linearAlgebra/svd
diff options
context:
space:
mode:
Diffstat (limited to '2.3-1/src/c/linearAlgebra/svd')
-rw-r--r--2.3-1/src/c/linearAlgebra/svd/dsvda.c126
-rw-r--r--2.3-1/src/c/linearAlgebra/svd/zsvda.c173
2 files changed, 299 insertions, 0 deletions
diff --git a/2.3-1/src/c/linearAlgebra/svd/dsvda.c b/2.3-1/src/c/linearAlgebra/svd/dsvda.c
new file mode 100644
index 00000000..e6af3008
--- /dev/null
+++ b/2.3-1/src/c/linearAlgebra/svd/dsvda.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: Sandeep Gupta
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+
+ */
+
+/*Funtion to find singular value decomposition of given matrix */
+
+#include "lapack.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include "string.h"
+#include <math.h>
+#include "svd.h"
+#include "matrixTranspose.h"
+
+int min(int a,int b);
+int max(int a,int b);
+
+extern double dgesvd_(char*,char*,int*,int*,double*,int*,double*,double*,int*,\
+ double*,int*,double *,int*,int*);
+
+/* 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<M;j++){
+ for(k=0;k<N;k++){
+ if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j);
+ else *((out2+j*(min(M,N)))+k) = 0;
+ }
+ }
+ dtransposea(VT,LDVT,N,out3);
+ }
+ else{
+ memcpy(out1,U,LDU*min(row,col)*sizeof(double));
+ for(j=0;j<min(M,N);j++){
+ for(k=0;k<min(M,N);k++){
+ if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j);
+ else *((out2+j*(min(M,N)))+k) = 0;
+ }
+ }
+ dtransposea(VT,LDVT,N,out3);
+ }
+}
+
+int min(int a,int b){
+ if(a > b) return b;
+ return a;
+}
+
+int max(int a,int b){
+ if(a > b) return a;
+ return b;
+}
diff --git a/2.3-1/src/c/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 <stdio.h>
+#include <stdlib.h>
+#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<Min(M,N);i++){
+ out2[i] = DoubleComplex(S[i],0);
+ //out2[i] = S[i];
+ //out2[i] = 0;
+ }
+ out1 = NULL;
+ out3 = NULL;
+ //for(i=0;i<Min(M,N);i++) printf("%lf ",S[i]);
+ //free(S);
+ }
+ else if(nout == 3){
+ if(in2 == 0 || M == N){
+ JOBU = 'A';
+ JOBVT = 'A';
+ LDA = M;
+ LDU = M;
+ LDVT = N;
+ doubleComplex *U = malloc(LDU*M*sizeof(doubleComplex));
+ doubleComplex *VT = malloc(LDVT*N*sizeof(doubleComplex));
+
+ /*doubleComplex wopt;
+ LWORK = -1;
+ zgesvd_(&JOBU,&JOBVT,&M,&N,buf,&LDA,S,U,&LDU,VT,&LDVT,&wopt,&LWORK,RWORK,&INFO);*/
+
+ //LWORK = (int)zreals(wopt);
+
+ WORK = (doubleComplex *)malloc(LWORK*sizeof(doubleComplex));
+ zgesvd_(&JOBU,&JOBVT,&M,&N,buf,&LDA,S,U,&LDU,VT,&LDVT,WORK,&LWORK,RWORK,&INFO);
+
+ memcpy(out1,U,LDU*Min(M,N)*sizeof(doubleComplex));
+ //memcpy(out3,VT,N*N*sizeof(doubleComplex));
+ for(i=0;i<N;i++){
+ for(j=i;j<N;j++){
+ out3[i+j*N] = zconjs(VT[j+i*N]);
+ out3[j+i*N] = zconjs(VT[i+j*N]);
+ }
+ }
+ //ztransposea(VT,LDVT,Min(M,N),out3);
+ /*for(i=0;i<N;i++){
+ for(j=0;j<N;j++){
+ printf("[ %lf %lf]",zreals(VT[i*N+j]),zimags(VT[i*N+j]));
+ }
+ printf("\n");
+ }*/
+ //free(U);
+ //free(VT);
+ }
+ else{
+ LDA = M;
+ LDU = M;
+ if(M > 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<Min(M,N);i++){
+ for(j=0;j<N;j++){
+ out3[j+i*N] = zconjs(VT[i+j*Min(M,N)]);
+ }
+ }
+ //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);
+ }
+ }
+ }
+}
generated by cgit (git 2.34.1) at 2024-12-21 04:47:37 +0000