Least Square Regression for IDF curve equation

!Least Square Regression for IDF curve function I=a/(t+b)^c where a,b and c are constants
    dimension Q(100,100),R(100,100),RTranspose(100,100),RTransposeR(100,100),RTransposeRInverse(100,100),RTransposeQ(100,100),U(100,100)
    dimension iIntensity(100),Times(100)
    integer qrow,qcol,rrow,rcol
    real intensity,t,b,tmp,a,c,sum

    !here, Q vector corresponds to rainfall intensities
    ![Q]_nx1 = [R]_nx2 [X]_nx1
    ![X]=[lna -c]'_2x1

    qrow=7 !7 rows corresponding to 7 I's
    qcol=1 !1 column of I matrix coz it's a vector
    rrow=7 !7 rows of R matrix
    rcol=2 !2 cols of R matrix

    open(6,file="Result.txt")

    do b=8.0,14.0,0.5
        open(5,file='I.txt')
        do i=1,qrow
            read (5,*) intensity
            Q(i,1)=log(intensity)
        end do
        open(5,file='t.txt')
        do i=1,rrow
            read (5,*) t
            R(i,1)=1
            R(i,2) =log(t+b)
        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)

        a=exp(U(1,1))
        c=-U(2,1)


        !read the given observed intensities into the iIntensity array
        open(5,file='I.txt')
        do i=1,qrow
            read (5,*) iIntensity(i)
        end do
        !read the given observed times into the Times array
        open(5,file='t.txt')
        do i=1,qrow
            read (5,*) Times(i)
        end do

        !find squared error
        sum=0
        do i=1,qrow
            tmp=iIntensity(i)-a/(Times(i)+b)**c
            tmp=tmp**2
            sum=sum+tmp
        end do

        write (6,*) "b = ",b
        write (6,*) "a = ",a
        write (6,*) "c = ",c
        write (6,*) "Del I squared = ",sum
        write (6,*)
    end do

    close(6)

    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,rrow,rcol,brow,bcol)
integer rrow,rcol,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,rrow !k runs from 1 to rrow
    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,rcol !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
        else
            XI(i,j)=0 !this is for good measure coz by default, XI may have been changed by the caller itself though it is 0 by default
        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