***********************************************************************
********************                        ***************************
**************             MLS Ver.1.10          **********************
********************                        ***************************
***********************************************************************
C
      PROGRAM MLS
      PARAMETER (MAX=100)
      DIMENSION A(MAX,MAX),B(MAX,MAX),C(MAX,MAX),R(MAX,MAX)
      DIMENSION D(MAX,MAX),EN(MAX),Z(MAX),XX(MAX),XS(MAX),V(MAX)
      INTEGER X(MAX),JCH(MAX)
      CHARACTER QT*1,FILE*25
 1    ID=0
      ICOUNT=0
      CALL START(MENU)
**************************  FOR MENU 5 ( START UP )
      IF(MENU.EQ.5) THEN
         ICH=2
 2       WRITE(6,*) ' DATA 数を入力してください'
         READ(5,*,ERR=2) N
         IF(N.LE.0) GOTO 2
 3       WRITE(6,*) ' パラメータ数を入力してください'
         READ(5,*,ERR=3) KP
         WRITE(6,*) ' DATA FILE 名を入力してください'
         READ(5,301) FILE
         OPEN(1,FILE=FILE)
         DO 14 I=1,N
            READ(1,*) (D(I,J),J=1,KP),XX(I)
 14      CONTINUE
         CLOSE(1)
         CALL COPYMTX(XX,XS,N,1,MAX)
         GOTO 15
 111  END IF
**************************  END
      IF(MENU.NE.2) THEN
         CALL MAKE(A,XS,M,N,MENU,MAX,ICH)
      ELSE
**************************  FOR MENU 2 ( BIWEIGHT METHOD )
 18      WRITE(6,104)'DATA 数を入力して下さい'
         READ(5,*,ERR=18) M
         IF(M.LE.0) THEN
            WRITE(6,*) ' 入力しなおしてください！'
            GOTO 18
 24      END IF
         WRITE(6,*)'DATA の入力方法を指定してください'
 73      WRITE(6,*)'1:直接入力  2:DATA FILEによる入力'
         READ(5,*,ERR=73) IDATA
         IF(IDATA.LE.0.OR.IDATA.GE.3) GOTO 73
 74      IF(IDATA.EQ.2) THEN
            WRITE(6,*)'FILE の名前を入力してください'
            READ(5,301) FILE
            OPEN(1,FILE=FILE)
         ELSE
            WRITE(6,*)' DATA を X,Y の順に入力して下さい'
         END IF
         DO 4 I=1,M
         IF(IDATA.EQ.1) THEN
               READ(5,*) XX(I),EN(I)
         ELSE
               READ(1,*) XX(I),EN(I)
               IF(I.EQ.M) CLOSE(1)
         END IF
 4       CONTINUE
         WRITE(6,105)'反復数、傾きの順に出力します'
         CALL BEW(XX,EN,XS,V,Z,M)
         GOTO 68
**************************  END OF BIWEIGHT METHOD
 5       CONTINUE
      END IF
 15   IDA=0
      K=0
      IF(ICH.EQ.1) THEN
         K=M-1
      END IF
 11   K=K+1
      ICOUNT=ICOUNT+1
 12   CONTINUE
**************************  FOR MENU 5
      IF(MENU.EQ.5) THEN
         K=KP
         CALL MAKE2(A,K,N,D,JCH,MAX,ICOUNT)
      END IF
