SUBROUTINE REGRE(X,N,M,A,S,Q,U,RR,F)
! CHUONG TRINH NAY TINH CAC HE SO HOI QUI CUA PHUONG TRINH
! HOI QUI GIUA X1 VA CAC X2,...,XM VA XAC DINH CAC DAI LUONG
! CHO PHAN TICH PHUONG SAI
! DANG PHUONG TRINH HQ: X1=A1+A2*X2+...+AM*XM
!
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M LUU TRU SO LIEU
!            QUAN TRAC CUA CAC X1..XM (LUU THEO COT: x1,x2,...,xn la cot 1; xn+1,..., x2n la cot 2,...)
!            Cot 1 la bien phu thuoc X1, cac cot 2..M la M-1 bien doc lap
!        + N DUNG LUONG MAU
!        + M SO BIEN (1 BIEN PHU THUOC (X1) VA M-1 BIEN DOC LAP X2..XM))
! Output:+ A MANG MOT CHIEU KICH THUOC M LUU CAC HE SO
!            A1, A2,..., AM tuong ung la He so tu do va cac he so HQ cua cac bien X2,..., XM
!        + S SAI SO CHUAN (CHUAN SAI THANG DU)
!        + Q TONG BINH PHUONG CAC BIEN SAI THANG DU
!        + U TONG BINH PHUONG CAC BIEN SAI HOI QUI
!        + RR HE SO TUONG QUAN BOI
!        + F BIEN CO PHAN BO FISHER
!
DIMENSION X(N*M),A(M)
DIMENSION R(M*M),RXX(M*M),RY(M),TB(M),RX(M*M)
Real :: Tmp,Yhq
CALL COVAR(X,N,M,R,TB)
L=0
J=0
DO I=M+1,M*M
   IF (MOD(I,M).NE.1) THEN
      L=L+1
      RXX(L)=R(I)
   ELSE
      J=J+1
      RY(J)=R(I)
   ENDIF
ENDDO
RX=RXX
CALL MINVERT(RX,M-1)
CALL MULTIP(RX,RY,A,M-1,M-1,1)
TMP=0.0
DO J=1,M-1
   TMP=TMP+A(J)*TB(J+1)
ENDDO
TMP=TB(1)-TMP
DO J=M,2,-1
   A(J)=A(J-1)
ENDDO
A(1)=TMP

RX=R
D=DET(RX,M)
RX=R
D11=PPDS(RX,M,1,1)
TMP=D/D11

RR=SQRT(1-TMP/R(1))

Q=0.0
DO I=1,N
   YHQ=0.0
   DO J=2,M
      IJ=(J-1)*N+I
      YHQ=YHQ+A(J)*X(IJ)
   ENDDO
   YHQ=YHQ+A(1)
   TMP=X(I)-YHQ
   Q=Q+TMP*TMP
ENDDO
U=REAL(N)*R(1)-Q
F=U*REAL(N-M)/(Q*REAL(M-1))
S=SQRT(Q/REAL(N-M))

END SUBROUTINE REGRE

SUBROUTINE MINVERT(AA,N)
! CHUONG TRINH NAY TINH MA TRAN NGHICH DAO CUA MA TRAN A
! INPUT: + AA MANG MOT CHIEU KICH THUOC N*N LUU THEO COT CUA A
!        + N SO HANG/COT CUA MA TRAN A
! OUTPUT: AA MA TRAN NGHICH DAO CUA A
!
DIMENSION AA(N*N), A(N,N)

K=0
DO J=1,N
   DO I=1,N
      K=K+1
      A(I,J)=AA(K)
   ENDDO
ENDDO
M=N-1
DO 10 I=1,N
   R=A(I,1)
   DO K=1,M
      A(I,K)=A(I,K+1)/R
   ENDDO
   A(I,N)=1.0/R
   DO 10 J=1,N
      IF (I.NE.J) THEN
         R=A(J,1)
         DO K=1,M
            A(J,K)=A(J,K+1)-R*A(I,K)
         ENDDO
         A(J,N)=-R*A(I,N)
      ENDIF
  10  CONTINUE
K=0
DO J=1,N
   DO I=1,N
      K=K+1
      AA(K)=A(I,J)
   ENDDO
ENDDO
END SUBROUTINE MINVERT

SUBROUTINE MULTIP(A,B,C,N,M,L)
! CHUONG TRINH NAY TINH TICH CUA HAI MA TRAN
! INPUT: + A MANG MOT CHIEU KICH THUOC N*M (N HANG, M COT)
!            LUU MA TRAN A THEO COT
!        + B MANG MOT CHIEU KICH THUOC M*L (M HANG, L COT)
!            LUU MA TRAN B THEO COT
!        + N SO HANG CUA MA TRAN A
!        + M SO COT CUA A VA LA SO HANG CUA B
!        + L SO COT CUA B
! OUTPUT: + C MANG MOT CHIEU KICH THUOC N*L (N HANG, L COT)
!             LUU MA TRAN C THEO COT

