         pROGRAM Ala22cgGS3  
C CCC for ALA potential N=66 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.09.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c test for script version - modular structure 
c GS cycle for RP: NewPoint calculation by projected gradient
c  goes along reduced gradient to chain point
c
c  corrector step  by (modified by wq) CG+ of Nocedal et al.  
C  the actual point, to start the correction, is Y(LA,..)
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      Integer istatus,ITALL,Natoms(22),N3 
      Integer LA,NDIM,NDIM1,L  
      PARAMETER(NDIM=66,L=12,NDIM1=67)
      DOUBLE PRECISION XX(NDIM1),G(NDIM1),DDOT,XTOL 
      double precision d(ndim1),gold(ndim1),w(ndim1)
      double precision Fun,EPS,TLEV,X(NDIM),Proj(NDIM)
c        double precision time1,time2,tottime   
      REAL*8 Y(L+1,NDIM),VV(NDIM),D1,Yold(NDIM) 
      REAL*8 Rs(NDIM),energy,Cx,gnorm,Pena
      REAL*8 ProMat(NDIM,NDIM)
      logical      finish
      integer iprint(2),iflag,icall,n,method,iter,nfun,IREST
      Integer JJ,I,J,Np,Np1,IterMax
      Character*1 CharAtoms(22)
      COMMON /RUNINF/ITER,NFUN
      Data CharAtoms/'h','c','h','c','h','o','n','c','h','c','h','c',
     1               'o','n','c','h','h','h','h','h','h','h'/
      Data Natoms/1,6,1,6,1,8,7,6,1,6,1,6,8,7,6,1,1,1,1,1,1,1/
c
       OPEN(9,FILE='ala22chain.dat')  
       OPEN(10,FILE='ala22Rs.dat')
       OPEN(11,FILE='ala22param.dat') 
       OPEN(12,FILE='ala22proj.dat') 
       OPEN(13,FILE='ala22point.dat') 
       OPEN(14,FILE='ala22ener.dat') 
       OPEN(15,FILE='ala22grad.dat')
       OPEN(16,FILE='ala22engs.weg')            
       OPEN(17,FILE='ala22lamb.dat')
       OPEN(19,FILE='ala22grno.dat')
       OPEN(21,FILE='iflag.dat')
       OPEN(22,FILE='cgplFi.dat')
       OPEN(23,FILE='cggsW.dat')
       OPEN(25,FILE='ala22dddd.dat') 
       OPEN(26,FILE='ala22gold.dat') 
       OPEN(28,FILE='ala22yold.dat') 
       OPEN(32,FILE='alaatoms.dat',status='old')
       OPEN(33,FILE='alapoig.dat')
       OPEN(38,FILE='ala22chain.mol',status='unknown',access='append')
       OPEN(44,FILE='protocol.txt',status='unknown',access='append')
c
c   Read problem input information
      Pena= 1.0000001d+0 
      N  = NDIM
      N1 = N+1
      N3 = N/3

C N1  for Lagrange formulation of one additional constraint
C N3: number of atoms
      IPRINT(1)= 0
      IPRINT(2)= 0   
      XTOL= 11.11111D-11
      ITERMAX= 2*N   
c Note: EPS1= 1.0D-5 of CG+ is very to small here !!!
C CG+ 
      method =    3
      IREST  =    1
      rewind 21
      read(21,*) IFLAG
      write( 6,*) '        iflag is   ',  IFLAG
c
      if(IFLAG.ne.0)  then
         rewind 22
         read(22,*) FINISH,ICALL,NFUN
         rewind 23
         read(23,50) (W(J),J=1,N) 
         rewind 25
         read(25,50) (D(J),J=1,N1)
         rewind 26
         read(26,50) (Gold(J),J=1,N1)
         rewind 28
         read(28,50) (Yold(J),J=1,N)
         rewind 17
         read(17,50) XX(N1)
       endif 
c
       istatus=0
       rewind 11
       read(11,*) JJ, LA, eps ,ITALL, NFUN
cc       itGS=itGS+1 -->>  other count: by NFUN of CG+
cc       JJ=L    !! 
        rewind 11
        write(11,*) L, LA, eps ,ITALL, NFUN
        write(44,*)' L=',L,' LA= ',LA,' ITall = ',ITALL,'+',NFUN
c Print parameters
      if (iprint(1).ge.0 .AND. LA.lt.3 .AND. IFLAG.eq.0)then
         write(*,820)
         write(*,840) N1, method, irest
         write(44,820)
         write(44,840) N1, method, irest
      endif
c 20    CONTINUE : NOTE: cycle replaced by script organisation
        If(NFUN .gt. 2*N1) then
         write(6,*) ' Cancel Loop: ERROR , no Convergence at '
         write(6,*) ' LA=  ' ,LA, ' NFUN  =  ', NFUN
         write(44,*)' Cancel Loop: ERROR , no Convergence at '
         write(44,*)' LA=  ' ,LA, ' NFUN  =  ', NFUN
         istatus=1
         iflag=0
         ITALL = ITALL+NFUN
         NFUN=0 
         goto 111
        Endif
c Call the function and gradient values here
       rewind 10
       READ(10,50) (Rs(I),I=1,N)
       rewind 12
       Do 3 J=1,N
3      READ(12,50) (ProMat(J,I),I=1,N)
       rewind 13
       READ(13,*) (X(I),I=1,N)
       rewind 14
       READ(14,*) energy
       WRITE(44,*)'energy=',energy 
       rewind 9
       Do 4, J=1,L+1
4     Read(9,50) (Y(J,I),I=1,N)
       rewind 15
ccc       READ(15,50) (VV(I),I=1,N)
       Do 1, J=1,N3 
1       READ(15,*,END=111) J1, JJ, (VV( 3*(J-1)+I ),I=1,3)
cc g3       call PERROR("I/O processing")
cc g3       write(6,*) ' for gradient '
       IF(J1 .ne. N3 ) WRITE(44,*) ' ERROR in gradient reading !!! '
       Do 5 I=1,N
5      VV(I)=-VV(I)
c using the special file with Chartesian Grad output
c from Gaussian03 !! note: - sign from gaussian gradient for VV !!
c
       if(IFLAG.eq.0)  then
         Do 211 I=1,N
