Recursive Least Square Technique for IDF curve equation

!recursive Least Square technique for IDF curve equation I=a/(t+b)^c

integer k,j
    real t,intensity
    real X(100,100) !the parameter vector
    integer n !number of data
    real A(100,100) !design/coefficient matrix
    real I(100,100) !output vector
    real b
    real ATranspose(100,100),ATransposeA(100,100),ATransposeAInverse(100,100),ATransposeI(100,100)
    real ARow(100,100),ARowTranspose(100,100),ARowTransposeARow(100,100),P_n_plus1(100,100),P_n_plus1_ARowTranspose(100,100),totalCorrection(100,100)
    real ARowX(100,100)
    real tmp(100,100)
    b=12
    n=0

    do while(.TRUE.)
        write(*,*) "t"
        read(*,*) t
        write(*,*) "I"
        read(*,*) intensity
        n=n+1

        I(n,1)=log(intensity)
        A(n,1)=1
        A(n,2)=log(t+b)



        if(n.EQ.2) then
            !solve the system only once, when the number of data provided=number of unknowns(here, 2)
            call matrixTranspose(A,ATranspose,n,2)
            call matrixMul(ATranspose,A,ATransposeA,2,n,n,2)
            call Inverse(ATransposeA,ATransposeAInverse,2)
            call matrixMul(ATranspose,I,ATransposeI,2,n,n,1)
            call matrixMul(ATransposeAInverse,ATransposeI,X,2,2,2,1)
        elseif(n.GT.2) then
            !use recursive Least Square formula
            !add elements of previous ATransposeA with (ARow'ARow) where '=transpose operation, ARow=row matrix
            ![1 ln(t_n+b)], t_n being the current(new) time data value
            ARow(1,1)=1
            ARow(1,2)=log(t+b)
            call matrixTranspose(ARow,ARowTranspose,1,2)
            call matrixMul(ARowTranspose,ARow,ARowTransposeARow,2,1,1,2)
            do k=1,2
                do j=1,2
                    tmp(k,j)=ATransposeA(k,j) + ARowTransposeARow(k,j)
                end do
            end do
            call Inverse(tmp,P_n_plus1,2)
            !update ATransposeA
            do k=1,2
                do j=1,2
                    ATransposeA(k,j)=tmp(k,j)
                end do
            end do
            call matrixMul(P_n_plus1,ARowTranspose,P_n_plus1_ARowTranspose,2,2,2,1) !P_n_plus1_ARowTranspose is the same dimension as X
            call matrixMul(ARow,X,ARowX,1,2,2,1) !ARowX(1,1) is the output, which is a scalar
            do k=1,2
                totalCorrection(k,1)=P_n_plus1_ARowTranspose(k,1)*(log(intensity)-ARowX(1,1))
            end do
            !update X vector
            do k=1,2
                X(k,1)=X(k,1)+totalCorrection(k,1)
            end do
        endif

        write(*,*) "a = ",exp(X(1,1))
        write(*,*) "c = ",-X(2,1)
        write(*,*)


    end do

    end program


subroutine matrixTranspose(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 matrixTranspose


subroutine matrixMul(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 matrixMul

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

originalX=X

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

!restore X to original X since upon inversion completion, X has been reduced to an Identity matrix. good practice
X=originalX

end subroutine