      PROGRAM ADM
*****************************************************************
*                                                               *
*     APPLICATION OF THE ALTERNATING DIRECTION METHOD OF        *
*     MULTIPLIERS TO SEPARABLE CONVEX PROGRAMMING PROBLEMS      *
*                                                               *
*                                  Coded by Masao Fukushima     *
*****************************************************************
      PARAMETER (NMAX=10, NVMAX=10, MMAX=30)
*
      REAL*8 X(NMAX,NVMAX), P(NMAX,MMAX), Z(NMAX,MMAX),
     -       Y(MMAX), R, EPS, Y0(MMAX), DYMAX, 
     -       FUN, CON(MMAX), FUN1,  
     -       INVHES(NVMAX,NVMAX), XX(NVMAX), OBJFUN, OBJGRA(NVMAX), 
     -       DIREC(NVMAX), NEWX(NVMAX), OLDX(NVMAX), 
     -       OLDGRA(NVMAX), NEWGRA(NVMAX),
     -       SS(NVMAX), YY(NVMAX), WW(NVMAX),CONS(MMAX),
     -       G(NMAX,NVMAX,NVMAX*(MMAX+1)),Q(NMAX,NVMAX*(MMAX+1)),
     -       B(NMAX,MMAX)
      INTEGER  AR(2*NVMAX,NVMAX)
*
      LOGICAL VALID,GRAD
      CHARACTER*12 OFNAME
      DATA EPS/1.0D-5/
      DATA KMAX/2500/
      DATA LEVEL/0/
      DATA GRAD/.TRUE./
*
*-----------------------------------------------------------------------
*-- 1 ----- Read in problem data.
*-----------------------------------------------------------------------
*
      WRITE(*,*) 'NUMBER OF BLOCKS: N = ?'
      READ(*,101) N
      WRITE(*,*) 'NUMBER OF VARIABLES IN EACH BLOCK: NV = ?'
      READ(*,101) NV
      WRITE(*,*) 'NUMBER OF CONSTRAINTS: M = ?'
      READ(*,101) M
      WRITE(*,*) 'SEED FOR RANDOM NUMBER GENERATOR: ISEED = ?'
      READ(*,101) ISEED
*
      WRITE(*,*) 'PENALTY PARAMETER: R = ?'
      READ(*,*) R
  101 FORMAT(I3)
*
*---------- Specify output file.
      WRITE(*,*) 'NAME OF OUTPUT FILE?'
      READ(*,100) OFNAME
  100 FORMAT(A)
*
      OPEN(6,FILE=OFNAME,STATUS='UNKNOWN')
*
      WRITE(6,105)
  105 FORMAT(
     - ' *** THE ALTERNATING DIRECTION METHOD OF MULTIPLIERS ***'/)
      WRITE(6,115) N,NV,M,ISEED
  115 FORMAT('    PROBLEM DATA (N,M,NV,ISEED) = (',I4,I4,I4,I6,' )')
*
*---------- Read in problem data.
      CALL INPUT1(N,NV,M,NMAX,MMAX,NVMAX,G,Q,B,AR,ISEED,VALID)
      IF(.NOT.VALID) THEN
	WRITE(6,120)
  120   FORMAT(' - INPUT DATA INVALID')
	STOP
      ENDIF
*
*---------- Print out input data
      WRITE(6,130) N
  130 FORMAT(/' NUMBER OF BLOCKS = ',I3)
      WRITE(6,132) NV
  132 FORMAT(/' NUMBER OF VARIABLES IN EACH BLOCK = ',I3)
      WRITE(6,134) M
  134 FORMAT(/' NUMBER OF CONSTRAINTS = ',I3)
