program f90_PCA
    !********************************************
    ! Principal Component Analysis
    !********************************************
    implicit none
    integer::knor                  ! Index for normalization (0)
    integer::kvoc                  ! Index for Var(0) or Cor(1)
    character::strcom*100          ! Comment
    integer::mcol                  ! Number of variables
    integer::ndata                 ! Number of data sets
    real(8),allocatable::xd(:,:)   ! Array for input data
    real(8),allocatable::varm(:,:) ! Array for variance-covariance matrix
    real(8),allocatable::corm(:,:) ! Array for correlation matrix
    real(8),allocatable::ary(:,:)  ! Array for eigen value analysis
    real(8),allocatable::ev(:)     ! Eigenvalue
    real(8),allocatable::vec(:,:)  ! Eigenvector
    real(8),allocatable::zd(:,:)   ! Principal component score
    real(8),allocatable::tr        ! Sum total of eigenvalues

    integer::i,j,k,nt
    character(len=50)::fnameR,fnameW
    character(len=100)::dummy
    character(len=1000)::linebuf

    call getarg(1,dummy);read(dummy,*) knor
    call getarg(2,dummy);read(dummy,*) kvoc
    call getarg(3,fnameR)
    call getarg(4,fnameW)

    !Data input
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) mcol,ndata
    nt=mcol
    allocate(xd(1:ndata,1:nt),zd(1:ndata,1:nt))
    allocate(varm(1:nt,1:nt),corm(1:nt,1:nt),ary(1:nt,1:nt))
    allocate(ev(1:nt),vec(1:nt,1:nt))
    do k=1,ndata
        read(11,*) (xd(k,j),j=1,mcol)
    end do
    close(11)

    !Normalization
    if(knor==0)call STAND(ndata,mcol,xd)

    !Variance-covariance matrix, Correlation matrix
    call VARCOR(ndata,mcol,xd,varm,corm)
    if(kvoc==0)then
        do i=1,mcol
            do j=1,mcol
                ary(i,j)=varm(i,j)
            end do
        end do
    else
        do i=1,mcol
            do j=1,mcol
                ary(i,j)=corm(i,j)
            end do
        end do
    end if

    !Sum total of diagonal elements
    tr=0.0D0
    do i=1,mcol
        tr=tr+ary(i,i)
    end do

    !Eigenvalue, eigenvector
    call JACOBI(mcol,ary,ev,vec,nt)

    !Principal component score
    do k=1,ndata
        do i=1,mcol
            zd(k,i)=0.0D0
            do j=1,mcol
                zd(k,i)=zd(k,i)+xd(k,j)*vec(j,i)
            end do
        end do
    end do

    do i=1,mcol
        write(6,'("ev(",i0,")=",e15.7)') i,ev(i)
    end do

    !Print out
    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)

    write(dummy,*) 'No.'
    linebuf=trim(dummy)
    do j=1,mcol
        write(dummy,*) j
        linebuf=trim(linebuf)//','//trim(dummy)
    end do
    call del_spaces(linebuf)
    write(12,*) trim(linebuf)
    write(dummy,*) 'Eigen_value'
    linebuf=trim(dummy)
    do j=1,mcol
        write(dummy,*) ev(j)
        linebuf=trim(linebuf)//','//trim(dummy)
    end do
    call del_spaces(linebuf)
    write(12,*) trim(linebuf)
    write(dummy,*) 'Proportion'
    linebuf=trim(dummy)
    do j=1,mcol
        write(dummy,*) ev(j)/tr
        linebuf=trim(linebuf)//','//trim(dummy)
    end do
    call del_spaces(linebuf)
    write(12,*) trim(linebuf)
    do i=1,mcol
        write(dummy,*) 'Eigen_vec(',i,')'
        linebuf=trim(dummy)
        do j=1,mcol
            write(dummy,*) vec(i,j)
            linebuf=trim(linebuf)//','//trim(dummy)
        end do
        call del_spaces(linebuf)
        write(12,*) trim(linebuf)
    end do
    do i=1,ndata
        write(dummy,*) 'Score(',i,')'
        linebuf=trim(dummy)
        do j=1,mcol
            write(dummy,*) zd(i,j)
            linebuf=trim(linebuf)//','//trim(dummy)
        end do
        call del_spaces(linebuf)
        write(12,*) trim(linebuf)
    end do
    close(12)

stop

