diff options
Diffstat (limited to '2.3-1/src/c/linearAlgebra/qr/dqra.c')
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/qr/dqra.c | 298 |
1 files changed, 0 insertions, 298 deletions
diff --git a/2.3-1/src/c/linearAlgebra/qr/dqra.c b/2.3-1/src/c/linearAlgebra/qr/dqra.c deleted file mode 100644 index bae4bc27..00000000 --- a/2.3-1/src/c/linearAlgebra/qr/dqra.c +++ /dev/null @@ -1,298 +0,0 @@ -/* 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; -} |