211      Yold(I)= Y(LA,I)
         rewind 28
         write(28,50) (Yold(I),I=1,N)
       endif  
c   norm of gradient
       gnorm= DSQRT(DDOT(N,VV,1,VV,1))
       WRITE(6,*) '  PES full gradient norm ',gnorm
       WRITE(44,*) '  PES full gradient norm ',gnorm
c
      Do 13 J=1,N
        Proj(J)=VV(J)
13      XX(J)=X(J)   
cc      if(iflag.eq.0) :    no , do not use this!
cc   enforce before every CG+ step the "true" lambda   
        XX(N1)= gnorm/Pena
c estimation for lambda of Lagrange ansatz (with + sign):
c turn direction of RS into direction of gradient: 
        If(DDOT(N,Rs,1,VV,1) .lt. 0.0d0) then 
         Do 14 J=1,N  
14        Rs(J)= -Rs(J)
        endif 
cccccccccccccccccccccccccccccccccccccccc
C   Newton trajectory with RGF Projector
c           and  Projector of rotations:   
C   Projektion of Gradient    
cccccccccccccccccccccccccccccccccccccccc
       call PROJECTION(N,Proj,ProMat)

       D1= DSQRT(DDOT(N,Proj,1,Proj,1))
       WRITE(6,*) ' projected   grad   norm ', D1
       WRITE(44,*)'  projected  grad   norm ', D1

cc     WRITE(44,*)'  lambda guess           ', XX(N1)
cc     WRITE(44,*)'  Penalty factor         ', Pena
      If(D1 .lt. eps) THEN
        WRITE(44,*)'at IteGS=',NFUN,'  ok - proj norm less then', 
     1   '  tolerance ', eps
        WRITE(44,*)' successful convergence (no ERROR) '
        WRITE(6, *)' IteGS=',NFUN,' == ok - proj norm ',D1
        ITALL = ITALL+NFUN
        NFUN=0 
        IFLAG=0 
        istatus=1
        GOTO 111
      ENDIF
c
c CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+ CG+
c  program for running the conjugate gradient methods described in
c  Gilbert, J.C. and Nocedal, J. (1992).
c  "Global Convergence Properties of Conjugate Gradient Methods"
c  SIAM Journal on Optimization, Vol.2, pp. 21-42.
c  Written by G. Liu, J. Nocedal and R. Waltz, October 1998
c  IFLAG=0 indicates an initial entry
c wq:
C   one constraint: C(x)=0= Sum_i^n rs(i) * (x(i)-Yold(i) )
C   Lagrange function:  F(x,Lambda)=fcn(x)-Lambda*c(x)
C   where Lambda = xx(n+1) is the Lagrange parameter  
c   Pena is an additional penalty factor, put to 1.0
cc
        Cx=0.0d0
        Do 23 J=1,N
23     Cx= Cx + (X(J)-Yold(J))*Rs(J)*Pena
        write(44,*) ' C(x) = ' , Cx
c   Energy and Gradient for constrained problem:
cc 
        Fun= energy - XX(N1)*Cx
        Do 18 J=1,N
18      G(J)= VV(J) - XX(N1)*Rs(J)*Pena
        G(N1)=  -Cx
c  TESTs
c        write(44,*) ' G(J)   = ' , (G(J),J=1,N1)
c          D1=0.0d0 
c          Do 1388 J=1,N
c1388      D1= D1+ G(J)*Rs(J)
c         write(44,*) 'sum G(J)*Rs(J) = ' ,D1
c          D1=0.0d0 
c          Do 1381 J=1,N
c1381      D1= D1+ VV(J)*Rs(J)
c         write(44,*) 'sum VV(J)*Rs(J)= ' ,D1
c
30      CONTINUE
c   Call the main optimization code
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        CALL CGFAM(N1,XX,Fun,G,D,GOLD,IPRINT,W,
     *            IFLAG,IREST,METHOD,FINISH,ProMat,Rs)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       IFLAG=
c              0 : successful termination
c              1 : return to evaluate F and G
c              2 : return with a new iterate, try termination test
cccc          -i : error
c        write(44,*) ' NFUN= ', NFUN,' at IFLAG',IFLAG,
c     1              ': lambda given by CG+ ',XX(N1) 
        rewind 11
        write(11,*) L, LA, eps ,ITALL, NFUN
        Do 201,J=1,N
 201    X(J) = XX(J) 
        if(IFLAG.ne.2)  write(44,*)  
c actual point approximation:
        rewind 13
        WRITE(13,50) (XX(I),I=1,N)
        rewind 17
        WRITE(17,50)  XX(N1)  

ccc   for gaussian input file
c        rewind 32
c        READ(32,*) (Natoms(I),I=1,N3)
        rewind 33       
        Do 7, J=1,N3    
 7      WRITE(33,*) Natoms(J), (X(3*(J-1)+I ),I=1,3)
        WRITE(33,*) '  ' 
c  store CG+ data for script organisation:
         rewind 23
         write(23,50) (W(J),J=1,N1)
         rewind 25
         write(25,50) (D(J),J=1,N1)
         rewind 26
         write(26,50) (Gold(J),J=1,N1)
         rewind 21
         write(21,*) IFLAG
cc       write(44,*)  '  iflag nach CG+ ',  IFLAG
         write( 6,*)  '  iflag nach CG+ ',  IFLAG
         rewind 22
         write(22,*) FINISH,ICALL,NFUN

      IF(IFLAG.LE.0.OR.ICALL.GT.ITERMAX) GO TO 501

       IF(IFLAG.EQ.1) THEN
         ICALL=ICALL + 1
         WRITE(6,*) '   Line search LOOP,   ICALL=  ',ICALL
         rewind 22
         write(22,*) FINISH,ICALL,NFUN
         rewind 11
         write(11,*) L, LA, eps ,ITALL, NFUN
            istatus=0
            call exit(istatus)
            STOP
cc         GO TO 20
        ENDIF    

        IF(IFLAG.EQ.2) THEN
c Termination Test. The user may replace it by some other test;
c the parameter 'FINISH' must be set to 'TRUE' when the test is satisfied.
         TLEV= EPS
cc wq              *(1.0d0 + DABS(F))
          I=0
 41       I=I+1
          IF(I.GT.N) THEN