C      WRITE(6,140)
C  140 FORMAT(/' INITIAL POINT:')
C      DO 1 J=1,N
C        WRITE(6,150) J,(X(J,JV),JV=1,NV)
C    1 CONTINUE
C  150 FORMAT(' ',5X,'X(',I2,') = ',10F12.5)
*
C      WRITE(6,160) R
C  160 FORMAT(/'     PENALTY PARAMETER R = ',F12.5)
*
*---------- Initialize dual variables (Step 1)
      DO 2 J=1,N
      DO 2 I=1,M
	P(J,I)=0.0
	Z(J,I)=0.0
    2 CONTINUE
*
      K=0
      ITER=0
      K1=0
      K2=0
      K3=0
      K4=0
*
*---------- Compute vector Y(K+1) (Step 2)
*
 1000 CONTINUE
*
C      WRITE(6,200) K
C  200 FORMAT(/'ITERATION * ',I5)
C      DO 3 J=1,N
C        WRITE(6,210) J,(P(J,I),I=1,M)
C    3 CONTINUE
C  210 FORMAT(' ',5X,'P(',I2,') = ',10F12.5)
C      DO 4 J=1,N
C        WRITE(6,220) J,(Z(J,I),I=1,M)
C    4 CONTINUE
C  220 FORMAT(' ',5X,'Z(',I2,') = ',10F12.5)
      IF(K.NE.0) THEN
	DO 30 I=1,M
	  Y0(I)=Y(I)
   30   CONTINUE
      ENDIF
      DO 6 I=1,M
	Y(I)=0.0
	DO 5 J=1,N
	  Y(I)=Y(I)+Z(J,I)-P(J,I)/R
    5   CONTINUE
	Y(I)=Y(I)/N
    6 CONTINUE
C      WRITE(6,230) (Y(I),I=1,M)
C  230 FORMAT(' ',5X,'Y     = ',10F12.5)
*
*---------- Convergence test
      IF(K.NE.0) THEN
	DYMAX=0.0
	DO 31 I=1,M
	  IF(DYMAX.LT.DABS(Y(I)-Y0(I))) DYMAX=DABS(Y(I)-Y0(I))
   31   CONTINUE
	IF(DYMAX.LT.EPS) THEN
	  GOTO 500
	ELSE
C
	  IF(DYMAX.LT.1.0D-1 .AND. K1.EQ.0) K1=K
	  IF(DYMAX.LT.1.0D-2 .AND. K2.EQ.0) K2=K
	  IF(DYMAX.LT.1.0D-3 .AND. K3.EQ.0) K3=K
	  IF(DYMAX.LT.1.0D-4 .AND. K4.EQ.0) K4=K
C
	  GOTO 600
	ENDIF
      ELSE
	GOTO 600
      ENDIF
*
  500 CONTINUE
C      WRITE(6,290) 
  290 FORMAT(
     -' *****************************************************')
      WRITE(6,300) K
  300 FORMAT(/' * OPTIMAL SOLUTION WAS OBTAINED AT ITERATION'
     -     ,I5/
C     - ' +++ ALTERNATING DIRECTION METHOD OF MULTIPLIERS +++'/
     -       )
C      WRITE(6,290) 
      WRITE(6,301) EPS
  301 FORMAT('    TOLERANCE: EPS = ',D7.1/)
*
      GOTO 450
*
  400 WRITE(6,290)
      WRITE(6,410) KMAX
  410 FORMAT(/' NUMBER OF ITERATIONS EXCEEDS  '
     -    ,I5/
     -  ' +++ ALTERNATING DIRECTION METHOD OF MULTIPLIERS +++'/)
      WRITE(6,290)
*
  450 CONTINUE
*
      DO 32 J=1,N
	WRITE(6,310) J,(X(J,JV),JV=1,NV)
  310   FORMAT('   X(',I2,') = ',5F12.5/(11X,5F12.5))
   32 CONTINUE
*
C      DO 33 J=1,N
C        WRITE(6,320) J,(P(J,I),I=1,M)
C  320   FORMAT('   P(',I2,') = ',5F12.5/11X,5F12.5)
C   33 CONTINUE
*
C      DO 34 J=1,N
C        WRITE(6,330) J,(Z(J,I),I=1,M)
C  330   FORMAT('   Z(',I2,') = ',5F12.5/11X,5F12.5)
C   34 CONTINUE
*
      WRITE(6,340) (Y(I),I=1,M)
  340 FORMAT('   Y     = ',5F12.5/(11X,5F12.5))
