      PROGRAM SNEW2
************************************************************************
*                                                                      *
*     Smoothing Newton Method for NCP (Normal Map formuration)         *
*     (CP)    F(X).GE.0  X.GE.0  F(X)^X=0.                             *
*     ---------------------------------                                *
*     N   : Number of dimensions                                       *
*     ---------------------------------                                *
*     KOUT   : Number of iterations                                    *
*     KOUT3  : Number of STEP 3                                        *
*     KOUT3A : Total number of STEP 3                                  *
*     KIN    : Number of line searches in STEP 2 or STEP 3-2           *
*     ---------------------------------                                *
*     XNEW(N)  : Current point                                         *
*     XOLD(N)  : Previous point                                        *
*     MUNEW    : Current smoothing parameter (epsilon)                 *
*     MUOLD    : Previous smoothing parameter (epsilon)                *
*     F(N)     : Mapping F(X)                                          *
*     PHI      : Smoothing function for MAX operator                   *
*     PHIMU    : Derivative of PHI w.r.t. epsiron                      *
*     PHIX     : Derivative of PHI w.r.t. x                            *
*     KAPPA    : KAPPA=(MU+1)*(MU+1)-1                                 *
*     APH(N)   : Approximate function of H                             *
*     PI       : Norm function H(X)*APH(X)/2                           *
*     THETA    : THETA=PI+KAPPA*KAPPA/2                                *
*     ---------------------------------                                *
*     Algorithm parameter                                              *
*       ETA=100                                                        *
*       LAMBDA=0.5                                                     *
*       SIGMA=0.0001                                                   *
*       ZERO=1.0D-14 : Algorithm zero                                  *
*       CC=1.0D-6    : Convergence Criterion                           *
*     ---------------------------------                                *
*                                                                      *
************************************************************************
*
      PARAMETER (N=10)
*
      REAL*8 F(N),LAMBDA,ETA,SIGMA,MUNEW,MUOLD,KAPPA,PHI,PHIMU,PHIX
      REAL*8 T(N,N),Q(N),APH(N)
      REAL*8 XNEW(N),XOLD(N),XTMP(N),YTMP(N),Z,DMU
      REAL*8 THETAN,THETAO,CC,ZERO,PIOLD,PINEW
      INTEGER KOUT,KOUT3,KOUT3A,KIN,KMAX,KOUTMAX,KMAX3,KSD
*
      CHARACTER*10 FNAME
*
      DATA ZERO,CC/1.0D-14,1.0D-6/
*
 1000 FORMAT(A)
*
*---------- Set output level
      WRITE(*,*) 'Output level (-1,0,1,2)'
      READ(*,*) PLEV
*
*---------- Set output file name
      WRITE(*,*) 'Output file name?'
      READ(*,1000) FNAME
      OPEN(7,FILE=FNAME,STATUS='UNKNOWN')
*
***********************************************************************
*                                                                     *
*    STEP 0 Initialization                                            *
*                                                                     *
***********************************************************************
*
      KOUTMAX=100
      KMAX3=1000
      KMAX=20
*
      KOUT=0
      KOUT3A=0
      KSD=0
*
*---------- Set algorithm parameter
*
      LAMBDA=0.5
      ETA=100
      SIGMA=0.0001
*
*---------- Input initial point
*
*      DO 11 I=1,N
*         WRITE(*,*) 'X(',I,')= ?'
*         READ(*,*) XOLD(I)
*   11 CONTINUE
*
         WRITE(*,*) 'X(I)= ?'
         READ(*,*) XOLD(1)
      DO 11 I=1,N
         XOLD(I)=XOLD(1)
   11 CONTINUE
*
   12 WRITE(*,*) 'epsilon_0= ?'
      READ(*,*) MUOLD
      IF(MUOLD.LE.ZERO) GOTO 12
*
*---------- Output parameters and initial points
      WRITE(7,*) 'LAMBDA=',LAMBDA
      WRITE(7,*) 'ETA=',ETA
      WRITE(7,*) 'SIGMA=',SIGMA
*
      WRITE(7,*) 'Initial point'
      DO 13 I=1,N
         WRITE(7,1010) I,XOLD(I)
   13 CONTINUE
      WRITE(7,*) 'epsilon_0=', MUOLD
*
 1010 FORMAT('X(',I3,')=',F13.6)