cc        IF(I.GT.N1) THEN  
cc wq 
             IF(DABS(G(N1))/Pena.GT.TLEV) go to 30
             FINISH = .TRUE.
             write(44,*) '            finish = true   '
             ITALL=ITALL+NFUN
             NFUN=0 
             GO TO 30
           ENDIF
cc wq         IF(DABS(G(I)).GT.TLEV) THEN
          Do 2020,J=1,N
2020      Proj(J) = G(J)
          call PROJECTION(N,Proj,ProMat)
          IF(DABS(Proj(I)).GT.TLEV) THEN 
            GO TO 30
          ELSE
            GO TO 41
          ENDIF
        ENDIF
c 
 501  continue

      IF(ICALL.eq.ITERMAX) then
        istatus=1
        ITALL=ITALL+NFUN
        NFUN=0 
        write(44,*)'Stop Corr: ICALL=ITERMAX under IFLAG=',IFLAG
       endif
      IF(IFLAG.eq.0) then
        write(44,*)'ICALL= ',ICALL,' under succsessful IFLAG=',IFLAG
      else
        write(44,*)'ERROR ? ICALL= ',ICALL,' under IFLAG=',IFLAG
      endif

      rewind 11
      write(11,*) L, LA, eps ,ITALL, NFUN
      write(44,*)' L=',  L,' LA= ', LA ,' ITall = ',ITALL, ' + ',NFUN
c
      if (iprint(1).ge.0.and.iflag.ge.0) then
         write (*,890) Fun
         write (*,*) ' at LA    ', LA
      end if
           write(44,890) energy
           write(44,50)  (x(i),i=1,N)
           write(44,*) ' C(x)= ', Cx/Pena
c           write(44,*) 'lambda=?gnorm? = ', XX(N1)*Pena
c     *  ,'            at LA    ', LA  
           write(44,*) 
111    continue
        IF(LA.le.L) write(44,*) ' new step: '
          write(38,*) '    22   '
          write(38,*)
c          rewind 32
c          READ(32,*) (Natoms(I),I=1,N3)
          Do 16, J=1,N3    
16         WRITE(38,*) CharAtoms(J),'  ',(X(3*(J-1)+I ),I=1,3)
cccc
        If(IFLAG.ne.1) then
          Do 17, J=1,N3    
17         WRITE(44,*) CharAtoms(J),'  ',(X(3*(J-1)+I ),I=1,3)
        endif

      IF(IFLAG.eq.0) then
        istatus=1
        ITALL = ITALL+NFUN
        NFUN=0 
      endif   

      IF(IFLAG.lt.0) then
       write(6,*)'  Cancel Loop: no Convergence in CG+ '
c       write(44,*)'ERROR : Corr no Convergence in CG+ '
       rewind 11
       write(11,*) L, LA, eps ,ITALL, NFUN
       write(44,*)' L=',  L,' LA= ', LA ,' ITall = ',ITALL, ' + ',NFUN
       istatus=1
       ITALL = ITALL+NFUN
       NFUN=0 
       call exit(istatus)
       STOP
      endif   
c
       rewind 21
       write(21,*) iflag
c   Istatus controls the script loops 
      If(LA.gt.L) then
            ITall = ITall +NFUN
            NFUN=0
            istatus=10 
      endif
c  new chain point becomes:
        IF(finish) then
          Do 117 I=1,N
117       Y(LA,I)= XX(I)
          rewind 9
          Do 118, J=1,L+1
118       write(9,50) (Y(J,I),I=1,N)

          WRITE(6,*)' energy= ',energy
          WRITE(16,*) '{ ' 
          IF(LA.le.L)  WRITE(16,*) '  ',  energy,' , '
          IF(LA.eq.L+1)WRITE(16,*) '  ',  energy,' } '
          istatus=1
          call exit(istatus)

          STOP
        endif
       rewind 11    
       write(11,*) L, LA, eps ,ITALL, NFUN
       write(44,*)' at the end of frame of CGFAM '
       write(44,*)' L=',  L,' LA= ', LA ,' ITall = ',ITALL, ' + ',NFUN
       
       call exit(istatus)
       STOP
 50    FORMAT(1X,3F20.10)
 820  format ( ' Conjugate Gradient Minimization Routine')
 840  format ( ' n   =', i6, 
     *         '     method =', i6,
     *         '     irest  =', i6)
 850  format (/,'  Error: negative N value')
 860  format (/,'  Error: N too large, increase parameter ndim')
 890  format (/,' f(x*) =', 1pd16.8)
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       SUBROUTINE PROJECTION(N,VEC,ProMat)
c  projection of vector VEC of dim N 
c  by Projection Matrix ProMat of dim(N,N)
      PARAMETER(Ndim=66)
      REAL*8 VEC(N),ProMat(N,N),Proj(Ndim)
      INTEGER N,J,JJ
      Do 1 J=1,N
      Proj(J)=0.0d0
      Do 2 JJ=1,N   
2     Proj(J)=Proj(J) + Vec(JJ)*ProMat(J,JJ)
1     Continue
      Do 3 J=1,N 
3     Vec(J) = Proj(J)     
      RETURN
      END  
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
      SUBROUTINE CVSMOD(N,X,F,G,S,STP,FTOL,GTOL,XTOL,
     *           STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA,DGINIT,DGOUT)
      INTEGER N,MAXFEV,INFO,NFEV
      DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
      DOUBLE PRECISION X(N),G(N),S(N),WA(N)