**************************  END
      CALL COPYMTX(A,B,N,K,MAX)
      CALL COPYMTX(XS,XX,N,1,MAX)
 46   CALL QR(B,R,X,EN,N,K,MAX,ID)
      IF(ID.EQ.1) GOTO 7
 25   CALL TRANSMTX(B,C,X,N,K,MAX)
      CALL TENMTX(C,B,N,K,MAX)
      CALL PRDMTX(B,XX,Z,K,N,1,MAX)
      CALL BACK(R,XX,Z,K,MAX)
      CALL TRANS(XX,V,X,K)
      CALL COPYMTX(V,XX,K,1,MAX)
      IF(MENU.EQ.1.AND.K.EQ.M) GOTO 6
      IF(IDA.EQ.1.OR.ICH.EQ.1) GOTO 6
      CALL PRDMTX(A,XX,V,N,K,1,MAX)
      CALL AIC(V,XS,K,N,AI)
      IF(MENU.NE.1) THEN
         IF(ICOUNT.EQ.1) THEN
            IAI=ICOUNT
            SAI=AI
         ELSE IF(AI+1.0.LT.SAI) THEN
            IAI=ICOUNT
            SAI=AI
         END IF
      END IF
      IF(MENU.EQ.5) THEN
         IF(ICOUNT.EQ.2**KP-1) THEN
            IDA=1
            ICOUNT=IAI
            GOTO 12
         ELSE
            GOTO 11
         END IF
      END IF
      IF(K.NE.M) GOTO 11
 17   K=IAI
      IDA=1
      IF(MENU.EQ.3) THEN
         L=K
      ELSE IF(MENU.EQ.4) THEN
         L=K-1
      END IF
      IF(ICH.EQ.2) THEN
         WRITE(6,201) L
      END IF
      GOTO 12
 6    CONTINUE
      IF(MENU.EQ.5) THEN
         CALL SOLVE5(XX,K,JCH)
      ELSE
         CALL SOLVE(XX,X,K,MENU)
      END IF
 68   WRITE(6,106)'以上の結果になりました'
 7    WRITE(6,103)'再び実行しますか？ ( Y or N )'
      READ(5,108) QT
      IF(QT.EQ.'N'.OR.QT.EQ.'n') THEN
         CALL MES
         STOP
      ELSE IF(QT.EQ.'Y'.OR.QT.EQ.'y') THEN
         GOTO 1
 8    ELSE
         GOTO 7
 9    END IF
      STOP
 103  FORMAT(2X,A34)
 104  FORMAT(/,2X,A27)
 105  FORMAT(/,2X,A32)
 106  FORMAT(/,5X,A35)
 108  FORMAT(A1)
 201  FORMAT( /,10X,I3,' 次式が最も良い近似です')
 301  FORMAT(A25)
      END
************************************************************************
*                  SUBROUTINE TO OUTPUT THE MESSAGES                   *
************************************************************************
      SUBROUTINE MES
         WRITE(6,*)'THANK YOU FOR USING MLS!  SEE YOU AGAIN!!'
         WRITE(6,*)'MERCI BEAUCOUP!!          AU REVOIR!!'
         WRITE(6,*)'DANKE!!                   AUF WIEDERSEHEN!!'
         WRITE(6,*)'謝謝!!                    再見!!'
         WRITE(6,*)'GRAZIE!!                  ARRIVEDERCI!!'
         WRITE(6,*)'GRACIAS!!                 HASTA LUEGO!!'
         WRITE(6,*)'OBRIGADO!!                ATE LOGO!!'
         WRITE(6,*)'СПАСИБО!!'
         WRITE(6,*)'БАЯРЛАЛАА!!'
      RETURN
      END
************************************************************************
*                 SUBROUTINE TO CHOOSE THE PROGRAM                     *
************************************************************************
      SUBROUTINE START(MENU)
      CHARACTER QT*1
      WRITE(6,101) '最小二乗法プログラム'
 19   WRITE(6,102)'近似する方法を選んで下さい'
      WRITE(6,*)
     *'1:多項式近似  2:比例式近似  3:連立方程式  4:重回帰モデル'
      READ(5,*,ERR=19) MENU
      IF(MENU.LE.0.OR.MENU.GE.5) GOTO 19
      IF(MENU.EQ.4) THEN
          MENU=5
          RETURN
      END IF
      MENU=4-MENU
      IF(MENU.EQ.3) THEN
 2       WRITE(6,103)'定数項は入れますか？( Y or N )'
         READ(5,108) QT
         IF(QT.EQ.'Y'.OR.QT.EQ.'y') THEN
            MENU=MENU+1
         ELSE IF(QT.NE.'N'.AND.QT.NE.'n') THEN
            GOTO 2
 3       END IF
      END IF
      RETURN
 101  FORMAT(/,5X,A46)
 102  FORMAT(/,2X,A30)
 103  FORMAT(2X,A34)
 108  FORMAT(A1)
      END
