#include <iostream.h>
#include <math.h>

double pl[100]; // D510PL

const int ndcv = 10000; // D5120SI
int na, indflg[5];
double Z[ndcv], G[100], da[100], endflg, DF[100], X[10];

const int mxdat = 4500;  //D510UI
int ned[2];
double A[100], pl0[100], amx[100], amn[100], sigma[100], exda[mxdat];

const int ndcv0=10000; // D510UO
const int ierm=5000;

double R[100], error[ierm],Z0[ndcv0];

void fumili (double s, int m, int n1, int n2, double eps, 
	     double akappa, double alambd, int n3, bool it, int mc)
{
  
  double rp=1.e-14; // Relative precision on CDC6000 :))

  indflg[2]=0;
  if (it)
    cout<<"1 FUNCTION MINIMISATION BY SUBROUTINE FUMILI/LIKE,\n"
	<<"LM/0 IN THE FOLLOWING PRINT-OUT \n" 
	<< "0   S ="<<s<<" VALUE OF OBJECTIVE FUNCTION,\n"
	<<" EC = EXPECTED CHANGE IN S DURING NEXT ITERATION\n"
	<< "0 KAPPA ="<<akappa<<" ESTIMATED DISTANCE TO MINIMUM,\n"
	<< " LAMBDA ="<<alambd << " M="<<m
	<< "STEP LENGTH MODIFIER"<<endl;


  for(int i=0;i<100;i++)
    {
      amx[i]=1.e37;
      amn[i]=-1.e37;
    }
  endflg=0.;
  indflg[1]=0;
  for(int i=0;i<m;i++)
    {
      if(eps>0.) sigma[i]=0.;
      pl[i]=pl0[i];
    }

  int nn2=0,nn3=0,k,k1,i1,k2,l,l1;
  double t=0;
  double olds=0, gt=0;
  // Start new iteration
  int nn1 = 1;  // label 3....
  double t1=0.; 
  // Repeat iteration  with smaller step
  //label 4
 m4:
  s = 0.;
  int n0=0;
  for (int i=0;i<m;i++)
    {
      G[i]=0.;
      if(pl0[i]>0.) {
	n0++;
	if(pl[i]>0.) pl0[i]=pl[i];
      }
    }

  int nn0=n0*(n0+1)/2;
  for(int i=0;i<nn0;i++) Z[i]=0.;
  na = m;
  indflg[0]=0;
  // Calculate Objective function
  // call sgz(M,S)
  double sp = rp*fabs(s);
  for(int i=0;i<nn0;i++) Z0[i]=Z[i];
  //label 11
  //label 12
  if(nn3<0 && (nn1-n1)<0){
    t=2.*(s-olds-gt);
    if (indflg[0] == 0){
      if (fabs(s-olds)<=sp && -gt<=sp) goto m19;
      if (.59*t+gt>=0) t=-gt/t; else  goto m19;
      if (t-.25<0) t=.25;
    } else t=.25;
    gt=gt*t;
    t1=t1*t;
    nn2=0;
    for(int i=0;i<m;i++)
      if(pl[i]>0){
	A[i] = A[i]-da[i];
	pl[i] = pl[i]*t;
	da[i] = da[i]*t;
	A[i] = A[i]+da[i];
      }
    nn1++;
    goto m4;
  }
  m19:
  //label 19
  //  REMOVE CONTRIBUTION OF FIXED PARAMETERS FROM Z
  
  if (indflg[0]!=0)
    {
      endflg=-4.;
      goto m85;
    }
  i1=k1=k2=1;
  for (int i=0;i<m;i++)
    {
      if(pl0[i]>0){
	if (pl[i]==0) pl[i]=pl0[i];
	if (pl[i]<=0) k1 += i1;
	else {
	  if((A[i]>=amx[i] && G[i]<0) || 
	     (A[i]<=amn[i] && G[i]>0)) 
	    {
	      pl[i]=0.; 
	      k1 += i1;
	    } else
	    {
	      for (int j=0;j<i;j++)
		if(pl0[j]>0){
		  if(pl[j]>0) Z[k2++] = Z0[k1];
		  k1++;		  
		}

	    }
	}
	i1++;
      }
    }

  // Invert Z
  i1=l=1;
  
  for (int i=1;i<m;i++)
    if (pl[i]>0) {
      R[i]=Z[l];
      i1++;
      l += i1;
    }
  n0 = i1-1;
  // call mconv(n0)
  if (indflg[0] != 0) {
    indflg[0] = 0;
    indflg[1] = 1;
    goto m49;
  }
  // CALCULATE THEORETICAL STEP TO MINIMUM
  i1=1;
  for (int i=0;i<m;i++) {
    da[i]=0;
    if(pl[i]>0) {
      l1=1;
      for (l=0;l<m;l++)
	if(pl[l]>0){
	  if(i1-l1<=0) k = l1*(l1-1)/2+i1;
	  else         k = i1*(i1-1)/2+l1;
	  da[i]-=G[l]*Z[k];
	  l1++;
	}
      i1++;
    }
  }
  // CHECK FOR PARAMETERS ON BOUNDARY
  
 m49:
 m85:
  return;
}
  /*
      SUBROUTINE FUMILI (S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
      INTEGER I,J,L,L1,N0,NN0,NN1,NN2,NN3,N,IFIX1,K,K1,K2,I1,IFIX,IMAX
      REAL*8  FIXFLG,FI,T1,SP,RP,T,OLDS,GT,AFIX,SIGI,AKAP,BM,ABI,ABM
      REAL*8  BI,AL,AIMAX,AMB
CC
C... CHECK FOR PARAMETERS ON BOUNDARY
C
      AFIX = 0.D0
      IFIX = 0
      I1   = 1
      L    = I1
      DO I = 1,N
       IF (PL(I).GT.0.D0) THEN
        SIGI = DSQRT(DABS(Z(L)))
        R(I) = R(I)*Z(L)
        IF (EPS.GT.0.D0) SIGMA(I) = SIGI
        IF ((A(I).LT.AMX(I).OR.DA(I).LE.0.).AND.(A(I).GT.AMN(I).OR.
     1        DA(I).GE.0.)) GO TO 46
        AKAP = DABS(DA(I)/SIGI)
        IF (AKAP-AFIX.GT.0.D0) THEN
         AFIX  = AKAP
         IFIX  = I
         IFIX1 = I
        ENDIF
 46     I1 = I1 + 1
        L  = L  + I1
       ENDIF
      ENDDO
C
      IF (IFIX) 48,50,48
 48   PL(IFIX)=-1.D0
 49   FIXFLG=FIXFLG+1.D0
      FI = 0.D0
C
C... REPEAT CALCULATION OF THEORETICAL STEP AFTER FIXING EACH PARAMETER
C
      GO TO 19
C
C... CALCULATE STEP CORRECTION FACTOR
C
 50   ALAMBD = 1.D0
      AKAPPA = 0.D0
      IMAX   = 0
      DO I = 1,N
       IF (PL(I).GT.0.D0) THEN 
        BM  = AMX(I) - A(I)
        ABI = A(I)   + PL(I)
        ABM = AMX(I)
        IF (DA(I).LE.0.D0) THEN 
         BM  = A(I)-AMN(I)
         ABI = A(I)-PL(I)
         ABM = AMN(I)
        ENDIF
        BI = PL(I)
        IF (BI-BM.GT.0.D0) THEN  
         BI  = BM
         ABI = ABM
        ENDIF
        IF (DABS(DA(I))-BI.GT.0.D0) THEN 
         AL = DABS(BI/DA(I))
         IF (ALAMBD-AL.GT.0.D0) THEN 
           IMAX   = I
           AIMAX  = ABI
           ALAMBD = AL
         ENDIF
        ENDIF
        AKAP = DABS(DA(I)/SIGMA(I))
        IF (AKAP-AKAPPA.GT.0.D0) AKAPPA = AKAP
       ENDIF
      ENDDO
C
C... CALCULATE NEW CORRECTED STEP
C
      GT  = 0.D0
      AMB = 1.D18
      IF (ALAMBD.GT.0.D0) AMB = 0.25D0/ALAMBD
 
      DO I = 1,N
       IF (PL(I).GT.0.D0) THEN 
        IF (NN2-N2.GE.0.AND.DABS(DA(I)/PL(I))-AMB.GE.0.D0) THEN 
          PL(I) = 4.D0*PL(I)
          T1 = 4.
        ENDIF
        DA(I) = DA(I)*ALAMBD
        GT    = GT + DA(I)*G(I)
       ENDIF
      ENDDO
C
C... CHECK IF MINIMUM ATTAINED AND SET EXIT MODE
C
      IF (-GT.GT.SP.OR.T1.GE.1..OR.ALAMBD.GE.1.) GO TO 68
      ENDFLG=-1.D0
 68   IF (ENDFLG) 85,69,69
 69   IF (AKAPPA-DABS(EPS)) 70,75,75
 70   IF (FIXFLG) 72,71,72
 71   ENDFLG=1.D0
      GO TO 85
 72   IF (ENDFLG) 85,77,73
 73   IF (IFIX1) 85,85,76
 74   IF (FI-FIXFLG) 76,76,77
 75   IF (FIXFLG) 74,76,74
 76   FI=FI+1.D0
      ENDFLG=0.D0
 85   IF(ENDFLG.EQ.0.D0.AND.NN3.GE.N3) ENDFLG=-3.
      IF(ENDFLG.GT.0.D0.AND.INDFLG(2).GT.0) ENDFLG=-2.
      CALL MONITO (S,M,NN3,IT,EPS,GT,AKAPPA,ALAMBD)
      IF (ENDFLG) 83,79,83
C
C... CHECK IF FIXING ON BOUND IS CORRECT
C
 77   ENDFLG = 1.D0
      FIXFLG = 0.D0
      IFIX1  = 0
      DO I = 1,M
       PL(I) = PL0(I)
      ENDDO
      INDFLG(2) = 0
      GO TO 19
C
C... NEXT ITERATION
C
 79   ENDFLG = 0.D0
      DO I = 1,N
       A(I) = A(I) + DA(I)
      ENDDO
      IF (IMAX.GT.0) A(IMAX) = AIMAX
      OLDS = S
      NN2  = NN2 + 1
      NN3  = NN3 + 1
      GO TO 3
C
 83   MC = ENDFLG
      RETURN
C
C...  ENTRY FOR MAXIMUM LIKLEHOOD
C     ENTRY LIKELM
      ENTRY LIKELM(S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
      INDFLG(3)=1
      GO TO 1
C
 84   FORMAT('1',43X,'FUNCTION MINIMISATION BY SUBROUTINE FUMILI/LIKE',
     +'LM'/'0',55X,'IN THE FOLLOWING PRINT-OUT'/
     +     '0',27X,'S = VALUE OF OBJECTIVE FUNCTION,',
     + 'EC = EXPECTED CHANGE IN S DURING NEXT ITERATION'/
     +   '0',34X,'KAPPA = ESTIMATED DISTANCE TO MINIMUM,  LAMBDA =',
     + 'STEP LENGTH MODIFIER'///)
      END
C
  */