*
************************************************************************
*                                                                      *
*     Main routine                                                     *
*                                                                      *
************************************************************************
*
*---------- Calculate pi
*
      KAPPA=MUOLD*(MUOLD+2.0)
      PIOLD=0.0
      DO 14 I=1,N
        CALL SMOOTH(PHI,MUOLD,XOLD(I))
        XTMP(I)=PHI
   14 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 15 I=1,N
        APH(I)=F(I)+XOLD(I)-XTMP(I)
        PIOLD=APH(I)*APH(I)+PIOLD
   15 CONTINUE
      PIOLD=0.5*PIOLD
*
***********************************************************************
*                                                                     *
*    STEP 1 |min(x,F(x)|.le.CC                                        *
*                                                                     *
***********************************************************************
*
    1 CONTINUE
*
      THETAO=PIOLD+0.5*KAPPA*KAPPA
*
*---------- Convergence judgement
*
      DO 16 I=1,N
        XTMP(I)=DMAX1(XOLD(I),0D1)
   16 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 18 I=1,N
        IF(DABS(F(I)+XOLD(I)-XTMP(I)).GT.CC) GOTO 19
   18 CONTINUE
*
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
        WRITE(7,*) 'TOTAL # OF STEP 3', KOUT3A
        WRITE(7,*) '# OF SEARCH FOR SD DIRECTION', KSD
        WRITE(7,*) 'X IS A SOLUTION'
        STOP
*
*---------- OVER 100 ITERATIONS
*
   19 IF(KOUT.GE.KOUTMAX) THEN
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
        WRITE(7,*) 'TOTAL # OF STEP 3', KOUT3A
        WRITE(7,*) '# OF SEARCH FOR SD DIRECTION', KSD
        WRITE(7,*) 'OVER 100 ITERATIONS'
        STOP
      ENDIF
*
*---------- Output current point
*
      IF(PLEV.GE.0) THEN
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
      ENDIF
*
************************************************************************
*                                                                      *
*     STEP 2   T=(I-D+NABF^*D)   Q=(-APH-DMU*HMU^)                     *
*                                                                      *
************************************************************************
*
*---------- Caluculate DMU
*
      DMU=(MUOLD+2.0)/(MUOLD+1.0)
      DMU=-0.5*DMU*MUOLD
*
      CALL NABMAP(N,T,XTMP)
*
*---------- Set Q
*
      DO 20 I=1,N
        CALL SMOOTHM(PHIMU,MUOLD,XOLD(I))
        YTMP(I)=PHIMU
   20 CONTINUE
*
      DO 21 I=1,N
        Q(I)=-YTMP(I)
        DO 21 J=1,N
          Q(I)=Q(I)+T(I,J)*YTMP(J)
   21 CONTINUE
*
      DO 22 I=1,N
        Q(I)=-APH(I)-DMU*Q(I)
   22 CONTINUE
*
*---------- Set T
*
      DO 23 I=1,N
        CALL SMOOTHX(PHIX,MUOLD,XOLD(I))
        DO 24 J=1,N
           T(J,I)=PHIX*T(J,I)
   24   CONTINUE
        T(I,I)=T(I,I)-PHIX+1.0
   23 CONTINUE
*
*---------- Solve Newton equation
*
      CALL GAUSS(N,T,Q,INDEX)
      KOUT=KOUT+1
      IF(INDEX.NE.0) THEN
        KOUT3=0
        KOUT3A=KOUT3A+1
        GOTO 3
      ENDIF
      MUNEW=MUOLD+DMU
      IF(MUNEW.LT.1.0D-20) THEN
        MUNEW=1.0D-20
      ENDIF
*
      DO 25 I=1,N
        XNEW(I)=XOLD(I)+Q(I)
   25 CONTINUE
*
*---------- Output search direction
*
      IF(PLEV.GE.1) THEN
        WRITE(7,1000) 'Search direction'
        DO 26 I=1,N
          WRITE(7,1030) I,Q(I)
   26   CONTINUE
      ENDIF
*
 1030 FORMAT('D(',I3,')=',F13.6)
*
***********************************************************************
*                                                                     *
*    STEP 2 Line search                                               *
*                                                                     *
***********************************************************************
*
      KIN=0
*
*---------- Caluculate the righthand side of (16)
*
      Z=2.0*SIGMA*THETAO
*
*---------- Goto Step 3-3
*
    2 CONTINUE
      IF(KIN.GT.KMAX) THEN
        KOUT3=0
        KOUT3A=KOUT3A+1
        GOTO 3
      ENDIF
*
*---------- Caluculate theta
*
      KAPPA=MUNEW*(MUNEW+2.0)
      PIOLD=0.0
      DO 31 I=1,N
        CALL SMOOTH(PHI,MUNEW,XNEW(I))
        XTMP(I)=PHI
   31 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 32 I=1,N
        APH(I)=F(I)+XNEW(I)-XTMP(I)
        PIOLD=APH(I)*APH(I)+PIOLD
   32 CONTINUE
*
      PIOLD=0.5*PIOLD
      THETAN=PIOLD+0.5*KAPPA*KAPPA
*
*---------- Output the next point
*
      IF(PLEV.GE.2) THEN
        CALL PRT(N,F,MUNEW,XNEW,PIOLD,THETAN,KOUT,KIN)
      ENDIF
*
*---------- Convergence judgement
*
      IF(THETAO.LT.(THETAN+Z))THEN
        DO 33 I=1,N
          XNEW(I)=LAMBDA*XNEW(I)+(1-LAMBDA)*XOLD(I)
   33   CONTINUE
          MUNEW=LAMBDA*MUNEW+(1-LAMBDA)*MUOLD
        Z=LAMBDA*Z
        KIN=KIN+1
        GOTO 2
      ENDIF
*
      DO 35 I=1,N
        XOLD(I)=XNEW(I)
   35 CONTINUE
      MUOLD=MUNEW
*
      IF(KIN.EQ.0) GOTO 1
*
***********************************************************************
*                                                                     *
*    STEP 3-1 Convergence judgement                                   *
*                                                                     *
***********************************************************************
*
*---------- Convergence judgement
*
      KOUT3=0
*
    6 IF(PLEV.GE.1) THEN
        WRITE(7,*) 'TOTAL # OF STEP 3', KOUT3
        WRITE(7,*) '# OF SEARCH FOR SD DIRECTION', KSD
        WRITE(7,*) 'PI=', PIOLD
        WRITE(7,*) 'epsilon=', MUOLD
        WRITE(7,*) 'eta*epsilon=', ETA*MUOLD
      ENDIF
*
      IF(PIOLD.LE.(ETA*MUOLD)) GOTO 1
*
************************************************************************
*                                                                      *
*     STEP 3-2 Solve (17)  T=(I-D+NABF^*D)   Q=(-APH)                  *
*                                                                      *
************************************************************************
*
*---------- Over 1000 iterations at STEP 3
*
      IF(KOUT3.GE.KMAX3) THEN
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
        WRITE(7,*) 'OVER 1000 ITERATIONS AT STEP 3'
        WRITE(7,*) 'TOTAL # OF STEP 3', KOUT3A
        WRITE(7,*) 'TOTAL # OF SD', KSD
        STOP
      ENDIF
*
*---------- Output current point
*
      IF(PLEV.GE.0) THEN
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
      ENDIF
*
*---------- Set Q
*
      DO 45 I=1,N
        Q(I)=-APH(I)
   45 CONTINUE
*
*---------- Set T
*
      CALL NABMAP(N,T,XTMP)
      DO 47 I=1,N
        CALL SMOOTHX(PHIX,MUOLD,XOLD(I))
        DO 48 J=1,N
           T(J,I)=PHIX*T(J,I)
   48   CONTINUE
        T(I,I)=T(I,I)-PHIX+1.0
   47 CONTINUE
*
*---------- Solve (17)
*
      CALL GAUSS(N,T,Q,INDEX)
      KOUT3=KOUT3+1
      KOUT3A=KOUT3A+1
      IF(INDEX.NE.0) GOTO 3
      DO 49 I=1,N
        XNEW(I)=XOLD(I)+Q(I)
   49 CONTINUE
*
*---------- Output search direction
*
      IF(PLEV.GE.1) THEN
        WRITE(7,1000) 'Search direction'
        DO 46 I=1,N
          WRITE(7,1030) I,Q(I)
   46   CONTINUE
      ENDIF
*
*---------- Line search
*
      KIN=0
*
*---------- Calculate the right hand side of (18)
*
      Z=2.0*SIGMA*PIOLD
*
*---------- Over iterations
*
    8 CONTINUE
      IF(KIN.GT.KMAX) GOTO 3
*
*---------- Caluculate pi
*
      PINEW=0.0
      DO 51 I=1,N
        CALL SMOOTH(PHI,MUOLD,XNEW(I))
        XTMP(I)=PHI
   51 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 52 I=1,N
        APH(I)=F(I)+XNEW(I)-XTMP(I)
        PINEW=APH(I)*APH(I)+PINEW
   52 CONTINUE
*
      PINEW=0.5*PINEW
*
*---------- Output the next point
*
      IF(PLEV.GE.2) THEN
        CALL PRT(N,F,MUOLD,XNEW,PINEW,THETAO,KOUT,KIN)
      ENDIF
*
*---------- Convergence judgement
*
      IF(PIOLD.LT.(PINEW+Z))THEN
        DO 53 I=1,N
          XNEW(I)=LAMBDA*XNEW(I)+(1-LAMBDA)*XOLD(I)
   53   CONTINUE
        Z=LAMBDA*Z
        KIN=KIN+1
        GOTO 8
      ENDIF
*
      DO 55 I=1,N
        XOLD(I)=XNEW(I)
   55 CONTINuE
      PIOLD=PINEW
      GOTO 6
*
************************************************************************
*                                                                      *
*     STEP 3-3 Search for SD                                           *
*                                                                      *
************************************************************************
*
*---------- Output Current point
*
    3 IF(PLEV.GE.0) THEN
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
      ENDIF
*
      KSD=KSD+1
*
      PIOLD=0.0
      DO 61 I=1,N
        CALL SMOOTH(PHI,MUOLD,XOLD(I))
        XTMP(I)=PHI
   61 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 62 I=1,N
        APH(I)=F(I)+XOLD(I)-XTMP(I)
        PIOLD=APH(I)*APH(I)+PIOLD
   62 CONTINUE
      PIOLD=0.5*PIOLD
*
*---------- Set T
*
      CALL NABMAP(N,T,XTMP)
      DO 63 I=1,N
        CALL SMOOTHX(PHIX,MUOLD,XOLD(I))
        DO 64 J=1,N
           T(J,I)=PHIX*T(J,I)
   64   CONTINUE
        T(I,I)=T(I,I)-PHIX+1.0
   63 CONTINUE
*
*---------- Search direction
*
      DO 65 I=1,N
        Q(I)=0.0
        DO 65 J=1,N
           Q(I)=Q(I)-T(J,I)*APH(J)
   65 CONTINUE
*
      DO 66 I=1,N
        XNEW(I)=XOLD(I)+Q(I)
   66 CONTINUE
*
*---------- Output search direction
*
      IF(PLEV.GE.1) THEN
        WRITE(7,1000) 'Search direction'
        DO 67 I=1,N
          WRITE(7,1030) I,Q(I)
   67   CONTINUE
      ENDIF
*
*---------- Line search
*
      KIN=0
*
*---------- Over iterations
*
    9 CONTINUE
      IF(KIN.GT.KOUTMAX) THEN
        WRITE(7,*) 'Over iterations at STEP 3'
        CALL MAP(N,F,XOLD)
        CALL PRT(N,F,MUOLD,XOLD,PIOLD,THETAO,KOUT,KIN)
        WRITE(7,*) '# of STEP 3', KOUT3
        WRITE(7,*) 'Total # of STEP 3', KOUT3A
        WRITE(7,*) 'Total # of SD', KSD
        STOP
      ENDIF
*
*---------- Calculate pi
*
      PINEW=0.0
      DO 68 I=1,N
        CALL SMOOTH(PHI,MUOLD,XNEW(I))
        XTMP(I)=PHI
   68 CONTINUE
*
      CALL MAP(N,F,XTMP)
      DO 69 I=1,N
        APH(I)=F(I)+XNEW(I)-XTMP(I)
        PINEW=APH(I)*APH(I)+PINEW
   69 CONTINUE
*
      PINEW=0.5*PINEW
*
*---------- Output the next point
*
      IF(PLEV.GE.2) THEN
        CALL PRT(N,F,MUOLD,XNEW,PINEW,THETAO,KOUT,KIN)
      ENDIF
*
*---------- Convergence Judgement
*
      IF(PIOLD.LE.PINEW)THEN
        DO 70 I=1,N
          XNEW(I)=LAMBDA*XNEW(I)+(1-LAMBDA)*XOLD(I)
   70   CONTINUE
        KIN=KIN+1
        GOTO 9
      ENDIF
*
      DO 71 I=1,N
        XOLD(I)=XNEW(I)
   71 CONTINUE
      PIOLD=PINEW
      GOTO 6
*
      END
*