***********************************************************************
*      SUBROUTINE TO IMPUT THE DATA AND TO MAKE THE MATRIX            *
***********************************************************************
      SUBROUTINE MAKE(A,X,M,N,MENU,MAX,ICH)
      DIMENSION A(MAX,MAX),X(1)
      CHARACTER FILE*25
      IF(MENU.EQ.1) THEN
 11      WRITE(6,*)'パラメーターの数を入力して下さい'
         READ(5,*,ERR=11) M
         IF(M.LE.0) THEN
            WRITE(6,*)' 入力しなおしてください！'
            GOTO 11
 34      END IF
         N=M
         WRITE(6,*)'各列毎に、左辺の係数、右辺の数の順に入力して下さい'
         DO 2 I=1,M
            DO 1 J=1,M
               READ(5,*) A(I,J)
 1          CONTINUE
            READ(5,*) X(I)
 2       CONTINUE
      ELSE
 12      CONTINUE
         WRITE(6,*) ' DATA 数を入力して下さい'
         READ(5,*,ERR=12) N
         IF(N.LE.0) THEN
            WRITE(6,*) ' 入力しなおしてください！'
            GOTO 12
 29      END IF
 30      WRITE(6,*)' 最高次数決定の方法を選択して下さい'
         WRITE(6,*)' 1:直接入力  2:自動決定'
         READ(5,*,ERR=30) ICH
         IF(ICH.LE.0.OR.ICH.GE.3) THEN
            WRITE(6,*) ' 選択し直して下さい'
            GOTO 30
 31      END IF
         IF(ICH.EQ.2) THEN
            M=INT(2*SQRT(REAL(N)))
         ELSE
 9          WRITE(6,*) ' 求める最高次数を入力して下さい'
            READ(5,*,ERR=9) M
            IF(M.GE.N) THEN
               WRITE(6,*)'次数を入力しなおして下さい'
               GOTO 9
 7          END IF
            IF(MENU.EQ.4) THEN
               M=M+1
            END IF
         END IF
         WRITE(6,*)'DATA の入力方法を指定してください'
 23      WRITE(6,*)'1:直接入力  2:DATA FILEによる入力'
         READ(5,*,ERR=23) IDATA
         IF(IDATA.LE.0.OR.IDATA.GE.3) GOTO 23
 24      IF(IDATA.EQ.2) THEN
            WRITE(6,*)'FILE の名前を入力してください'
            READ(5,101) FILE
            OPEN(1,FILE=FILE)
         ELSE
            WRITE(6,*)' X,Y の順に入力して下さい'
         END IF
         DO 5 J=1,N
            IF(IDATA.EQ.1) THEN
               READ(5,*) XDATA,YDATA
            ELSE
               READ(1,*) XDATA,YDATA
            END IF
            IF(MENU.EQ.4) THEN
               DO 3 I=1,M
                  A(J,I)=XDATA**(I-1)
 3             CONTINUE
            ELSE
               DO 4 I=1,M
                  A(J,I)=XDATA**I
 4             CONTINUE
            END IF
            X(J)=YDATA
 5       CONTINUE
      END IF
      IF(IDATA.EQ.2) CLOSE(1)
      RETURN
 101  FORMAT(A25)
      END
***********************************************************************
*                    SUBROUTINE TO EXECUTE QR-DIVISION                *
***********************************************************************
      SUBROUTINE QR(A,R,X,EN,NA,AM,MAX,ID)