*
      FUN1=0.0
      DO 37 I=1,M
	CONS(I)=0.0
   37 CONTINUE
      DO 36 J=1,N
	DO 35 JV=1,NV
	  XX(JV)=X(J,JV)
   35   CONTINUE
	CALL EVAL(J,NV,NVMAX,M,NMAX,MMAX,YI,P,R,XX,
     -          OBJFUN,OBJGRA,FUN,CON,G,Q,B,GRAD)
C        WRITE(6,350) J,(CON(I),I=1,M)
C  350   FORMAT('   C(',I2,') = ',5F12.5/11X,5F12.5)
	FUN1=FUN1+FUN
	DO 38 I=1,M
	  CONS(I)=CONS(I)+CON(I)
   38   CONTINUE
   36 CONTINUE
      WRITE(6,360) FUN1
  360 FORMAT(/'  OBJECTIVE = ',F12.5)
      WRITE(6,375) (CONS(I),I=1,M)
  375 FORMAT(/'  CONSTRAINTS = ',5F11.5/(16X,5F11.5))
      WRITE(6,370) R
  370 FORMAT(/'  PENALTY PARAMETER: R = ',F10.5)
      WRITE(6,380) ITER
  380 FORMAT(/'  CUMULATIVE NO. OF MINOR ITERATIONS = ',I4)
      BITER=FLOAT(ITER)/N
      WRITE(6,385) BITER
  385 FORMAT(/'  NO. OF MINOR ITERATIONS / BLOCK = ',F12.5)
      WRITE(6,390) K1,K2,K3,K4
  390 FORMAT(/'  K1 = ',I3,'  K2 = ',I3,
     -        '  K3 = ',I3,'  K4 = ',I3)
      STOP
*
*---------- Solve subproblems (Step 3)
  600 DO 50 J=1,N
	DO 10 JV=1,NV
	  XX(JV)=X(J,JV)
   10   CONTINUE
*
	 CALL BFGS(J,NV,NVMAX,M,NMAX,MMAX,Y,P,R,XX,
     -             OBJFUN,OBJGRA,FUN,CON,INVHES,DIREC,NEWX,
     -             OLDX,OLDGRA,NEWGRA,SS,YY,WW,GRAD,G,Q,B,IT,INDEX)
*
	IF(INDEX.NE.0) THEN
	  WRITE(6,240) INDEX, J
  240     FORMAT(/'   INDEX = ',I3,
     -           ' * ABNORMAL TERMINATION OCCURS AT SUBPROBLEM',I3)
	  STOP
	ELSE
	  ITER=ITER+IT
C          WRITE(6,250) J, IT
C  250     FORMAT(/'   SOLUTION OF SUBPROBLEM',I3,3X,
C     -             '(NO. OF ITERATIONS =',I3,')')
C          WRITE(6,260) (XX(I),I=1,NV)
C  260     FORMAT('   * VARIABLES X = ',10F12.5)
C          WRITE(6,270) (OBJGRA(I),I=1,NV)
C  270     FORMAT('   *    OBJGRA   = ',10D12.5)
	ENDIF
*
	DO 11 JV=1,NV
	  X(J,JV)=XX(JV)
   11   CONTINUE
*
	DO 12 I=1,M
	  Z(J,I)=Y(I)+(P(J,I)+CON(I))/R
	  IF(Z(J,I).LT.0.0) Z(J,I)=0.0
   12   CONTINUE
   50 CONTINUE
*
*---------- Update P(J) (Step 4)
      DO 13 J=1,N
	DO 13 I=1,M
	  P(J,I)=P(J,I)+R*(Y(I)-Z(J,I))
   13 CONTINUE
