program f90_MRA
    !*********************************************
    ! Multiple Regression Analysis
    !*********************************************
    implicit none
    character::strcom*100        !Comment
    integer::mcol                !Number of explanatory variables
    integer::ndata               !Number of data sets
    real(8),allocatable::xd(:,:) !Array for explanatory variables
    real(8),allocatable::yd(:)   !Array for response variables
    real(8),allocatable::ary(:,:)!Array for normal equation
    real(8),allocatable::bb(:)   !Array for partial regression coefficients
    real(8),allocatable::ye(:)   !Array for regression estimates
    real(8)::rr                  !Maltiple correlation coefficient
    real(8)::ydmean,yemean,y1,y2,y3
    integer::i,j,k,n,nt
    character::fnameR*50,fnameW*50,linebuf*1000

    call getarg(1,fnameR)
    call getarg(2,fnameW)

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) mcol,ndata
    nt=mcol+1
    n=nt
    allocate(xd(1:ndata,1:nt),yd(1:ndata))
    allocate(ary(1:nt,1:nt+1))
    allocate(bb(1:nt),ye(1:ndata))
    do k=1,ndata
        do j=1,nt
            xd(k,j)=1.0D0
        end do
    end do
    do k=1,ndata
        read(11,*) yd(k),(xd(k,j),j=2,nt)
    end do
    close(11)

    !Normal equation
    do i=1,nt
        do j=1,nt
            ary(i,j)=0.0D0
            do k=1,ndata
                ary(i,j)=ary(i,j)+xd(k,i)*xd(k,j)
            end do
        end do
    end do
    do i=1,nt
        ary(i,nt+1)=0.0D0
        do k=1,ndata
            ary(i,nt+1)=ary(i,nt+1)+xd(k,i)*yd(k)
        end do
    end do
    !To solve simultaneous linear equations
    call STEGJ(n,ary,nt)
    do i=1,nt
        bb(i)=ary(i,nt+1)
    end do
    !Multiple correlation coefficient
    ydmean=0.0D0
    yemean=0.0D0
    do k=1,ndata
        ye(k)=0.0D0
        do j=1,nt
            ye(k)=ye(k)+xd(k,j)*bb(j)
        end do
        ydmean=ydmean+yd(k)
        yemean=yemean+ye(k)
    end do
    ydmean=ydmean/dble(ndata)
    yemean=yemean/dble(ndata)
    y1=0.0D0
    y2=0.0D0
    y3=0.0D0
    do k=1,ndata
        y1=y1+(yd(k)-ydmean)*(ye(k)-yemean)
        y2=y2+(yd(k)-ydmean)*(yd(k)-ydmean)
        y3=y3+(ye(k)-yemean)*(ye(k)-yemean)
    end do
    rr=y1/sqrt(y2*y3)

    do i=1,nt
        write(6,'("B(",i0,")=",e15.7)') i-1,bb(i)
    end do
    write(6,'("rr=",f10.3)') rr

    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(linebuf,*) 'mcol=,',mcol
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(linebuf,*) 'ndata=,',ndata
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    do i=1,nt
        write(linebuf,*) 'B(',i-1,')=,',bb(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(linebuf,*) 'rr=,',rr
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write (12,'(a)') 'org.data,estimated'
    do k=1,ndata
        write(linebuf,*) yd(k),',',ye(k)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do

stop

contains

    subroutine STEGJ(n,ary,nt)
        !Gauss-Jordan
        !by Mr. Kasai
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt+1)

        integer::i,j,k,s
        real(8)::p,d,max,dumy

        do k=1,n
            max=0.0D0
            s=k
            do j=k,n
                if(abs(ary(j,k))>max)then
                    max=abs(ary(j,k))
                    s=j
                end if
            end do
            do j=1,n+1
                dumy=ary(k,j)
                ary(k,j)=ary(s,j)
                ary(s,j)=dumy
            end do
            p=ary(k,k)
            do j=k,n+1
                ary(k,j)=ary(k,j)/p
            end do
            do i=1,n
                if(i/=k)then
                    d=ary(i,k)
                    do j=k,n+1
                        ary(i,j)=ary(i,j)-d*ary(k,j)
                    end do
                end if
            end do
        end do
    end subroutine STEGJ


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


end program f90_MRA