C     MODIFIED GRAM-SCHMIDT METHOD
      INTEGER AM
      DIMENSION A(MAX,MAX),R(MAX,MAX),EN(1)
      INTEGER X(1)
      DO 1 I=1,NA
      X(I)=I
      DO 1 J=1,AM
      R(I,J)=0.0
 1    CONTINUE
      DO 2 J=1,AM
         CALL NORM(A,EN,X,J,NA,AM,MAXNORM,MAX)
         IF(J.NE.AM) THEN
            ISTORE=X(J)
            X(J)=X(MAXNORM)
            X(MAXNORM)=ISTORE
            IF(J.GE.2) THEN
               IF(X(J).NE.X(IX)) THEN
                  DO 7 I=1,J-1
                     STORE=R(I,J)
                     R(I,J)=R(I,MAXNORM)
                     R(I,MAXNORM)=STORE
 7                CONTINUE
               END IF
            END IF
         END IF
         R(J,J)=EN(X(J))
         IF(R(J,J).LT.1.0E-6) THEN
             WRITE(*,*)'ランク落ちがおこりました'
             WRITE(*,*)'もう一度やり直してください'
             ID=1
             RETURN
         END IF
         DO 3 I=1,NA
            A(I,X(J))=A(I,X(J))/R(J,J)
 3       CONTINUE
         DO 4 K=J+1,AM
            DO 5 L=1,NA
               R(J,K)=R(J,K)+A(L,X(J))*A(L,X(K))
 5          CONTINUE
            DO 6 L=1,NA
               A(L,X(K))=A(L,X(K))-R(J,K)*A(L,X(J))
 6          CONTINUE
 4       CONTINUE
 2    CONTINUE
      RETURN
      END
***********************************************************************
*            SUBROUTINE TO CALCULATE THE EUCLID NORM                  *
***********************************************************************
      SUBROUTINE NORM(A,EN,X,IS,MA,NA,MAXNORM,MAX)
      DIMENSION A(MAX,MAX),EN(1)
      INTEGER X(1)
      DO 1 JJ=IS,NA
         J=X(JJ)
         EN(J)=0.0
         DO 2 I=1,MA
            EN(J)=EN(J)+A(I,J)*A(I,J)
 2       CONTINUE
         EN(J)=SQRT(EN(J))
         IF(JJ.EQ.IS) THEN
            MAXNORM=JJ
         ELSE
            IF(EN(J).GT.EN(X(MAXNORM))) THEN
               MAXNORM=JJ
            END IF
         END IF
 1    CONTINUE
      RETURN
      END
***********************************************************************
*                SUBROUTINE TO MAKE TRANSPOSED-MATRIX                 *
***********************************************************************
      SUBROUTINE TENMTX(A,B,M,N,MAX)
      DIMENSION A(MAX,MAX),B(MAX,MAX)
      DO 1 I=1,N
         DO 2 J=1,M
            B(I,J)=A(J,I)
 2       CONTINUE
 1    CONTINUE
      RETURN
      END
***********************************************************************
*        SUBROUTINE TO CALCCULATE THE PRODUCTION OF MATRIX A,B        *
***********************************************************************
      SUBROUTINE PRDMTX(A,B,C,L,M,N,MAX)
      DIMENSION A(MAX,MAX),B(MAX,MAX),C(MAX,MAX)
      DO 30 I=1,L
         DO 20 J=1,N
            C(I,J)=0.0
            DO 10 K=1,M
               C(I,J)=C(I,J)+A(I,K)*B(K,J)
 10         CONTINUE
 20      CONTINUE
 30   CONTINUE
      RETURN
      END
***********************************************************************
*               SUBROUTINE TO TRANSLATE THE MATRIX                    *
***********************************************************************
      SUBROUTINE TRANSMTX(A,B,X,M,N,MAX)
      DIMENSION A(MAX,MAX),B(MAX,MAX)
      INTEGER X(1)
      DO 1 I=1,M
         DO 2 J=1,N
            B(I,J)=A(I,X(J))
 2       CONTINUE
 1    CONTINUE
      RETURN
      END