*
      K=K+1
      IF(K.GT.KMAX) GOTO 400
      GOTO 1000
*
      END
*
      SUBROUTINE INPUT1(N,NV,M,NMAX,MMAX,NVMAX,G,Q,B,AR,ISEED,VALID)
************************************************************************
*                                                                      *
*     Read in initial data from file.                                  *
*     ------------------------                                         *
*     N: number of blocks                                              *
*     NV: number of variables in each block.                           *
*     M: number of constraints.                                        *
*     X: array containing initial values of the variables.             *
*     VALID: logical variable; if the input value of N is invalid,     *
*       then set VALID=.FALSE. and return to the main program.         *
*                                                                      *
************************************************************************
*
      REAL*8 G(NMAX,NVMAX,NVMAX*(MMAX+1)),Q(NMAX,NVMAX*(MMAX+1)),
     -       B(NMAX,MMAX)
      INTEGER AR(2*NVMAX,NVMAX)
      INTEGER*4 ISEED
      LOGICAL VALID
*
      VALID=.TRUE.
*
*---------- Number of blocks.
      IF(N.LE.0 .OR. N.GT.NMAX) THEN
	WRITE(6,100)
  100   FORMAT(' INVALID N')
      ENDIF
*
*---------- Number of variables in each block.
      IF(NV.LE.0 .OR. NV.GT.NVMAX) THEN
	WRITE(6,200)
  200   FORMAT(' INVALID NV')
      ENDIF
*
*---------- Number of constraints.
      IF(M.LE.0 .OR. M.GT.MMAX) THEN
	WRITE(6,210)
  210   FORMAT(' INVALID M')
      ENDIF
*
*
*---------- Coefficients of the problem functions
*
*
      DO 5 JB=1,N
	DO 5 IC=0,M
	  DO 4 L=1,2*NV
	  DO 4 J=1,NV
	    IDATA=INT(10.0*RAND(ISEED)-5.0)
	    AR(L,J)=IDATA
    4     CONTINUE
	  DO 5 I=1,NV
	  DO 5 J=1,NV
	    G(JB,I,NV*IC+J)=0.0
	    DO 5 L=1,2*NV
	      G(JB,I,NV*IC+J)=G(JB,I,NV*IC+J)
     -         +AR(L,I)*AR(L,J)
    5 CONTINUE
*
      DO 10 JB=1,N
	DO 7 J=1,NV
	  IDATA=INT(200.0*RAND(ISEED)-100.0)
	  Q(JB,J)=IDATA
    7   CONTINUE
	DO 8 J=1,NV*M
	  IDATA=INT(2.0*RAND(ISEED)-1.0)
	  Q(JB,NV+J)=IDATA
    8   CONTINUE
	DO 9 I=1,M
	  IDATA=INT(10.0*RAND(ISEED)+10.0)
	  B(JB,I)=IDATA
    9   CONTINUE
   10 CONTINUE
*
      RETURN
      END
*
*
      REAL FUNCTION RAND(IR)
******************************************************************
*     Random number generator                                    *
******************************************************************
*
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (ML=1664501)
      PARAMETER (IA=1229, IC=351750)
      PARAMETER (ANORM=1.0/ML)
*
      IR=MOD(IA*IR+IC,ML)
      RAND=IR*ANORM
*
      RETURN
      END
*
*
      SUBROUTINE BFGS(JB,NV,NVMAX,M,NMAX,MMAX,YI,P,R,X,
     -             OBJFUN,OBJGRA,FUN,CON,INVHES,DIREC,
     -             NEWX,OLDX,OLDGRA,NEWGRA,S,Y,W,GRAD,G,Q,B,IT,
     -             INDEX)