C
C     *** This is a modification of More's line search routine **
C                   * * * * * * 
C     THE PURPOSE OF CVSMOD IS TO FIND A STEP WHICH SATISFIES
C     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
C     THE USER MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE
C     FUNCTION AND THE GRADIENT.
C
C     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
C     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
C     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
C     MINIMIZER OF THE MODIFIED FUNCTION
C
C          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
C
C     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
C     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
C     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
C     CONTAINS A MINIMIZER OF F(X+STP*S).
C
C     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
C     THE SUFFICIENT DECREASE CONDITION
C
C           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
C
C     AND THE CURVATURE CONDITION
C
C           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
C
C     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
C     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
C     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
C     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
C     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
C     SATISFIES THE SUFFICIENT DECREASE CONDITION.
C
C     THE SUBROUTINE STATEMENT IS
C
C        SUBROUTINE CVSMOD(N,X,F,G,S,STP,FTOL,GTOL,XTOL,
C                   STPMIN,STPMAX,MAXFEV,INFO,NFEV,WA,DG,DGOUT)
C     WHERE
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF VARIABLES.
C       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
C         X + STP*S.
C       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
C         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
C       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
C         OF F AT X + STP*S.
C       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
C         SEARCH DIRECTION.
C       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
C         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
C         STP CONTAINS THE FINAL ESTIMATE.
C       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION
C         OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE
C         DIRECTIONAL DERIVATIVE CONDITION ARE SATISFIED.
C       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
C         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C         IS AT MOST XTOL.
C       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
C         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP.
C       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
C         MAXFEV BY THE END OF AN ITERATION.
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C
C         INFO = 0  IMPROPER INPUT PARAMETERS.
C         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
C         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
C                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
C         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C                   IS AT MOST XTOL.
C         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
C         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
C         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
C         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
C                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
C                   TOLERANCES MAY BE TOO SMALL.
C
C       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C         CALLS TO FCN.
C       WA IS A WORK ARRAY OF LENGTH N.
C
C       *** The following two parameters are a modification to the code
C       DG IS THE INITIAL DIRECTIONAL DERIVATIVE (IN THE ORIGINAL CODE
C                 IT WAS COMPUTED IN THIS ROUTINE0
C       DGOUT IS THE VALUE OF THE DIRECTIONAL DERIVATIVE WHEN THE WOLFE
C             CONDITIONS HOLD, AND AN EXIT IS MADE TO CHECK DESCENT.
C
C     SUBPROGRAMS CALLED
C       CSTEPM 
C       FORTRAN-SUPPLIED...ABS,MAX,MIN
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C
      INTEGER INFOC,J
      LOGICAL BRACKT,STAGE1
      DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO,DGOUT
      DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
      save
      OPEN(44,FILE='protocol.txt',status='unknown',access='append')
      OPEN(31,FILE='cvsRest.dat')
C wq    for script organisation 
       IF(INFO.ne.0) then
         rewind31 
         read(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT,
     *       INFOC, BRACKT,STAGE1
       endif
C end wq
      IF(INFO.EQ.-1) GOTO 45
      IF(INFO.EQ.1)  GOTO 321
      INFOC = 1
C wq
      if(STP .LT. ZERO) then
       write(44,*) 'ERROR - warning: STP < 0 ??? '
       goto 1    
      endif
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
     *    GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
     *    .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) Then
        write(44,*)  '  return out off CVSMOD by ERROR   '  
        RETURN
      endif 

1     continue 

C
C     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
C     AND CHECK THAT S IS A DESCENT DIRECTION.
C
      IF (DGINIT .GE. ZERO) then

c         rewind31
c         write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
c     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
c     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT,  
c     *       INFOC, BRACKT,STAGE1
c        RETURN
         write(44,*) 'ERROR - warning: DGINIT.GE.ZERO 0 ??? '
      endif
C
C     INITIALIZE LOCAL VARIABLES.
C
      BRACKT = .FALSE.
      STAGE1 = .TRUE.
      NFEV = 0
      FINIT = F
      DGTEST = FTOL*DGINIT
      WIDTH = STPMAX - STPMIN
      WIDTH1 = WIDTH/P5
      DO 20 J = 1, N
         WA(J) = X(J)
   20    CONTINUE
C
C     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
C     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
C     THE INTERVAL OF UNCERTAINTY.
C     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
C     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
C
      STX = ZERO
      FX = FINIT
      DGX = DGINIT
      STY = ZERO
      FY = FINIT
      DGY = DGINIT
C
C     START OF ITERATION.
C
   30 CONTINUE
C
C        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
C        TO THE PRESENT INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            STMIN = MIN(STX,STY)
            STMAX = MAX(STX,STY)
         ELSE
            STMIN = STX
            STMAX = STP + XTRAPF*(STP - STX)
            END IF
C
C        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
C
         STP = MAX(STP,STPMIN)
         STP = MIN(STP,STPMAX)
C
C        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
C        STP BE THE LOWEST POINT OBTAINED SO FAR.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
     *      .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
C
C        EVALUATE THE FUNCTION AND GRADIENT AT STP
C        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
C
         DO 40 J = 1, N
            X(J) = WA(J) + STP*S(J)
   40       CONTINUE
C        Return to compute function value
         INFO=-1
         rewind31
         write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, 
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT,  
     *       INFOC, BRACKT,STAGE1
         RETURN
C
   45    INFO=0
         NFEV = NFEV + 1
         DG = ZERO
         DO 50 J = 1, N
            DG = DG + G(J)*S(J)
   50       CONTINUE
         FTEST1 = FINIT + STP*DGTEST
C
C        TEST FOR CONVERGENCE.
C
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
     *      .OR. INFOC .EQ. 0) INFO = 6
         IF (STP .EQ. STPMAX .AND.
     *       F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
         IF (STP .EQ. STPMIN .AND.
     *       (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
         IF (NFEV .GE. MAXFEV) INFO = 3
         IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
C        More's code has been modified so that at least one new
C        function value is computed during the line search (enforcing
C        at least one interpolation is not easy, since the code may
C        override an interpolation)
         IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)
     *       .AND. NFEV.GT.1) INFO = 1
C
C        CHECK FOR TERMINATION.
C
         IF (INFO .NE. 0)THEN
            DGOUT=DG
         rewind31
         write(31,*)DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
     *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, 
     *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,DGOUT,  
     *       INFOC, BRACKT,STAGE1
            RETURN
         ENDIF
 321     continue
C
C        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
C
         IF (STAGE1 .AND. F .LE. FTEST1 .AND.
     *       DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
C
C        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
C        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
C        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
C        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
C        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
C
         IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
C
C           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
C
            FM = F - STP*DGTEST
            FXM = FX - STX*DGTEST
            FYM = FY - STY*DGTEST
            DGM = DG - DGTEST
            DGXM = DGX - DGTEST
            DGYM = DGY - DGTEST
C
C           CALL CSTEPM TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL CSTEPM(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
     *                 BRACKT,STMIN,STMAX,INFOC)
C
C           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
C
            FX = FXM + STX*DGTEST
            FY = FYM + STY*DGTEST
            DGX = DGXM + DGTEST
            DGY = DGYM + DGTEST
         ELSE
C
C           CALL CSTEPM TO UPDATE THE INTERVAL OF UNCERTAINTY
C           AND TO COMPUTE THE NEW STEP.
C
            CALL CSTEPM(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
     *                 BRACKT,STMIN,STMAX,INFOC)
         END IF
C
C        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
C        INTERVAL OF UNCERTAINTY.
C
         IF (BRACKT) THEN
            IF (ABS(STY-STX) .GE. P66*WIDTH1)
     *         STP = STX + P5*(STY - STX)
            WIDTH1 = WIDTH
            WIDTH = ABS(STY-STX)
            END IF
C
C        END OF ITERATION.
C
         GO TO 30
C
C     LAST CARD OF SUBROUTINE CVSMOD.
C
      END
C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE CSTEPM(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
     *                 STPMIN,STPMAX,INFO)
      INTEGER INFO
      DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
      LOGICAL BRACKT,BOUND
C     **********
C
C     SUBROUTINE CSTEPM
C
C     THE PURPOSE OF CSTEPM IS TO COMPUTE A SAFEGUARDED STEP FOR
C     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
C     A MINIMIZER OF THE FUNCTION.
C
C     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
C     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
C     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
C     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
C     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
C     WITH ENDPOINTS STX AND STY.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE CSTEPM(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
C                        STPMIN,STPMAX,INFO)
C
C     WHERE
C
C       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
C         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
C         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
C         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
C
C       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
C         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
C         UPDATED APPROPRIATELY.
C
C       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
C         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
C         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
C         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
C
C       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
C         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
C         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
C         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
C
C       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
C         AND UPPER BOUNDS FOR THE STEP.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C         IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
C         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
C         INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
C
C     FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C     JORGE J. MORE', DAVID J. THUENTE
C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
      DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
      INFO = 0
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.
     *    STP .GE. MAX(STX,STY))) .OR.
     *    DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