***********************************************************************
*          SUBROUTINE TO EXECUTE THE BACKWARD SUBSTITUTION            *
***********************************************************************
      SUBROUTINE BACK(R,A,Z,M,MAX)
      DIMENSION R(MAX,MAX),A(1),Z(1)
      DO 10 I=1,M
         A(I)=0.0
 10   CONTINUE
      A(M)=Z(M)/R(M,M)
      DO 1 I=M-1,1,-1
         DO 2 J=I+1,M
            A(I)=A(I)+R(I,J)*A(J)
 2       CONTINUE
         A(I)=(Z(I)-A(I))/R(I,I)
 1    CONTINUE
      RETURN
      END
***********************************************************************
*                    SUBROUTINE TO COPY MATRIX                        *
***********************************************************************
      SUBROUTINE COPYMTX(A,B,M,N,MAX)
      DIMENSION A(MAX,MAX),B(MAX,MAX)
      DO 1 I=1,M
         DO 2 J=1,N
            B(I,J)=A(I,J)
 2       CONTINUE
 1    CONTINUE
      RETURN
      END
***********************************************************************
*                      SUBROUTINE TO SOLVE THE MLS                    *
***********************************************************************
      SUBROUTINE SOLVE(A,X,M,MENU)
      DIMENSION A(1)
      INTEGER X(1)
      WRITE(*,203)
      WRITE(*,*)
      DO 10 LOOP=1,M
         DO 20 J=1,M
            IF(X(J).EQ.LOOP) THEN
               I=J
               GOTO 30
 15         END IF
 20      CONTINUE
 30      CONTINUE
         IF(MENU.EQ.1) THEN
            WRITE(6,200) LOOP,A(I)
         ELSE IF(MENU.EQ.3) THEN
            WRITE(6,201) LOOP,A(I)
         ELSE
            WRITE(6,202) LOOP-1,A(I)
         END IF
 200     FORMAT('           X(',I3,') =',F13.3)
 201     FORMAT('          ',I3,'次の係数 =',F13.6)
 202     FORMAT('          ',I3,'次の係数 =',F13.6)
 203     FORMAT( /,10X,'*****  結 果 の 出 力  *****')
 10   CONTINUE
      RETURN
      END
***********************************************************************
*                   SUBROUTINE TO CALCULATE MENU 2                   *
***********************************************************************
      SUBROUTINE BEW(X,Y,W,SUM,Z,M)
      DIMENSION X(1),Y(1),Z(1),W(1),SUM(1)
      DO 1 I=1,M
         W(I)=1.0
         SUM(I)=W(I)
 1    CONTINUE
      DO 2 ICOUNT=1,10
         XY=0.0
         X2=0.0
         DO 3 I=1,M
            XY=XY+W(I)*X(I)*Y(I)
            X2=X2+W(I)*X(I)*X(I)
 3       CONTINUE
         V2=0.0
         A=XY/X2
         WRITE(*,101) ICOUNT,A
         DO 7 I=1,M
            WRITE(6,102) I,W(I)
 7       CONTINUE
         DO 4 I=1,M
            W(I)=ABS(Y(I)-A*X(I))
            V2=V2+W(I)*W(I)
 4       CONTINUE
         CALL MEDIAN(W,M,Z,S)
         C=20.8-15.0*EXP(-1.3730E-3*(S**1.5))
         DO 5 I=1,M
            IF(W(I).LT.S*C) THEN
               W(I)=(1-(W(I)/(C*S))**2)**2
            ELSE
               W(I)=0.0
            END IF
 5       CONTINUE
         HS=0.0
         HB=0.0
         DO 6 I=1,M
            HS=HS+ABS(W(I)-SUM(I))
            HB=HB+W(I)
            SUM(I)=W(I)
 6       CONTINUE
         IF(HS.LT.HB*1.0E-4) RETURN
 2    CONTINUE
      RETURN
 101  FORMAT(/,10X,I4,'回目の傾きの値 ＝',F13.6)
 102  FORMAT(14X,I4,'番目の重み ＝',F13.6)
      END
