Least Square Method for UH from DRH and ERH

dimension Q(100,100),R(100,100),RTranspose(100,100),RTransposeR(100,100),RTransposeRInverse(100,100),RTransposeQ(100,100),U(100,100)
    integer qrow,qcol,rrow,rcol
    qrow=16 !16 non-zero DRH ordinates
    qcol=1
    rrow=16 !16 rows of Rainfall matrix
    rcol=14 !14 rows of Rainfall matrix (qrow-(no. of ERH ordinates-1)) = 16-(3-1) = 16-2 = 14
    open(5,file='Q.txt')
    do i=1,qrow
        read (5,*)(Q(i,j),j=1,qcol)
    end do
    open(5,file='R.txt')
    do i=1,rrow
        read (5,*)(R(i,j),j=1,rcol)
    end do

    call Transpose(R,RTranspose,rrow,rcol)
    call MatMul(RTranspose,R,RTransposeR,rcol,rrow,rrow,rcol)
    call Inverse(RTransposeR,RTransposeRInverse,rcol)
    call MatMul(RTranspose,Q,RTransposeQ,rcol,rrow,qrow,qcol)
    call MatMul(RTransposeRInverse,RTransposeQ,U,rcol,rcol,rcol,1)

    open(5,file='U.txt')
    do i=1,rcol
        write (5,*) U(i,1)
    end do

    end program


subroutine Transpose(A,R,rows,cols) !rows,cols = of A
integer rows,cols
dimension A(100,100),R(100,100) !CAUTION: This declaration is crucial
    do i=1,cols !i goes from 1 to cols
        do j=1,rows !and j from 1 to rows because transpose should have the dimension: cols x rows
            R(i,j)=A(j,i) !fills from column of A to row of R
        end do
    end do
    return
end subroutine Transpose


subroutine MatMul(A,B,C,arow,acol,brow,bcol)
integer arow,acol,brow,bcol
!To see how this works, try writing the elements of the product matrix in sigma notation
dimension A(100,100),B(100,100),C(100,100) !CAUTION: This declaration is crucial
do k=1,arow !k runs from 1 to arow
    do l=1,bcol !and l from 1 to bcol because that's the output matrix dimension and we _have_ to populate each element of it
        sum=0
        do i=1,acol !this loop is the essence of the aforementioned sigma notation for the (k,l) element of the product matrix
            sum=sum+A(k,i)*B(i,l)
        end do
        C(k,l)=sum
    end do
end do
return
end subroutine MatMul

subroutine Inverse(X,XI,size)
dimension X(100,100),XI(100,100)
integer rowindex,colindex,size,tmpi,tmpj,i,j
real pivotelement

do i=1,size
    do j=1,size
        if(i.EQ.j) then
            XI(i,j)=1
        endif
    end do
end do


do rowindex=1,size !loop over rows of X

    pivotelement=X(rowindex,rowindex) !this is crucial!

    !divide the row by pivot element
    do colindex=1,size
        X(rowindex,colindex)=X(rowindex,colindex)/pivotelement
    end do
    do colindex=1,size
        XI(rowindex,colindex)=XI(rowindex,colindex)/pivotelement
    end do


    !perform row operations to make non-pivot elements of this column zero
    do i=1,size !i gives the row number being manipulated
        pivotelement=X(i,rowindex) !pivot element for the currently manipulated row. this i crucial

        if(i.NE.rowindex) then !row operations are for non-current rows only

            do j=1,size !column index for ith row of X matrix
                X(i,j)=X(i,j)-pivotelement*X(rowindex,j)
            end do
            do j=1,size !column index for ith row of XI matrix
                XI(i,j)=XI(i,j)-pivotelement*XI(rowindex,j)
            end do


        endif
    end do



end do


end subroutine