    MODULE interface_definitions
    interface
    subroutine UV_f77(U,V,m,n,s,Us,lda,ldb,jpvt,rowindex,work,lwork,rwork)    

 implicit none
    
    
    integer(4) m,n,s,lda,ldb,lwork
    complex(4), dimension(:,:)::U
    complex(4) V(ldb,*)
    complex(4) US(lda,*)
    complex(4), dimension(:)::work
    real(4), dimension(:)::rwork
    integer(4), dimension(:)::rowindex,jpvt
    real(4) rcond
    integer(4) rank,info
    
    integer(4) i,j

    END
    
    
    subroutine calculate_interaction(Q,m,n)
        implicit none
        integer(4) m,n
        complex(4) Q(:,:)
        integer(4) i,j
        real(4) dx,dy,r,x,y
    end
    
    end interface
    end MODULE interface_definitions

   
    subroutine UV_f77(U,V,m,n,s,Us,lda,ldb,jpvt,rowindex,work,lwork,rwork) 
    
    implicit none
    include 'mkl.fi'
   
    
    integer(4) m,n,s,lda,ldb,lwork
    complex(4), dimension(:,:)::U
    complex(4) V(ldb,*)
    complex(4) US(lda,*)
    complex(4), dimension(:)::work
    real(4), dimension(:)::rwork
    integer(4), dimension(:)::rowindex,jpvt
    real(4) rcond
    integer(4) rank,info
    
    integer(4) i,j
    
    do i=1,s
        do j=1,s
            Us(j,i)=U(rowindex(j),i) 
        enddo
    enddo    
    
    
    jpvt(1:s)=0
    rcond=1e-5
    
    work=cmplx(0)
    rwork=0
    call cgelsy(s, s, n, Us, lda, V, s, jpvt, rcond, rank, work, lwork, rwork, info)
    end
    
    module test_types
    type QR_pair_test
        complex(4), pointer::QLE(:,:),RLE(:,:)
    end type 
    end
    
    subroutine calculate_interaction(Q,m,n)
        implicit none
        integer(4) m,n
        complex(4) Q(:,:)
        integer(4) i,j
        real(4) dx,dy,r,x,y
        dx=0.2
        dy=0.2
        
            do j=1,n
                do i=1,m
                x=i*dx+10
                y=j*dy
                r=sqrt(x*x+y*y)
                Q(i,j)=exp(cmplx(0,1)*2*3.14*r)/r
            enddo
        enddo
        
        
    end
    
    subroutine test_QR_fact()
    use interface_definitions
    use test_types
    implicit none
    integer(4) num_QRS
    complex(4),target, allocatable::QLE(:,:)
    complex(4),target, allocatable::RLE(:,:)
    integer(4) rank
    integer(4),  allocatable::localindex(:)
    complex(4), target,allocatable::Us_temp1(:,:)
    complex(4), pointer::Ustemp1(:,:)
    integer(4), target,allocatable::jpvt(:)
    integer(4), pointer::jpvt1(:)
    complex(4), allocatable,dimension(:)::work
    real(4), allocatable,dimension(:)::rwork
    integer(4) lwork
    type(QR_pair_test) QR(300)
    integer(4) maxnumr,maxnumc,maxrank,numr,numc,ii,iii
    logical(1) zeromatrix
    integer(4) number_threads,num_QRS_NUM_THREADS,lp1,lp2
    logical run_parallel
    
    
    
    
    
    number_threads=1
    run_parallel=.true.
    num_QRS=300 
    
    
    
    num_QRS_NUM_THREADS=num_QRS/number_threads
    maxnumr=200
    maxnumc=198
    maxrank=72
    
    
    
    
    !$OMP parallel IF(run_parallel .eq. .true.) &
    !$OMP private(numr,numc,rank,QLE,RLE, &
    !$OMP localindex,ii,lp1,&
    !$OMP Us_temp1,Ustemp1,jpvt,jpvt1,work,rwork,lwork) &
    !$OMP firstprivate(maxnumr,maxnumc,maxrank) &
    !$OMP NUM_THREADS(number_threads) 
    !$OMP do
    
    do lp2=1,number_threads
    
    if(run_parallel .eq. .true. .or. lp2 .eq. 1)then
        lwork=maxrank+2*maxrank+1+1000000
        allocate(Us_temp1(maxrank,maxrank),jpvt(maxnumc*3),work(lwork),rwork(2*maxrank))
        allocate(localindex(1:maxrank))
        do iii=1,maxrank
            localindex(iii)=iii
        enddo
    endif
    
    
    
    do lp1=1,num_QRS_NUM_THREADS+1 
    
    
    
    ii=(lp1-1)*number_threads+lp2
    
    
    if(ii .le. num_QRS)then
    
    
    numr=maxnumr/2
    numc=maxnumc/2
    rank=maxrank/2
    
    Ustemp1=>Us_temp1
    jpvt1=>jpvt
    
    allocate(QR(ii)%QLE(numr,rank),QR(ii)%RLE(rank,numc))
    call calculate_interaction(QR(ii)%QLE,numr,rank)
    QR(ii)%RLE=cmplx(100)
    
    
    call UV_f77(QR(ii)%QLE,QR(ii)%RLE,numr,numc,rank,Us_temp1,maxrank,rank,jpvt,localindex,work,lwork,rwork)
    
    
    
    
    end if
    enddo
    if(run_parallel .or. lp2 .eq. 1)then
        deallocate(Us_temp1,jpvt,work,rwork)           
    endif
    
    
    enddo 
    !$OMP end do
    !$OMP end parallel
    
    
    end
    
    program test_qr
        call test_QR_fact()
    end
    