***********************************************************************
*               SUBROUTINE TO SEARCH FOR THE MEDIAN                   *
***********************************************************************
      SUBROUTINE MEDIAN(W,M,X,S)
      DIMENSION W(1),X(1)
      DO 1 I=1,M
         X(I)=I
 1    CONTINUE
      ID=MOD(M,2)
      IS=(M-ID)/2+1
      DO 2 I=1,IS
         DO 3 J=1,M-I
            IF(W(X(J)).LT.W(X(J+1))) THEN
               LARGE=X(J)
               X(J)=X(J+1)
               X(J+1)=LARGE
            END IF
 3       CONTINUE
 2    CONTINUE
      IF(ID.EQ.1) THEN
         S=W(X(IS))
      ELSE
         S=(W(X(IS-1))+W(X(IS)))/2.0
      END IF
      RETURN
      END
************************************************************************
*              SUBROUTINE TO CALCULATE THE NUMBER OF AIC               *
************************************************************************
      SUBROUTINE AIC(V,XS,M,N,AI)
      DIMENSION XS(1),V(1)
      PAI=3.145927
      SV=0.0
      DO 1 I=1,N
         V(I)=XS(I)-V(I)
         SV=SV+V(I)*V(I)
 1    CONTINUE
      SV=SV/REAL(N)
      IF(SV.EQ.0.0) THEN
         WRITE(*,*)
         WRITE(*,*) '   残差０です！'
         AI=-1.0E7
         GOTO 5
 2    END IF
      AI=REAL(N)*(LOG(2*PAI*SV)+1.0)+2*(REAL(M)+1.0)
  5   RETURN
      END
************************************************************************
*        SUBROUTINE TO TRANSLATE THE MATRIX ( SINGLE DIMENSION )       *
************************************************************************
      SUBROUTINE TRANS(A,V,X,M)
      DIMENSION A(1),V(1)
      INTEGER X(1)
      DO 1 LOOP=1,M
         DO 20 J=1,M
            IF(X(J).EQ.LOOP) THEN
               I=J
               GOTO 30
 15         END IF
 20      CONTINUE
 30      V(LOOP)=A(I)
 1    CONTINUE
      RETURN
      END
************************************************************************
*                    SUBROUTINE TO EXECUTE MENU 5                      *
************************************************************************
      SUBROUTINE MAKE2(X,M,N,XDATA,JCH,MAX,JCOUNT)
      DIMENSION X(MAX,MAX),XDATA(MAX,MAX)
      INTEGER JCH(1)
      IST=M
         JID=0
         JN=JCOUNT
         DO 7 I=1,IST
            IR=IST-I+1
            JCH(IR)=MOD(JN,2)
            IF(JCH(IR).EQ.1) JID=JID+1
            JN=JN/2
 7       CONTINUE
         M=JID+1
         DO 22 I=1,N
             X(I,1)=1.0
 22       CONTINUE
          ICOUNT=1
          DO 11 J=1,N
             DO 10 I=2,M
 9              CONTINUE
                IF(JCH(ICOUNT).EQ.0) THEN
                   ICOUNT=ICOUNT+1
                   GOTO 9
                END IF
                X(J,I)=XDATA(J,ICOUNT)
                ICOUNT=ICOUNT+1
 10          CONTINUE
          ICOUNT=1
 11       CONTINUE
       RETURN
 101   FORMAT(A25)
       END
************************************************************************
*                   SUBROUTINE TO SOLVE THE MENU 5                     *
************************************************************************
      SUBROUTINE SOLVE5(A,M,JCH)
      DIMENSION A(1)
      INTEGER JCH(1)
      WRITE(*,203)
      WRITE(*,*)
      JCT=0
      DO 4 I=1,M
         IF(I.EQ.1) GOTO 3
 1       IF(JCH(JCT).EQ.0) THEN
            JCT=JCT+1
            GOTO 1
 2       END IF
 3       IF(I.EQ.1) THEN
            WRITE(6,102) A(I)
         ELSE
            WRITE(6,101) JCT,A(I)
         END IF
         JCT=JCT+1
 4    CONTINUE
 101  FORMAT(  '  説明因子',I3,'番目の係数 =',F13.6)
 102  FORMAT(  19X,'定数 =',F13.6)
 203  FORMAT( /,10X,'*****  結果の出力  ****')
      RETURN
      END