******************************************************************
*     Solve jth subproblem by BFGS method                        *
******************************************************************
      REAL*8 P(NMAX,MMAX),X(NVMAX),YI(MMAX),R,  
     -       OBJFUN,OBJGRA(NVMAX),FUN,CON(MMAX), 
     -       INVHES(NVMAX,NVMAX), 
     -       DIREC(NVMAX),NEWX(NVMAX),OLDX(NVMAX), 
     -       OLDGRA(NVMAX),S(NVMAX),Y(NVMAX),W(NVMAX),
     -       NEWFUN,OLDFUN,NEWGRA(NVMAX),
     -       GNORM,SNORM,YNORM,GNORM1,DNORM1, 
     -       SY,YHY,C1,C2,EPS1,EPS2,EPS3,BOUND,ETA,
     -       G(NMAX,NVMAX,NVMAX*(MMAX+1)),Q(NMAX,NVMAX*(MMAX+1)),
     -       B(NMAX,MMAX)
      LOGICAL GRAD
* 
      DATA EPS1,EPS2,EPS3/1.0D-4,1.0D-6,1.0D-8/
      DATA BOUND,ETA/1.0D+20,1.0D-4/
*
      N=NV
*
*---------- Set the initial inverse Hessian approximation
*           to be the identity matrix.
      DO 1 I=1,N
	INVHES(I,I)=1.0D0
	DO 1 J=1,I-1
	  INVHES(I,J)=0.0
	  INVHES(J,I)=0.0
    1 CONTINUE
*
      INDEX=0
      IT=0
      ITMAX=MAX0(100,20*N)
*
*---------- Iteration.
*
*---------- Compute search direction.
    2 CONTINUE
      CALL EVAL(JB,NV,NVMAX,M,NMAX,MMAX,YI,P,R,X,
     -          OBJFUN,OBJGRA,FUN,CON,G,Q,B,GRAD)
* 
C
      GNORM1=0.0
      DO 50 I=1,N
	IF(GNORM1.LT.DABS(OBJGRA(I))) GNORM1=DABS(OBJGRA(I))
   50 CONTINUE
      IF(GNORM1.LT.EPS1) RETURN
C
      DO 3 I=1,N
	DIREC(I)=0.0
	DO 3 J=1,N
	  DIREC(I)=DIREC(I)-INVHES(I,J)*OBJGRA(J)
    3 CONTINUE
C
      DNORM1=0.0
      DO 51 I=1,N
	IF(DNORM1.LT.DABS(DIREC(I))) DNORM1=DABS(DIREC(I))
   51 CONTINUE
      IF(DNORM1.LT.EPS1) RETURN
C
*
*---------- Perform line search.
      CALL LSEARC(JB,NV,NVMAX,M,NMAX,MMAX,YI,P,R,X,
     -     NEWX,OBJFUN,NEWFUN,OBJGRA,NEWGRA,FUN,CON,G,Q,B,GRAD,
     -     DIREC,EPS3,INDEX)
      IF(INDEX.EQ.5) RETURN
*
*---------- Compute vectors S and Y.
      OLDFUN=OBJFUN
      OBJFUN=NEWFUN
      DO 4 I=1,N
	OLDX(I)=X(I)
	X(I)=NEWX(I)
	S(I)=X(I)-OLDX(I)
	OLDGRA(I)=OBJGRA(I)
	OBJGRA(I)=NEWGRA(I)
    4 CONTINUE
*
      DO 5 I=1,N
	Y(I)=OBJGRA(I)-OLDGRA(I)
    5 CONTINUE
*
*---------- Check termination criteria.
      DO 6 I=1,N
	IF(DABS(X(I)).GT.BOUND) THEN
	  INDEX=4
	  RETURN
	ENDIF
    6 CONTINUE
*
      GNORM=0.0
      DO 7 I=1,N
	GNORM=GNORM+OBJGRA(I)*OBJGRA(I)
    7 CONTINUE
      GNORM=DSQRT(GNORM)
*
      IF (GNORM.LT.EPS1) RETURN
*
      IF(IT.GT.ITMAX) THEN
	IF(GNORM.LT.EPS1) THEN
	  INDEX=1
	ELSEIF(OLDFUN-OBJFUN.LT.EPS2*(1.0+DABS(OLDFUN))) THEN
	  INDEX=2
	ELSE
	  INDEX=3
	ENDIF
	RETURN
      ENDIF