contains

    !**************************************************
    ! Normalization
    !**************************************************
    subroutine STAND(ndata,mcol,xd)
        integer,intent(in)::ndata,mcol
        real(8),intent(inout)::xd(1:ndata,1:mcol)

        integer::i,k
        real(8)::xm(1:mcol),sd(1:mcol)

        !Average
        do i=1,mcol
            xm(i)=0.0D0
            do k=1,ndata
                xm(i)=xm(i)+xd(k,i)
            end do
            xm(i)=xm(i)/dble(ndata)
        end do
        !Standard deviation
        do i=1,mcol
            sd(i)=0.0D0
            do k=1,ndata
                sd(i)=sd(i)+(xd(k,i)-xm(i))*(xd(k,i)-xm(i))
            end do
            sd(i)=sqrt(sd(i)/dble(ndata-1))
        end do
        !Normalization
        do i=1,mcol
            do k=1,ndata
                xd(k,i)=(xd(k,i)-xm(i))/sd(i)
            end do
        end do
    end subroutine STAND
    !**************************************************
    ! Variance-covariance matrix, Correlation matrix
    !**************************************************
    subroutine VARCOR(nrow,ncol,xd,varm,corm)
        integer,intent(in)::nrow,ncol
        real(8),intent(in)::xd(1:nrow,1:ncol)
        real(8),intent(out)::varm(1:ncol,1:ncol),corm(1:ncol,1:ncol)
        integer::i,j,k,nn,mm
        real(8)::xm(1:ncol),wmat(1:nrow,1:ncol)

        nn=nrow
        mm=ncol
        do i=1,nn
            do j=1,mm
                wmat(i,j)=xd(i,j)
            end do
        end do
        do j=1,mm
            xm(j)=0.0D0
            do i=1,nn
                xm(j)=xm(j)+wmat(i,j)
            end do
            xm(j)=xm(j)/dble(nn)
        end do
        do i=1,mm
            do j=1,mm
                varm(i,j)=0.0D0
                do k=1,nn
                    varm(i,j)=varm(i,j)+(wmat(k,i)-xm(i))*(wmat(k,j)-xm(j))
                end do
                varm(i,j)=varm(i,j)/dble(nn-1)
            end do
        end do
        do i=1,mm
            do j=1,mm
                corm(i,j)=varm(i,j)/sqrt(varm(i,i))/sqrt(varm(j,j))
            end do
        end do
    end subroutine VARCOR
    !**************************************************
    ! Jacobi method for eigen value analysis
    !**************************************************
    subroutine JACOBI(n,array,ev,vec,nt)
        !Jacobi method by Dr.Okumura
        !Eigenvalue: ev(k), eigenvector: vec(k,i)

        integer,intent(in)::n,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        real(8),intent(out)::ev(1:nt),vec(1:nt,1:nt)

        integer::i,j,k,im
        real(8)::s,c,aoff,eps,t
        real(8)::aa,aj,ak,vj,vk,vmax,bb

        do i=1,n
            do j=1,n
                vec(i,j)=0.0
            end do
            vec(i,i)=1.0
        end do

        aoff=0.0
        do j=1,n-1
            do k=j+1,n
                aoff=aoff+array(j,k)*array(j,k)
            end do
        end do
        eps=aoff*2.0
        do j=1,n
            eps=eps+array(j,j)*array(j,j)
        end do
        eps=eps*0.00000000005

        do while(aoff>eps)
            do j=1,n-1
                do k=j+1,n
                    if(array(j,k)==0.0)then
                        t=0.0
                    else
                        aa=(array(k,k)-array(j,j))/(2.0*array(j,k))
                    end if
                    if(aa>=0)then
                        t=1.0/(aa+sqrt(1.0+aa*aa))
                    else
                        t=1.0/(aa-sqrt(1.0+aa*aa))
                    end if
                    c=1.0/sqrt(1.0+t*t)
                    s=t*c
                    array(j,j)=array(j,j)-t*array(j,k)
                    array(k,k)=array(k,k)+t*array(j,k)
                    array(j,k)=0.0
                    do i=1,j-1
                        aj=array(i,j)
                        ak=array(i,k)
                        array(i,j)=aj*c-ak*s
                        array(i,k)=aj*s+ak*c
                    end do
                    do i=j+1,k-1
                        aj=array(j,i)
                        ak=array(i,k)
                        array(j,i)=aj*c-ak*s
                        array(i,k)=aj*s+ak*c
                    end do
                    do i=k+1,n
                        aj=array(j,i)
                        ak=array(k,i)
                        array(j,i)=aj*c-ak*s
                        array(k,i)=aj*s+ak*c
                    end do
                    do i=1,n
                        vj=vec(i,j)
                        vk=vec(i,k)
                        vec(i,j)=vj*c-vk*s
                        vec(i,k)=vj*s+vk*c
                    end do
                end do
            end do
            aoff=0.0
            do j=1,n-1
                do k=j+1,n
                    aoff=aoff+array(j,k)*array(j,k)
                end do
            end do
        end do

        do k=1,n
            ev(k)=array(k,k)
        end do

        do k=1,n-1
            vmax=-1.0E+30
            do i=k,n
                if(ev(i)>=vmax)then
                    im=i
                    vmax=ev(i)
                end if
            end do
            aa=ev(k)
            ev(k)=ev(im)
            ev(im)=aa
            do i=1,n
                bb=vec(i,k)
                vec(i,k)=vec(i,im)
                vec(i,im)=bb
            end do
        end do
    endsubroutine JACOBI


    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_PCA