C
C     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
C
      SGND = DP*(DX/ABS(DX))
C
C     FIRST CASE. A HIGHER FUNCTION VALUE.
C     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
C     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
C     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
C
      IF (FP .GT. FX) THEN
         INFO = 1
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .LT. STX) GAMMA = -GAMMA
         P = (GAMMA - DX) + THETA
         Q = ((GAMMA - DX) + GAMMA) + DP
         R = P/Q
         STPC = STX + R*(STP - STX)
         STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX)
         IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN
            STPF = STPC
         ELSE
           STPF = STPC + (STPQ - STPC)/2
           END IF
         BRACKT = .TRUE.
C
C     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
C     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
C     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
C     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
C
      ELSE IF (SGND .LT. 0.0) THEN
         INFO = 2
         BOUND = .FALSE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = ((GAMMA - DP) + GAMMA) + DX
         R = P/Q
         STPC = STP + R*(STX - STP)
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN
            STPF = STPC
         ELSE
            STPF = STPQ
            END IF
         BRACKT = .TRUE.
C
C     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
C     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
C     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
C     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
C     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
C     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
C     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
C
      ELSE IF (ABS(DP) .LT. ABS(DX)) THEN
         INFO = 3
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
C
C        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
C        TO INFINITY IN THE DIRECTION OF THE STEP.
C
         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = (GAMMA + (DX - DP)) + GAMMA
         R = P/Q
         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
            STPC = STP + R*(STX - STP)
         ELSE IF (STP .GT. STX) THEN
            STPC = STPMAX
         ELSE
            STPC = STPMIN
            END IF
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (BRACKT) THEN
            IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
         ELSE
            IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
            END IF
C
C     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
C     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
C     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
C     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
C
      ELSE
         INFO = 4
         BOUND = .FALSE.
         IF (BRACKT) THEN
            THETA = 3*(FP - FY)/(STY - STP) + DY + DP
            S = MAX(ABS(THETA),ABS(DY),ABS(DP))
            GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S))
            IF (STP .GT. STY) GAMMA = -GAMMA
            P = (GAMMA - DP) + THETA
            Q = ((GAMMA - DP) + GAMMA) + DY
            R = P/Q
            STPC = STP + R*(STY - STP)
            STPF = STPC
         ELSE IF (STP .GT. STX) THEN
            STPF = STPMAX
         ELSE
            STPF = STPMIN
            END IF
         END IF
C
C     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
C     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
C
      IF (FP .GT. FX) THEN
         STY = STP
         FY = FP
         DY = DP
      ELSE
         IF (SGND .LT. 0.0) THEN
            STY = STX
            FY = FX
            DY = DX
            END IF
         STX = STP
         FX = FP
         DX = DP
         END IF
C
C     COMPUTE THE NEW STEP AND SAFEGUARD IT.
C
      STPF = MIN(STPMAX,STPF)
      STPF = MAX(STPMIN,STPF)
      STP = STPF
      IF (BRACKT .AND. BOUND) THEN
         IF (STY .GT. STX) THEN
            STP = MIN(STX+0.66*(STY-STX),STP)
         ELSE
            STP = MAX(STX+0.66*(STY-STX),STP)
            END IF
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE CSTEPM.
      END

C   --------------------------------------------------------------------
C   Conjugate Gradient methods for solving unconstrained nonlinear
C   optimization problems, as described in the paper:
C   Gilbert, J.C. and Nocedal, J. (1992). "Global Convergence Properties 
C   of Conjugate Gradient Methods", SIAM Journal on Optimization, Vol. 2,
C   pp. 21-42. 
C   A web-based Server which solves unconstrained nonlinear optimization
C   problems using this Conjugate Gradient code can be found at:
C   http://www-neos.mcs.anl.gov/neos/solvers/UCO:CGPLUS/
C   --------------------------------------------------------------------
C
      SUBROUTINE CGFAM(N,X,F,G,D,GOLD,IPRINT,W,
     *                 IFLAG,IREST,METHOD,FINISH,ProMat,Rs)