*
*---------- Update the inverse Hessian approximation.
      SNORM=0.0
      YNORM=0.0
      SY=0.0
      DO 8 I=1,N
	SNORM=SNORM+S(I)*S(I)
	YNORM=YNORM+Y(I)*Y(I)
	SY=SY+S(I)*Y(I)
    8 CONTINUE
      SNORM=DSQRT(SNORM)
      YNORM=DSQRT(YNORM)
*
      IF(SY.LT.ETA*SNORM*YNORM) THEN
	IF(LEVEL.EQ.2) THEN
	  WRITE(6,100)
  100     FORMAT(/' HESSIAN UPDATE IS SKIPPED AT ITERATION ',I5,
     -           /' TO PRESERVE POSITIVE DEFINITENESS')
	ENDIF
	GOTO 12
      ENDIF
*
*
      YHY=0.0
      DO 10 I=1,N
	W(I)=0.0
	DO 9 J=1,N
	  W(I)=W(I)+INVHES(I,J)*Y(J)
    9   CONTINUE
	YHY=YHY+Y(I)*W(I)
   10 CONTINUE
*
      C2=1.0D0/SY
      C1=(1.0D0+C2*YHY)*C2
*
      DO 11 I=1,N
	DO 11 J=1,I
	  INVHES(I,J)=INVHES(I,J)+C1*S(I)*S(J)-C2*(S(I)*W(J)+W(I)*S(J))
	  INVHES(J,I)=INVHES(I,J)
   11 CONTINUE
*
   12 IT=IT+1
      GOTO 2
*
      END
*
*
*
      SUBROUTINE LSEARC(JB,N,NVMAX,M,NMAX,MMAX,YI,P,R,X,
     -     NEWX,OBJFUN,NEWFUN,OBJGRA,NEWGRA,FUN,CON,G,Q,B,GRAD,
     -     DIREC,EPS,INDEX)
************************************************************************
*                                                                      *
*     Perform line search by safeguarded quadratic interpolation.      *
*     ----------------------------                                     *
*     X: array containing vector of current solution.                  *
*     DIREC: array containing vector of search direction.              *
*     DIRDER: directional derivative of the objective function at      *
*             point X with respect to direction DIREC.                 *
*     ST: step-size determined by line search.                         *
*     NEWX: array containing the vector X+ST*DIREC.                    *
*     EPS: minimum step-size.                                          *
*     BETA: parameter used in the line search criterion (default       * 
*       = 1.0D-3).                                                     *
*     LAMBDA, MU: safeguard parameters (default values LAMBDA=1.0D-1,  *
*       MU=5.0D-1).                                                    *
*                                                                      *
************************************************************************
*
      REAL*8 X(N),NEWX(N),OBJFUN,NEWFUN,OBJGRA(N),DIREC(N),
     -       P(NMAX,MMAX),YI(MMAX),R,FUN,CON(MMAX),  
     -       G(NMAX,NVMAX,NVMAX*(MMAX+1)),Q(NMAX,NVMAX*(MMAX+1)),
     -       B(NMAX,MMAX),NEWGRA(N),
     -       DIRDER,DNORM,C,BETA,LAMBDA,MU,EPS,MINST,ST1,ST2,ST
      LOGICAL GRAD
      DATA BETA,LAMBDA,MU/1.0D-3,1.0D-1,5.0D-1/
*
*---------- Compute the directional derivative.
      DIRDER=0.0
      DNORM=0.0
      DO 1 I=1,N
	DIRDER=DIRDER+OBJGRA(I)*DIREC(I)
	DNORM=DNORM+DIREC(I)*DIREC(I)
    1 CONTINUE
*
      DNORM=DSQRT(DNORM)