DIMENSION A(N*M),B(M*L),C(N*L)

DO I=1,N
DO K=1,L
   IK=(K-1)*N+I
   C(IK)=0.0
   DO J=1,M
      IJ=(J-1)*N+I   !(I-1)*M+J
      JK=(K-1)*M+J
      C(IK)=C(IK)+A(IJ)*B(JK)
   ENDDO
ENDDO
ENDDO

END SUBROUTINE MULTIP

SUBROUTINE COVAR(X,N,M,R,TB)
! CHUONG TRINH NAY TINH MA TRAN TUONG QUAN CUA M BIEN
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M CHUA SO LIEU BAN DAU
!          LUU THEO COT CUA MA TRAN XX(N,M)
!        + N DUNG LUONG MAU
!        + M SO BIEN
! OUTPUT: R MANG MOT CHIEU KICH THUOC M*M CHUA MA TRAN TUONG QUAN
!
DIMENSION X(N*M),XX(N,M),R(M*M),RR(M,M),TB(M)

K=0
DO J=1,M
   TB(J)=0.0
   DO I=1,N
      K=K+1
      TB(J)=TB(J)+X(K)
      XX(I,J)=X(K)
   ENDDO
   TB(J)=TB(J)/REAL(N)
ENDDO
DO J=1,M
   DO K=J,M
      RR(J,K)=0.0
      DO I=1,N
         RR(J,K)=RR(J,K)+(XX(I,J)-TB(J))*(XX(I,K)-TB(K))
      ENDDO
      RR(J,K)=RR(J,K)/REAL(N)
      RR(K,J)=RR(J,K)
   ENDDO
ENDDO
L=0
DO K=1,M
   DO J=1,M
      L=L+1
      R(L)=RR(J,K)
   ENDDO
ENDDO

END SUBROUTINE COVAR

FUNCTION PPDS(A,N,IR,JC)
!
! HAM NAY TINH PHAN PHU DAI SO CUA PHAN TU B(IR,JC)
! (HANG IR, COT JC) CUA MA TRAN B(NxN) MA CAC PHAN TU
! CUA NO DUOC LUU TRONG MANG A(N*N) THEO QUI CACH COT
! TRUOC DONG SAU
! INPUT: + MANG A(N*N) CHUA MA TRAN DAU VAO CUA B
!        + N KICH THUOC CUA B
!        + IR CHI SO HANG
!        + JC CHI SO COT
! OUTPUT: PHAN PHU DAI SO CUA PHAN TU HANG IR COT JC
!
DIMENSION A(N*N), B(N,N)
K=0
DO J=1,N
   DO I=1,N
      K=K+1
      B(I,J)=A(K)
   ENDDO
ENDDO

K=0
DO J=1,N
   IF (J.NE.JC) THEN
   DO I=1,N
      IF (I.NE.IR) THEN
         K=K+1
         A(K)=B(I,J)
      ENDIF
   ENDDO
   ENDIF
ENDDO

PPDS=DET(A,N-1)*(-1)**(IR+JC)

END FUNCTION PPDS


FUNCTION DET(X,N)
! HAM NAY TINH DINH THUC CUA MA TRAN A(N,N)
! INPUT: + MANG X DO DAI N*N CHUA CAC PHAN TU CUA MA TRAN A
!        + N KICH THUOC MA TRAN A
!   CHU Y: X(1)=A(1,1), X(2)=A(2,1),..., X(N)=A(N,1),X(N+1)=A(1,2),
!          X(N+2)=A(2,2),...
! OUTPUT: DINH THUC CUA A
!
PARAMETER (EP=1.0E-6)
DIMENSION X(N*N),A(N,N)
K=0
DO J=1,N
   DO I=1,N
      K=K+1
      A(I,J)=X(K)
   ENDDO
ENDDO
D=1.0
N1=N-1
DO 10 K=1,N1
   AM=0.0
   DO 11 I=K,N
      T=A(I,K)
      IF (ABS(T).GE.ABS(AM)) THEN
         AM=T
         J=I
      ENDIF
  11     CONTINUE
         IF (ABS(AM).LE.EP) THEN
      DT=0.0
      DET=DT
      RETURN
   ELSE
      IF (J.NE.K) THEN
         D=-D
         DO I=K,N
            T=A(J,I)
            A(J,I)=A(K,I)
            A(K,I)=T
         ENDDO
      ENDIF
   ENDIF
   M=K+1
   DO I=M,N
      T=A(I,K)/AM
      DO J=M,N
         A(I,J)=A(I,J)-T*A(K,J)
      ENDDO
   ENDDO
   D=D*A(K,K)
  10  CONTINUE
      DT=D*A(N,N)
DET=DT

END FUNCTION DET