C
C Subroutine parameters
C
      INTEGER N,oldN
      PARAMETER (oldN=66)
      DOUBLE PRECISION X(N),G(N),D(N),GOLD(N),W(N),F
     1  ,Rs(oldN),ProMat(oldN,oldN)
      INTEGER IPRINT(2),IFLAG,IREST,METHOD,IM,NDES
C
C     N      =  NUMBER OF VARIABLES
C     X      =  ITERATE
C     F      =  FUNCTION VALUE
C     G      =  GRADIENT VALUE
C     GOLD   =  PREVIOUS GRADIENT VALUE
C     IPRINT =  FREQUENCY AND TYPE OF PRINTING
C               IPRINT(1) < 0 : NO OUTPUT IS GENERATED
C               IPRINT(1) = 0 : OUTPUT ONLY AT FIRST AND LAST ITERATION
C               IPRINT(1) > 0 : OUTPUT EVERY IPRINT(1) ITERATIONS
C               IPRINT(2)     : SPECIFIES THE TYPE OF OUTPUT GENERATED;
C                               THE LARGER THE VALUE (BETWEEN 0 AND 3),
C                               THE MORE INFORMATION
C               IPRINT(2) = 0 : NO ADDITIONAL INFORMATION PRINTED
C           IPRINT(2) = 1 : INITIAL X AND GRADIENT VECTORS PRINTED
C           IPRINT(2) = 2 : X VECTOR PRINTED EVERY ITERATION
C           IPRINT(2) = 3 : X VECTOR AND GRADIENT VECTOR PRINTED 
C                       EVERY ITERATION 
C     EPS    =  CONVERGENCE CONSTANT !! not used ?? here -> deleted 
C     W      =  WORKING ARRAY
C     IFLAG  =  CONTROLS TERMINATION OF CODE, AND RETURN TO MAIN
C               PROGRAM TO EVALUATE FUNCTION AND GRADIENT
C               IFLAG = -3 : IMPROPER INPUT PARAMETERS
C               IFLAG = -2 : DESCENT WAS NOT OBTAINED
C               IFLAG = -1 : LINE SEARCH FAILURE
C               IFLAG =  0 : INITIAL ENTRY OR 
C                            SUCCESSFUL TERMINATION WITHOUT ERROR   
C               IFLAG =  1 : INDICATES A RE-ENTRY WITH NEW FUNCTION VALUES
C               IFLAG =  2 : INDICATES A RE-ENTRY WITH A NEW ITERATE
C     IREST  =  0 (NO RESTARTS); 1 (RESTART EVERY N STEPS)
C     METHOD =  1 : FLETCHER-REEVES 
C               2 : POLAK-RIBIERE
C               3 : POSITIVE POLAK-RIBIERE: BETA=MAX{BETA,0} 
C Local variables
C
      DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,ddot
      INTEGER MP,LP,ITER,NFUN,MAXFEV,INFO,I,NFEV,NRST,IDES
      LOGICAL NEW,FINISH
C
C     THE FOLLOWING PARAMETERS ARE PLACED IN COMMON BLOCKS SO THEY
C     CAN BE EASILY ACCESSED ANYWHERE IN THE CODE
C
C     MP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE REGULAR OUTPUT
C     LP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE ERROR OUPUT
      COMMON /CGDD/MP,LP
C
C     ITER: KEEPS TRACK OF THE NUMBER OF ITERATIONS
C     NFUN: KEEPS TRACK OF THE NUMBER OF FUNCTION/GRADIENT EVALUATIONS
      COMMON /RUNINF/ITER,NFUN
      SAVE
      DATA ONE/1.0D+0/,ZERO/0.0D+0/ 