C
C      WRITE(6,100) DIRDER, DNORM
C  100 FORMAT(/'     DIRDER    = ',D12.5,'    DNORM = ',D12.5)
C      WRITE(6,170) OBJFUN
C  170 FORMAT('     OBJFUN = ',F12.5)
C
      MINST=DMIN1(1.0D0,EPS/DNORM)
*
*---------- Set the initial step-size to be unity.
      ST=1.0D0
*
*---------- Check if the current step-size is too small.
    2 IF(ST.LT.MINST) THEN
	INDEX=5
	RETURN
      ENDIF
*
*---------- Test the line search criterion.
      DO 3 I=1,N
	NEWX(I)=X(I)+ST*DIREC(I)
    3 CONTINUE
* 
      CALL EVAL(JB,N,NVMAX,M,NMAX,MMAX,YI,P,R,NEWX,
     -          NEWFUN,NEWGRA,FUN,CON,G,Q,B,GRAD)
*
      IF(NEWFUN.LT.OBJFUN+BETA*DIRDER*ST) THEN
C
C        WRITE(6,110) ST
C  110   FORMAT('     STEP SIZE = ',D12.5)
C        WRITE(6,170) NEWFUN
C
*
	RETURN
      ENDIF
*
*---------- Calculate the new step-size.
      C=(NEWFUN-OBJFUN-DIRDER*ST)/ST**2
      ST1=-0.5*DIRDER/C
      ST2=DMIN1(ST1,MU*ST)
      ST=DMAX1(LAMBDA*ST,ST2)
*
      GOTO 2
*
      END

      SUBROUTINE EVAL(JB,NV,NVMAX,M,NMAX,MMAX,YI,P,R,XX,
     -          OBJFUN,OBJGRA,FUN,CON,G,Q,B,GRAD)
************************************************************************
*     Evaluate function/gradient values of problem functions           *
************************************************************************
*
      REAL*8 P(NMAX,MMAX),XX(NVMAX),YI(MMAX),R,  
     -       OBJFUN,OBJGRA(NVMAX),FUN,CON(MMAX), 
     -       G(NMAX,NVMAX,NVMAX*(MMAX+1)),Q(NMAX,NVMAX*(MMAX+1)),
     -       B(NMAX,MMAX),CCC
      LOGICAL GRAD
*
      FUN=0.0
      DO 1 J=1,NV
      DO 1 I=1,NV
	FUN=FUN+0.5*G(JB,I,J)*XX(I)*XX(J)
    1 CONTINUE
      DO 2 J=1,NV
	FUN=FUN+Q(JB,J)*XX(J)
    2 CONTINUE
*
      DO 100 IC=1,M
	CON(IC)=-B(JB,IC)
	DO 3 J=1,NV
	DO 3 I=1,NV
	  CON(IC)=CON(IC)+0.5*G(JB,I,(NV*IC)+J)*XX(I)*XX(J)
    3   CONTINUE
	DO 4 J=1,NV
	  CON(IC)=CON(IC)+Q(JB,(NV*IC)+J)*XX(J)
    4   CONTINUE
  100 CONTINUE
*
      OBJFUN=FUN
      DO 5 J=1,NV
	OBJGRA(J)=Q(JB,J)
	DO 5 I=1,NV
	  OBJGRA(J)=OBJGRA(J)+G(JB,I,J)*XX(I)
    5 CONTINUE
*
      DO 200 IC=1,M
	CCC=YI(IC)+(P(JB,IC)+CON(IC))/R
	IF(CCC.GT.0.0) THEN
	  OBJFUN=OBJFUN+0.5*R*CCC*CCC
	  DO 150 J=1,NV
	    OBJGRA(J)=OBJGRA(J)+CCC*Q(JB,(NV*IC)+J)
	    DO 150 I=1,NV
	      OBJGRA(J)=OBJGRA(J)+CCC*G(JB,I,(NV*IC)+J)*XX(I)
  150     CONTINUE
	ENDIF
  200 CONTINUE
*
      RETURN
      END
