summaryrefslogtreecommitdiff
path: root/2.3-1/src/c/linearAlgebra/qr
diff options
context:
space:
mode:
Diffstat (limited to '2.3-1/src/c/linearAlgebra/qr')
-rw-r--r--2.3-1/src/c/linearAlgebra/qr/dqra.c298
1 files changed, 298 insertions, 0 deletions
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;
+}
generated by cgit (git 2.34.1) at 2024-12-21 13:44:06 +0000