c wq  store for script organisation:
      OPEN(24,FILE='cgRest.dat')
       IF(IFLAG.eq.1) then  
         rewind 24
         read(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
        endif 
C
C IFLAG = 1 INDICATES A RE-ENTRY WITH NEW FUNCTION VALUES
      IF(IFLAG.EQ.1) GOTO 72
C
C IFLAG = 2 INDICATES A RE-ENTRY WITH A NEW ITERATE
      IF(IFLAG.EQ.2) GOTO 80
C
C     INITIALIZE
C     ----------
C
C     IM =   NUMBER OF TIMES BETAPR WAS NEGATIVE FOR METHOD 2 OR
C            NUMBER OF TIMES BETAPR WAS 0 FOR METHOD 3
C
C     NDES = NUMBER OF LINE SEARCH ITERATIONS AFTER WOLFE CONDITIONS
C            WERE SATISFIED
C
      ITER= 0
      IF(N.LE.0) GO TO 96
      NFUN= 1
      NEW=.TRUE.
      NRST= 0
      IM=0
      NDES=0
C
      DO 5 I=1,N
 5    D(I)= -G(I)
cc
cc wq:  here in projection, use old N = N-1
cc Projector use to enforce the C(x)=0 condition, automatically 
c
       call PROJECTION(oldN,D,ProMat)
c
      GNORM= DSQRT(DDOT(N,D,1,D,1))
cc   write(44,*) 'N1-norm of D ', GNORM 
cc   here is used the shortened length of projected D !!
      STP1= ONE/GNORM
C
C     PARAMETERS FOR LINE SEARCH ROUTINE
C     ----------------------------------
C
C     FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION
C       OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE
C       DIRECTIONAL DERIVATIVE CONDITION ARE SATISFIED.
C
C     XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
C       WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C       IS AT MOST XTOL.
C
C     STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
C       SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP.
C
C     MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C       OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
C       MAXFEV BY THE END OF AN ITERATION.

      FTOL= 1.0D-4
      GTOL= 1.0D-1
      IF(GTOL.LE.1.D-04) THEN
        IF(LP.GT.0) WRITE(LP,145)
        GTOL=1.D-02
      ENDIF
      XTOL=   1.0D-17
c     STPMIN= 1.0D-20
c     STPMAX= 1.0D+20
c     MAXFEV= 20
cc wq new: 
        stpmin= 0.23081947D-5
        stpmax= 2.50000023081947d+0
        MAXFEV= 2*N
cc wq 
C
      IF(IPRINT(1).GE.0) CALL CGBD(IPRINT,ITER,NFUN,
     *   GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA)
C
C     MAIN ITERATION LOOP
C    ---------------------
C
 8    ITER= ITER+1

C     WHEN NRST>N AND IREST=1 THEN RESTART
      NRST= NRST+1
      INFO=0
C
C     CALL THE LINE SEARCH ROUTINE OF MOR'E AND THUENTE
C     (modified for our CG method)
C     -------------------------------------------------
C       JJ Mor'e and D Thuente, "Linesearch Algorithms with Guaranteed
C       Sufficient Decrease". ACM Transactions on Mathematical
C       Software 20 (1994), pp 286-307.
C
      NFEV=0
      DO 70 I=1,N
 70    GOLD(I)= G(I)
      DG= DDOT(N,D,1,G,1)
      DGOLD=DG
      STP=ONE
C
C Shanno-Phua's Formula For Trial Step
C
      IF(.NOT.NEW) STP= DG0/DG
      IF (ITER.EQ.1) STP=STP1
      IDES=0
      new=.false.
  72  CONTINUE
C wq  ?????  Notbremse !!!
      if(stp.lt.0.0d0) then 
        stp= 1.0D-3
        write(44,*) ' ERROR ? -stp is changed to 10^-3 '
      endif
      write(6,*) 'step  = ', stp
      write(44,*) 'step  = ', stp
c
C Call to the line search subroutine
C
      CALL CVSMOD(N,X,F,G,D,STP,FTOL,GTOL,
     *            XTOL,STPMIN,STPMAX,MAXFEV,INFO,NFEV,W,DG,DGOUT)

C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C         INFO = 0  IMPROPER INPUT PARAMETERS.
C         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
C         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
C                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
C         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C                   IS AT MOST XTOL.
C         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
C         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
C         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
C         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
C                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
C                   TOLERANCES MAY BE TOO SMALL.
c
      IF (INFO .EQ. -1) THEN
C       RETURN TO FETCH FUNCTION AND GRADIENT
        IFLAG=1
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
        RETURN
      ENDIF
c wq shut down some test which are to strong:       INFO5
      IF(INFO.EQ.4 .or. INFO.EQ.5 .or. INFO.EQ.6)  goto 91
      IF(INFO.NE.1) GO TO 90
91     continue
c wq
C
C     TEST IF DESCENT DIRECTION IS OBTAINED FOR METHODS 2 AND 3
C     ---------------------------------------------------------
C
      GG= DDOT(N,G,1,G,1)
      GG0= DDOT(N,G,1,GOLD,1)
      BETAPR= (GG-GG0)/GNORM**2
      IF (IREST.EQ.1.AND.NRST.GT.N) THEN
        NRST=0
        NEW=.TRUE.
        GO TO 75
      ENDIF 
C
      IF (METHOD.EQ.1) THEN
        GO TO 75
      ELSE
        DG1=-GG + BETAPR*DGOUT
        IF (DG1.lt. 0.0d0 ) GO TO 75
        IF (IPRINT(1).GE.0) write(6,*) 'no descent'
        IF (IPRINT(1).GE.0) write(44,*) 'ERROR, no descent'
        IDES= IDES + 1
        IF(IDES.GT.5) GO TO 95
        GO TO 72
      ENDIF
C
C     DETERMINE CORRECT BETA VALUE FOR METHOD CHOSEN
C     ----------------------------------------------
C
C     IM =   NUMBER OF TIMES BETAPR WAS NEGATIVE FOR METHOD 2 OR
C            NUMBER OF TIMES BETAPR WAS 0 FOR METHOD 3
C
C     NDES = NUMBER OF LINE SEARCH ITERATIONS AFTER WOLFE CONDITIONS
C            WERE SATISFIED
C
  75  NFUN= NFUN + NFEV
      NDES= NDES + IDES
      BETAFR= GG/GNORM**2
      IF (NRST.EQ.0) THEN
        BETA= ZERO
      ELSE
        IF (METHOD.EQ.1) BETA=BETAFR
        IF (METHOD.EQ.2) BETA=BETAPR
        IF ((METHOD.EQ.2.OR.METHOD.EQ.3).AND.BETAPR.LT.0) IM=IM+1
        IF (METHOD.EQ.3) BETA=MAX(ZERO,BETAPR)
      ENDIF
C
C     COMPUTE THE NEW DIRECTION
C     --------------------------
C
      DO 78 I=1,N
c      W(I) =  G(I)
  78  D(I) = -G(I) +BETA*D(I) 

c wq  here, in projection, use is old N = N-1
      call PROJECTION(oldN,D,ProMat)
c
      DG0= DGOLD*STP
C
C     RETURN TO DRIVER FOR TERMINATION TEST
C     -------------------------------------
C 
       GNORM=DSQRT(DDOT(N,G,1,G,1))
c
      IFLAG=2
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
      RETURN
  80  CONTINUE
C
C Call subroutine for printing output
C
      IF(IPRINT(1).GE.0) CALL CGBD(IPRINT,ITER,NFUN,
     *     GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA)
      IF (FINISH) THEN
         IFLAG = 0
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
         RETURN
      END IF
      GOTO 8
C
C     ----------------------------------------
C     END OF MAIN ITERATION LOOP. ERROR EXITS.
C     ----------------------------------------
C
  90  IFLAG=-1
      IF(LP.GT.0) WRITE(LP,100) INFO
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
      RETURN
  95  IFLAG=-2
      IF(LP.GT.0) WRITE(LP,135) I
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
      RETURN
  96  IFLAG= -3
      IF(LP.GT.0) WRITE(LP,140)
         rewind 24
         write(24,*) GTOL,GNORM,STP1,FTOL,XTOL,STPMIN,
     .       STPMAX,STP,BETA,BETAFR,BETAPR,DG0,GG,GG0,DGOLD,
     .       DGOUT,DG,DG1,d1,
     .       MP,LP,ITER,NFUN,MAXFEV,INFO,NFEV,NRST,IDES,
     .       NEW,FINISH ,IM,NDES
      RETURN
C
 100  FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED - ERROR RETURN',
     .          ' OF LINE SEARCH: INFO= ',I2)
cc   .       ,/,' POSSIBLE CAUSE: FUNCTION OR GRADIENT ARE INCORRECT')
 135  FORMAT(/' IFLAG= -2',/' DESCENT WAS NOT OBTAINED')
 140  FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N',
     .       ' IS NOT POSITIVE)')
 145  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',
     .       / ' IT HAS BEEN RESET TO 1.D-02')
      END
