! 8.3.f90
! 2024-10-24
! $Id: 8.3.f90 1.1 2024/11/02 12:56:41 s Exp $

implicit none
integer i,j,m,n
real s1(10),s2(10),s3(10),s4(10),s(4,10)
real average(4),variance(4)
real covariance(4,4), a(4,4)
real eigenvalue(4),eigenvector(4,4)
real ave1,avei,avej

! input
m=4
n=10
s1 = [11.91, 18.37,  3.64, 24.37, 30.42, &
      -1.45, 20.11,  9.28, 17.63, 15.71]
s2 = [29.59, 15.25,  3.53, 17.67, 12.74,&
      -2.56, 25.46,  6.92,  9.73, 25.09]
s3 = [23.27, 19.47, -6.58, 15.08, 16.24,&
     -15.05, 17.80, 18.82,  3.05, 16.94]
s4 = [27.24, 17.05, 10.20, 20.26, 19.84,&
       1.51, 12.24, 16.12, 22.93,  3.49]
   
s(1,1:n) = s1(1:n)
s(2,1:n) = s2(1:n)
s(3,1:n) = s3(1:n)
s(4,1:n) = s4(1:n)

! calculate
! average
do i=1,m
  average(i)=sum(s(i,1:n))/n
end do

! variance
do i=1,m
  avei=average(i)
  variance(i)=dot_product(s(i,1:n)-avei,s(i,1:n)-avei)/(n-1)
end do

! covariance
do i=1,m
  do j=1,m
    avei=average(i)
    avej=average(j)
    covariance(i,j)=dot_product(s(i,1:n)-avei,s(j,1:n)-avej)/(n-1)
  end do
end do

! eigenvector
a=covariance
call JAC(a,eigenvalue,eigenvector,m,m)

! output
print *,'8.3.f90'
print '(a6)','s'
do i = 1,m
  print '(10f6.2)', s(i,1:n)
end do
print *,''

print '(a10,4f8.2)','mean',     average(1:m)
print '(a10,4f8.2)','variance', variance(1:m)
print *,''

print *,'covariance matrix'
do i=1,m
    print '(4f7.2)', covariance(i,1:m)
end do
print *,''

print '(a10,a15)','eigenvalue','eigenvector'
DO J=1,m
    PRINT '(f10.2,(4F10.2))',eigenvalue(j),(eigenvector(i,j), i=1,m)
END DO 
end

      SUBROUTINE JAC(A,E,V,N,NN) ! 出典 戸川隼人『数値計算』岩波書店, 1991, pp.122-123.
      INTEGER I,J,K,N,NN
      INTEGER P,Q,KAISUU,KMAX
      PARAMETER (KMAX=300)
      REAL A(NN,NN),E(NN),V(NN,NN)
      REAL EPS,BUNBO,R,T,S,C
      REAL APQ,APQMAX
      REAL AIP,AIQ,APJ,AQJ,VIP,VIQ
      EPS=0.0001
! 配列Vに単位行列を入れる
      DO 10 I=1,N
        DO 20 J=1,N
          V(I,J)=0.0
   20   CONTINUE
        V(I,I)=1.0
   10 CONTINUE
! ヤコビ法の反復
      DO 30 KAISUU=1,KMAX
        DO 40 P=1,N-1
          DO 50 Q=P+1,N
! 回転角の算出
            BUNBO=A(P,P)-A(Q,Q)
            IF(BUNBO.NE.0.0) THEN
              R=2.0*A(P,Q)/BUNBO
              T=0.5*ATAN(R)
            ELSE
              T=0.78539818
            END IF
            S=SIN(T)
            C=COS(T)
! 左からの乗算
            DO 60 J=1,N
              APJ=A(P,J)
              AQJ=A(Q,J)
              A(P,J)= APJ*C+AQJ*S
              A(Q,J)=-APJ*S+AQJ*C
   60       CONTINUE
! 右からの乗算
            DO 70 I=1,N
              AIP=A(I,P)
              AIQ=A(I,Q)
              A(I,P)= AIP*C+AIQ*S
              A(I,Q)=-AIP*S+AIQ*C
              VIP=V(I,P)
              VIQ=V(I,Q)
              V(I,P)= VIP*C+VIQ*S
              V(I,Q)=-VIP*S+VIQ*C
   70       CONTINUE
   50     CONTINUE
   40   CONTINUE
! 収束判定
        APQMAX=0.0
        DO 80 P=1,N-1
          DO 90 Q=P+1,N
            APQ=ABS(A(P,Q))
            IF(APQ.GT.APQMAX) APQMAX=APQ
   90     CONTINUE
   80   CONTINUE
        IF (APQMAX.LE.EPS) GOTO 100
   30 CONTINUE
  100 CONTINUE
! 固有値を移す
      DO 110 I=1,N
        E(I)=A(I,I)
  110 CONTINUE
      RETURN
      END
! eof