C
C     LAST LINE OF ROUTINE CGFAM
C
C**************************************************************************
      SUBROUTINE CGBD(IPRINT,ITER,NFUN,
     *           GNORM,N,X,F,G,STP,FINISH,NDES,IM,BETAFR,BETAPR,BETA)
C
C     ---------------------------------------------------------------------
C     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND AMOUNT
C     OF OUTPUT ARE CONTROLLED BY IPRINT.
C     ---------------------------------------------------------------------
C
      INTEGER N
      DOUBLE PRECISION X(N),G(N),F,GNORM,STP,BETAFR,BETAPR,BETA,QQ
      INTEGER IPRINT(2),ITER,NFUN,LP,MP,NDES,IM,Iqq,I
      LOGICAL FINISH
      COMMON /CGDD/MP,LP
C ??:
       QQ=betafr
       QQ=betapr
       Iqq=NDES
       Iqq= IM + Iqq* QQ
       QQ = Iqq  
c wq : ?? employ never used variables 
      IF (ITER.EQ.0)THEN
           PRINT*
           WRITE(MP,10)
           WRITE(MP,20) N
           WRITE(MP,30) F,GNORM
                 IF (IPRINT(2).GE.1)THEN
                     WRITE(MP,40)
                     WRITE(MP,50) (X(I),I=1,N)
                     WRITE(MP,60)
                     WRITE(MP,50) (G(I),I=1,N)
                 ENDIF
           WRITE(MP,10)
           WRITE(MP,70)
      ELSE
          IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
          IF (IPRINT(1).NE.0)THEN
               IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
                     IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70)
                     WRITE(MP,80)ITER,NFUN,F,GNORM,STP,BETA
               ELSE
                     RETURN
               ENDIF
          ELSE
               IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70)
               WRITE(MP,80)ITER,NFUN,F,GNORM,STP,BETA
          ENDIF
          IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
                  WRITE(MP,40)
                  WRITE(MP,50)(X(I),I=1,N)
              IF (IPRINT(2).EQ.3)THEN
                  WRITE(MP,60)
                  WRITE(MP,50)(G(I),I=1,N)
              ENDIF
          ENDIF
          IF (FINISH) WRITE(MP,100)
      ENDIF
C
 10   FORMAT('*************************************************')
 20   FORMAT(' N=',I5, / , ' INITIAL VALUES:')
 30   FORMAT(' F= ',1PD10.3,'   GNORM= ',1PD10.3)
 40   FORMAT(/,' VECTOR X= ')
 50   FORMAT(6(2X,1PD10.3/))
 60   FORMAT(' GRADIENT VECTOR G= ')
 70   FORMAT(/'   I  NFN',4X,'FUNC',7X,'GNORM',6X,
     *   'STEPLEN',4x,'BETA',/,
     *   ' ----------------------------------------------------')
 80   FORMAT(I4,1X,I3,2X,2(1PD10.3,2X),1PD8.1,2x,1PD8.1)
100   FORMAT(/' SUCCESSFUL CONVERGENCE (NO ERRORS).'
     *          ,/,' IFLAG = 0')
C
      RETURN
      END
C     LAST LINE OF CGBD
C*************************************************************************
C     DATA 
C     MP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE REGULAR OUTPUT
C     LP = UNIT NUMBER WHICH DETERMINES WHERE TO WRITE ERROR OUPUT
      BLOCK DATA CGCD
      COMMON /CGDD/MP,LP
      INTEGER LP,MP
      DATA MP,LP/44,44/
      END
C     LAST LINE OF CGFAM
c***********************************************************************
c BLAS  parts:
c***********************************************************************
c ddot
c daxpy
c dasum
c dcopy
c***********************************************************************
c matmult 
c***********************************************************************
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
c***********************************************************************
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
c***********************************************************************
      double precision function dasum(n,dx,incx)
c
c     takes the sum of the absolute values.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dtemp
      integer i,incx,m,mp1,n,nincx
c
      dasum = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dtemp = dtemp + dabs(dx(i))
   10 continue
      dasum = dtemp
      return
c
c        code for increment equal to 1
c        clean-up loop
c
   20 m = mod(n,6)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dabs(dx(i))
   30 continue
      if( n .lt. 6 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,6
        dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
     *  + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
   50 continue
   60 dasum = dtemp
      return
      end
c***********************************************************************
      subroutine  dcopy(n,dx,incx,dy,incy)
c
c     copies a vector, x, to a vector, y.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1)
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c        clean-up loop
c
   20 m = mod(n,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        dy(i) = dx(i)
        dy(i + 1) = dx(i + 1)
        dy(i + 2) = dx(i + 2)
        dy(i + 3) = dx(i + 3)
        dy(i + 4) = dx(i + 4)
        dy(i + 5) = dx(i + 5)
        dy(i + 6) = dx(i + 6)
   50 continue
      return
      end
c ***********************************************************
      subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:
c    m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
      real*8 m1,m2,mout,swap
      integer d1,d2,dim,i,j,k
      dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2)
      call mat_init(swap,d1,d2,0.0d0)
      do 10,i=1,d1
          do 20,j=1,d2
            do 30,k=1,dim
             swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j)
   30     continue
   20   continue
   10 continue
      call mat_copy(swap,mout,d1,d2)
      return
      end
cc
      subroutine mat_init(mat,nrow,ncol,value)
c   fills the NROW * NCOL matrix MAT with VALUE   
      real*8 mat,value
      integer nrow,ncol,i,j
      dimension mat(nrow,ncol)
      do 10,i=1,nrow
        do 20,j=1,ncol
         mat(i,j)=value
   20   continue
   10 continue  
      return
      end
cc
      subroutine mat_copy(matin,matout,nrow,ncol)
c copy MATIN(NROW,NCOL) to MATOUT(NROW,NCOL)
      real*8 matin,matout
      integer nrow,ncol,i,j
      dimension matin(nrow,ncol),matout(nrow,ncol)
      do 10,i=1,nrow
        do 20,j=1,ncol
         matout(i,j)=matin(i,j)
   20   continue
   10 continue  
      return
      end
