In this chapter we discuss:
Minimal descriptions of these user exit subroutines are as follows.
EKKBRNU, EKKCHNU, EKKEVNU and EKKSLVU can be used initially for diagnostic purposes to see what choices are being made. After noting the choices, you may decide to alter them with these subroutines. Note that EKKEVNU is called from the version of EKKSLVU supplied with the library, and so if a user-supplied version of EKKSLVU is called, then EKKEVNU may not be called.
EKKCUTU, EKKHEUU, and EKKNODU, do not provide information on the branch and bound strategy, but instead provide you with a way of altering the integer solution algorithm in complex ways. They should be used only by experienced users to solve problems whose solution by other means is impractical or even impossible.
Default user exit subroutines are provided. So if you don't want to gain control at an exit point, do nothing, and execution will continue normally. If, however, you do code your own subroutines with the names specified below, and link them before the Optimization Library when running your job, your user exit subroutine will gain control, rather than the stub subroutine provided.
When your subroutine has control, it could record statistics, report progress, intercept messages, alter or stop execution, allow execution to continue without change, and so on.
User exit subroutines may be coded in C or FORTRAN, but since the current library code is not re-entrant, these subroutines should not make calls that could result in recursion. Sample user exit routines, in C and FORTRAN, are included with the product.
An example main program that uses the "Save Node Information User Exit Subroutine (EKKNODU),", pseudo-code for the call to EKKBRNU, and also a sample EKKBRNU user exit subroutine are given below.
C*********************************************************************** C C EXMSLV2 C C This program solves a mixed-integer problem and uses the user C exit routines EXCHNU, EXBRNU, and EXEVNU to modify the node C selection, branching technique, and node evaluation technique C that would be selected by default by EKKMSLV. C C*********************************************************************** C PROGRAM MAIN C C Bring in include files with control variable definitions. IMPLICIT NONE INCLUDE (OSLI) INCLUDE (OSLR) C C Allocate dspace. INTEGER*4 MAXSPC,RTCOD PARAMETER (MAXSPC=1000000) REAL*8 DSPACE(MAXSPC) COMMON/BIG/DSPACE C C Describe application and specify that there is 1 model. CALL EKKDSCA(RTCOD,DSPACE,MAXSPC,1) IF (RTCOD.GT.0) CALL CHKRT('EKKDSCA',RTCOD) C C Describe the model as having 1 block. CALL EKKDSCM(RTCOD,DSPACE,1,1) IF (RTCOD.GT.0) CALL CHKRT('EKKDSCM',RTCOD) C C Read model data from MPS file at unit 98. C Unit 15 will used to store branch and bound information. CALL EKKMPS(RTCOD,DSPACE,98,2,15) IF (RTCOD.GT.0) CALL CHKRT('EKKMPS ',RTCOD) C C Scale the coefficient matrix. CALL EKKSCAL(RTCOD,DSPACE) IF (RTCOD.GT.0) CALL CHKRT('EKKSCAL',RTCOD) C C Create a vector copy of the matrix. CALL EKKNWMT(RTCOD,DSPACE,3) IF (RTCOD.GT.0) CALL CHKRT('EKKNWMT',RTCOD) C C Solve problem using simplex method. CALL EKKSSLV(RTCOD,DSPACE,1,2) IF (RTCOD.GT.0) CALL CHKRT('EKKSSLV',RTCOD) C C Modify a number of integer control variables. CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKIGET',RTCOD) C Switch off printing of simplex log. ILOGLEVEL=0 C Stop at first solution. IMAXSOLS=1 CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKISET',RTCOD) C C Choose a pessimistic target so that code will go as fast as C possible to any solution. Since RDEGSCALE is zero all estimates C will be scaled so that the continuous estimate is RTARGET. Since C RTARGET is pessimistic, each branch should be cheaper than expected, C and this can speed up the search for a feasible solution. C C Modify a number of real control variables. CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) IF (RTCOD.GT.0) CALL CHKRT('EKKRGET',RTCOD) RDEGSCALE = 0.0D0 RTARGET = 1000.0D0 CALL EKKRSET(RTCOD,DSPACE,OSLR,OSLRLN) IF (RTCOD.GT.0) CALL CHKRT('EKKRSET',RTCOD) C C Solve mixed-integer programming problem. CALL EKKMSLV(RTCOD,DSPACE,1,35,50) IF (RTCOD.GT.0) CALL CHKRT('EKKMSLV',RTCOD) C C Now that a solution has been found, the scaling can be reduced C so that the objective function is given more weight. CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) IF (RTCOD.GT.0) CALL CHKRT('EKKRGET',RTCOD) RDEGSCALE=1.0D-2*RDEGSCALE CALL EKKRSET(RTCOD,DSPACE,OSLR,OSLRLN) IF (RTCOD.GT.0) CALL CHKRT('EKKRSET',RTCOD) C C Do not stop until the end. CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKIGET',RTCOD) IMAXSOLS=999999 CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKISET',RTCOD) C C Solve mixed-integer programming problem. CALL EKKMSLV(RTCOD,DSPACE,2,35,50) IF (RTCOD.GT.0) CALL CHKRT('EKKMSLV',RTCOD) C C Modify the print mask. CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKIGET',RTCOD) IPRTINFOMASK=511 CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKISET',RTCOD) C C Print the solution. CALL EKKPRTS(RTCOD,DSPACE) IF (RTCOD.GT.0) CALL CHKRT('EKKPRTS',RTCOD) C STOP END C C*********************************************************************** C This subroutine prints the character string RTNAME and the return C code RTCOD and stops if RTCOD is large enough to indicate that an C OSL error or severe error has occurred. C*********************************************************************** C SUBROUTINE CHKRT(RTNAME,RTCOD) CHARACTER*7 RTNAME INTEGER*4 RTCOD C WRITE(6,9000) RTNAME,RTCOD IF (RTCOD.GE.200) STOP 16 RETURN 9000 FORMAT (1X,'********** ',A7,' return code of ',I4,' **********') END
C EKKBRNU is called while EKKMSLV is choosing a branch: C No sets have yet been processed. C (Sets refers to special ordered sets and C integer variable sets.) . CALL EKKBRNU(DSPACE,MSPACE, 1 ,IARRAY,RARRAY,USERRTCD) . C Loop over all the sets. DO 200 I = 1, "the number of sets" . . CALL EKKBRNU(DSPACE,MSPACE, 2 ,IARRAY,RARRAY,USERRTCD) C Here you could set USERRTCD to instigate a call to EKKBRNU C for each variable in set I (i.e., REASON=3 below). .
C Loop over each variable in a set. DO 100 J = 1, "the number of variables in set I" . . . CALL EKKBRNU(DSPACE,MSPACE, 3 ,IARRAY,RARRAY,USERRTCD) C Here you might do your own branch consideration. 100 CONTINUE . CALL EKKBRNU(DSPACE,MSPACE, 4 ,IARRAY,RARRAY,USERRTCD) C Here you may change the choice of variable to branch on. . 200 CONTINUE . CALL EKKBRNU(DSPACE,MSPACE, 1 ,IARRAY,RARRAY,USERRTCD) C Here you may change the choice of set to branch on. .
C*********************************************************************** C C EXBRNU C C This user exit routine modifies the default choice of the C branching variable. C C IREASN is the situation when EKKBRNU is called, where: C C IREASN=1: subroutine is being called before any set is processed; C MARRAY and DARRAY are not used. C =2: subroutine is being called before processing each set; C only set number, type of set, priority, and number in set C are valid in MARRAY and DARRAY. C =3: subroutine is being called after processing each C variable in the set; all of MARRAY and DARRAY are valid C and you may change pseudocosts (takes effect on the next C branch). C =4: subroutine is being called after processing each set; C all of MARRAY and DARRAY are valid and the user may C change variable number and priority (takes effect on C next branch). C =5: subroutine is being called after processing all sets; C all of MARRAY and DARRAY are valid and the set number in C OSL's chosen set, which you may now change. C C MARRAY(1): set number being processed. C MARRAY(2): type of set being processed. C MARRAY(3): priority (1 is highest). C MARRAY(4): number of variables in the set C MARRAY(5): variable column number of the current variable. C MARRAY(6): direction of the branching. C DARRAY(1): current value of variable in continuous solution. C DARRAY(2): down pseudocost. C DARRAY(3): up pseudocost for single variables. C For sets it is the reference row entry. C DARRAY(4): estimated degradation for the down branch. C DARRAY(5): estimated degradation for the up branch. C C NOTE THIS SAMPLE USER EXIT IS NOT VALID FOR MIP PROBLEMS WITH SOS C C*********************************************************************** C SUBROUTINE EKKBRNU(DSPACE,MSPACE,IREASN,MARRAY,DARRAY,JRTCOD) C C Bring in include file with integer control variables. IMPLICIT NONE INCLUDE (OSLI) C REAL*8 DSPACE(*),DARRAY(5),DVAL,DBEST,DBEST2 INTEGER*4 MSPACE(*),MARRAY(6),M2,IREASN,JRTCOD,RTCOD,ISEQ, + ISET,IF9 SAVE DBEST,DBEST2,ISEQ,ISET,IF9,DVAL C C First check to see if this problem has SOS -- not valid with C this example user exit. C M2 = MARRAY(2) IF((M2.EQ.1).OR.(M2.EQ.2).OR.(M2.EQ.3)) THEN PRINT *,'THIS EKKBRNU NOT FOR SOS' PRINT *,'STOPPING YOUR APPLICATION NOW' STOP ENDIF C GOTO (1000,2000,3000,4000,5000),IREASN C C Initialization C Using 0.499999 so that variables at 0.0 will be skipped C 1000 DBEST=0.499999D0 DBEST2=0.0D0 ISEQ=0 ISET=0 IF9=0 CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) WRITE (6,101) 'There are ',INUMSETS,' sets.' GOTO 6000 C C Before Set C 2000 WRITE (6,102) 'Set number ',MARRAY(1),MARRAY(2),MARRAY(3), + MARRAY(4) C C Set return code to force EKKBRNU calls. JRTCOD=0 GOTO 6000 C C In Set C Choose variable closest to 1.0 or to 0.5 if none above .9. C 3000 IF(DARRAY(1).GE.0.9D0) THEN IF(DARRAY(1).GT.1.000005D0) THEN PRINT *,'THIS EKKBRNU NOT FOR GENERAL INTEGER VARIABLES' PRINT *,'STOPPING YOUR APPLICATION NOW' STOP ENDIF IF(DARRAY(1).GT.DBEST2.AND.DARRAY(1).LE.0.99999D0) THEN ISEQ=MARRAY(5) DBEST2=DARRAY(1) DVAL=DARRAY(1) C Branch up. MARRAY(6)=1 ISET=MARRAY(1) IF9=1 ENDIF ELSEIF(ABS(DARRAY(1)-0.5D0).LT.DBEST.AND.IF9.EQ.0) THEN ISEQ=MARRAY(5) DBEST=ABS(DARRAY(1)-0.5D0) DVAL=DARRAY(1) C Branch up. MARRAY(6)=1 ISET=MARRAY(1) ENDIF GOTO 6000 C C At end of set -- point to the chosen one. C 4000 WRITE (6,103) 'OSL chooses ',MARRAY(5),MARRAY(6),DARRAY(1) WRITE (6,104) 'I chose ',ISEQ,1,DVAL IF(ISET.EQ.MARRAY(1)) THEN MARRAY(5)=-ISEQ ENDIF GOTO 6000 C 5000 MARRAY(1)=-ISET 6000 CONTINUE RETURN 101 FORMAT(1X,A10,I5,A5) 102 FORMAT(1X,A11,I5,I5,I5,I5) 103 FORMAT(1X,A11,I5,I5,F9.2) 104 FORMAT(1X,A8,I5,I5,F9.2) END
C EKKCHNU is called while EKKMSLV is in the process of choosing a node: C C No nodes have yet been processed. . CALL EKKCHNU(DSPACE,MSPACE, 1 ,IARRAY,RARRAY,USERRTCD) C Here you could set USERRTCD to have EKKCHNU called C for each active node (i.e. REASON=2 below). . C Loop over all the active nodes. DO 100 I = 1, "the number of active nodes" . . . CALL EKKCHNU(DSPACE,MSPACE, 2 ,IARRAY,RARRAY,USERRTCD) C Here you might do your own node choice evaluation. 100 CONTINUE . CALL EKKCHNU(DSPACE,MSPACE, 3 ,IARRAY,RARRAY,USERRTCD) C Here you may change the default choice of node.
C*********************************************************************** C C EXCHNU C C This user exit subroutine chooses the next node for the branch C and bound algorithm. C C IREASN=1 : Subroutine is being called before going through the list C of active nodes. C =2 : Subroutine is being called as each active node is being C processed. C =3 : Subroutine is being called after going through the list C of active nodes. C C MARRAY(1): If IREASN=1, MARRAY(1) is the number of active nodes; C the rest of MARRAY and DARRAY are not meaningful. C If IREASN=2, MARRAY(1) is the number of the current node. C If IREASN=3, MARRAY(1) is the number of the chosen node; C the rest of MARRAY and DARRAY are not meaningful. C C MARRAY(2): Number of unsatisfied variables. C C DARRAY(2): Solution value of parent node. C C*********************************************************************** C SUBROUTINE EKKCHNU(DSPACE,MSPACE,IREASN,MARRAY,DARRAY,JRTCOD) C C Bring in include file with real control variable definitions. IMPLICIT NONE INCLUDE (OSLR) C REAL*8 DSPACE(*),DARRAY(4),RXTARGET,DBEST,DVALOBJ INTEGER*4 MSPACE(*),IREASN,MARRAY(6),JRTCOD,RTCOD,NUMBER,IBEST SAVE DBEST,IBEST,NUMBER,RXTARGET C C Communication with EXEVNU(EKKEVNU). COMMON/USERCOM/DVALOBJ C C Choose node with smallest number of infeasibilities, C but try to stay above the target value. C GOTO (1000,2000,3000) ,IREASN C C Initialization C 1000 DBEST=1.0D30 IBEST=0 NUMBER=999999 C C Find target value. CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) RXTARGET=RTARGET C C Set return code to force EKKCHNU calls. JRTCOD=0 GOTO 4000 C C Ordinary call. C 2000 IF(DARRAY(2).LT.RXTARGET) THEN C Better than the target, so choose one with fewest C infeasibilities. IF(NUMBER.GT.MARRAY(2)) THEN C Fewer infeasibilities -- choose this one. DBEST=DARRAY(2) IBEST=MARRAY(1) NUMBER=MARRAY(2) ELSEIF(NUMBER.EQ.MARRAY(2).AND.DBEST.GT.DARRAY(2)) THEN C Same number, but better value. DBEST=DARRAY(2) IBEST=MARRAY(1) ENDIF ENDIF GOTO 4000 C C Done C If any node better than target, point to one with fewest C infeasibilities. C 3000 IF(IBEST.NE.0) MARRAY(1)=IBEST C C Store value of objective for EKKEVNU. DVALOBJ=DBEST 4000 CONTINUE C RETURN END
SUBROUTINE EKKEVNU(DSPACE,MSPACE,JRTCOD) INCLUDE (OSLI) REAL*8 DSPACE(*) INTEGER*4 MSPACE(*) INTEGER*4 JRTCOD CALL EKKSSLV(JRTCOD,DSPACE,2,1) CALL EKKIGET(IRTCOD,DSPACE,OSLI,OSLILN) C Set the user return code to 1 if the optimal solution was not found. IF(IPROBSTAT.NE.0) JRTCOD = 1 RETURN ENDThe following is a alternative, sample EKKEVNU user exit subroutine.
C*********************************************************************** C C EXEVNU C C This user exit routine determines whether the primal or the C dual simplex method will be used to evaluate the current node. C C The routine is written under the assumption that for the class C problems being solved, the primal algorithm is faster than the C dual algorithm. However, with the dual algorithm, the branch can C be cut off when the objective increases above the cutoff. So this C routine selects the primal algorithm unless the objective is near C the cutoff. It is assumed that EKKCHNU saved the value of the C objective function in DVALOBJ. C C Note that the return code must be set if the problem in infeasible. C C*********************************************************************** C SUBROUTINE EKKEVNU(DSPACE,MSPACE,JRTCOD) C C Bring in include files with control variable definitions. IMPLICIT NONE INCLUDE (OSLR) INCLUDE (OSLI) C REAL*8 DSPACE(*),DVALOBJ INTEGER*4 MSPACE(*),RTCOD,JRTCOD C C Communication with EXCHNU(EKKCHNU). COMMON/USERCOM/DVALOBJ C C Get cutoff. CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) IF(DVALOBJ+2000.0.GT.RBBCUTOFF) THEN C Near cutoff, so use dual algorithm. CALL EKKSSLV(RTCOD,DSPACE,2,1) ELSE C Not near cutoff, so use primal algorithm. CALL EKKSSLV(RTCOD,DSPACE,1,1) ENDIF CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) C C Set user return code if needed. IF(IPROBSTAT.NE.0) THEN JRTCOD = 1 ELSE JRTCOD = 0 ENDIF C RETURN END
To use EKKSLVU you need to know when and how the user exit subroutine is called. Pseudo-code for the call to EKKBRNU is given below.
The sample user exit program
EXSLVU, in FORTRAN,
shows how EKKSLVU can be used to solve a quadratic MIP problem.
C EKKSLVU is called to specify a solver for each node C of the branch and bound tree; C No nodes have been solved. . . C Make first call to solve the root node CALL EKKSLVU(DSPACE,MSPACE,DOBJVAL,JRTCOD, 1 ) . C Proceed with branch and bound analysis DO WHILE UNSOLVED NODES EXIST . . CALL EKKSLVU(DSPACE,MSPACE,DOBJVAL,JRTCOD, 2 ) . . ENDDO C C Make a final call to re-solve the optimal node CALL EKKSLVU(DSPACE,MSPACE,DOBJVAL,JRTCOD, 3 )
To use EKKNODU you need to know when and how the user exit subroutine
is called. Following is pseudo-code for the call to EKKNODU and a sample
EKKNODU user exit subroutine.
Solve the initial problem and other initialize structures. LOOP: Take the solved problem off list of active nodes. IF (solution possible(i.e. feasible and good value)) THEN Decide what variable would be branched on (at the end of this EKKBRNU is called) IF(integer solution or user says integer solution) THEN Action for solution e.g. change real control variable Rbbcutoff ELSE Call EKKNODU(rtcod,dspace,1...) to obtain any information the user wants to save about this node. The user can delete branching information. Save branching information, basis and EKKNODU information on the list of active nodes. ENDIF ENDIF Choose one node from the list of active nodes. (at the end of this EKKCHNU is called) IF(any active nodes) THEN Restore branching information (set bounds), basis and EKKNODU information Call EKKNODU(rtcod,dspace,2,...) so the user can use information which was saved to set the bounds in a more flexible way. Evaluate node (call EKKEVNU) IF (supernode processing) THEN action similar to EKKMPRE ENDIF GOTO LOOP ENDIF The search is finished.
C*********************************************************************** C C EXMSLV3 C C This driver finds the minimum value (approximately) of a polynomial C of degree 3 in two variables. x and y may vary between -2 and +2. C Although this example is linear with a nonlinear objective, the same C method may be applied with nonlinear constraints and the problem C need not be convex or even continuous. This driver uses a piece-wise C linear approximation of the objective function. It is intended for C use with the user exit sample subroutines EXNODU and EXBRNU2. C C*********************************************************************** C PROGRAM MAIN C C Bring in include files with control variable definitions. IMPLICIT REAL*8 (D) INCLUDE (OSLR) INCLUDE (OSLI) INCLUDE (OSLN) C C Allocate dspace and other arrays. PARAMETER (MAXSPC=1000000) REAL*8 DSPACE(MAXSPC) INTEGER*4 MSPACE(2*MAXSPC),RTCOD EQUIVALENCE(DSPACE,MSPACE) COMMON/BIG/DSPACE C C Coefficients for each term. C DCOEFF(1,0) is coefficient of linear X term etc. PARAMETER (NDEGREE=3) REAL*8 DCOEFF(0:NDEGREE,0:NDEGREE) C C Number of grid points for x and y; the example below is every 0.1 PARAMETER (NGRID=41) C We need NGRID for x and 2*NDEGREE*NGRID for y + 2 for result. PARAMETER (NCOL=NGRID+2*NDEGREE*NGRID+2) C We need two rows for x and 3*NDEGREE rows for y. PARAMETER (NROW=2+3*NDEGREE) C Arrays for bounds and costs. REAL*8 DRLO(NROW),DRUP(NROW) REAL*8 DCLO(NCOL),DCUP(NCOL),DCOST(NCOL) REAL*8 XVAL,YVAL C Array to generate names (not necessary but makes it easier C for user to see what is happening). CHARACTER*12 CROW(NROW),CCOL(NCOL) C Arrays for SOS variables and sets (one set per row). REAL*8 DREFRW(NCOL) C We can use same array for two purposes in IMDL. INTEGER*4 MSEQ(NCOL),MPRI(NROW),MSETS(NROW+1),MTYPE(NROW) C C This example is 2 X**2 Y + X Y**2 + 25 X**2 +100 Y**2 C + 20X - 200Y DATA DCOEFF/ 0.0D0, 2.0D1, 2.5D1, 0.0D0, + -2.0D2, 0.0D0, 2.0D0, 0.0D0, + 1.0D2, 1.0D0, 0.0D0, 0.0D0, + 0.0D0, 0.0D0, 0.0D0, 0.0D0/ C C Describe application and specify that there is 1 model. CALL EKKDSCA(RTCOD,DSPACE,MAXSPC,1) IF (RTCOD.GT.0) CALL CHKRT('EKKDSCA',RTCOD) C C Describe the model as having 1 block. CALL EKKDSCM(RTCOD,DSPACE,1,1) IF (RTCOD.GT.0) CALL CHKRT('EKKDSCM',RTCOD) C C Build up matrix. Put in as triplets directly into dspace. CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN) IF (RTCOD.GT.0) CALL CHKRT('EKKNGET',RTCOD) C C Compute maximum number of elements and where they may go. NELMAX=(NLASTFREE-NFIRSTFREE)/2 IE=NFIRSTFREE IC=IE+NELMAX C Will be using mspace so change. IC=2*IC-1 IR=IC+NELMAX CALL MATRIX(DSPACE(IE),MSPACE(IC),MSPACE(IR),CROW,CCOL, + DRLO,DRUP,DCLO,DCUP,DCOST,DCOEFF,NGRID,NELMAX,NEL) C C Load model into dspace. OSL will know that some arrays already C in dspace and move them down to low end of dspace. CALL EKKLMDL(RTCOD,DSPACE,1,NROW,NCOL,NEL,DCOST,DRLO,DRUP, + DCLO,DCUP,MSPACE(IR),MSPACE(IC),DSPACE(IE)) IF (RTCOD.GT.0) CALL CHKRT('EKKLMDL',RTCOD) C C Make column copy. CALL EKKNWMT(RTCOD,DSPACE,2) IF (RTCOD.GT.0) CALL CHKRT('EKKNWMT',RTCOD) C C Put in names - 12 characters CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKIGET',RTCOD) INUMCHAR=12 CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKSGET',RTCOD) CALL EKKNAME(RTCOD,DSPACE,NROW,CROW,1,NCOL,CCOL,1,0) IF (RTCOD.GT.0) CALL CHKRT('EKKNAME',RTCOD) C C Put variables into ordinary S1 sets. C (S2 might give slightly better answer) C Build up arrays - doing x and y in same loop. ICOL=0 DO ISET=1,2*NDEGREE+1 MTYPE(ISET)=1 C Do X first and then in degree order. MPRI(ISET)=ISET MSETS(ISET)=ICOL+1 C Put in correct reference entry (value of x or y) C (could use 1,2 etc BUT this way we can use irregular grid). DREF=-2.0D0 DINCR=(2.0D0-(-2.0D0))/DFLOAT(NGRID-1) DO I=1,NGRID ICOL=ICOL+1 MSEQ(ICOL)=ICOL DREFRW(ICOL)=DREF DREF=DREF+DINCR ENDDO ENDDO C fill in last MSETS MSETS(NROW+1)=ICOL+1 C C Describe Integer structure (default down pseudocosts). C (last two variables are x and y - continuous) CALL EKKIMDL(RTCOD,DSPACE,NCOL-2,MSEQ,1+2*NDEGREE,MTYPE, + MPRI,NCOL-2,MSETS,MSEQ,0.0D0,DREFRW) IF (RTCOD.GT.0) CALL CHKRT('EKKIMDL',RTCOD) C C Solve. CALL EKKSSLV(RTCOD,DSPACE,1,2) IF (RTCOD.GT.0) CALL CHKRT('EKKSSLV',RTCOD) C C Print solution (selectively). CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKIGET',RTCOD) ISOLMASK=6 IMIPLENGTH=20 CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN) IF (RTCOD.GT.0) CALL CHKRT('EKKISET',RTCOD) CALL EKKPRTS(RTCOD,DSPACE) IF (RTCOD.GT.0) CALL CHKRT('EKKPRTS',RTCOD) C C Branch and Bound (The user may find it instructive to run this C without EKKNODU and see what happens). CALL EKKMSLV(RTCOD,DSPACE,1,0,0) IF (RTCOD.GT.0) CALL CHKRT('EKKMSLV',RTCOD) C C Print solution again CALL EKKPRTS(RTCOD,DSPACE) IF (RTCOD.GT.0) CALL CHKRT('EKKPRTS',RTCOD) C C Print values of x and y CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN) IF (RTCOD.GT.0) CALL CHKRT('EKKNGET',RTCOD) XVAL=DSPACE(NCOLSOL+NCOL-2) YVAL=DSPACE(NCOLSOL+NCOL-1) CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) IF (RTCOD.GT.0) CALL CHKRT('EKKRGET',RTCOD) DO I=0,NDEGREE PRINT 1,(DCOEFF(J,I),J=0,NDEGREE) 1 FORMAT(6F12.4) ENDDO PRINT 2,XVAL,YVAL,ROBJVALUE 2 FORMAT(' Minimum at X=',F12.4,' ,Y=',F12.4,' Objective',F20.4) STOP END C C*********************************************************************** C Generate matrix C*********************************************************************** C SUBROUTINE MATRIX(DELS,MCOL,MROW,CROW,CCOL, + DRLO,DRUP,DCLO,DCUP,DCOST,DCOEFF,NGRID,NELMAX,NEL) * ---------------------------------------------------- IMPLICIT REAL*8 (D) INCLUDE (OSLR) INCLUDE (OSLI) INCLUDE (OSLN) C Coefficients for each term PARAMETER (NDEGREE=3) REAL*8 DCOEFF(0:NDEGREE,0:NDEGREE) C Arrays for bounds and costs REAL*8 DRLO(*),DRUP(*) REAL*8 DCLO(*),DCUP(*),DCOST(*) C Array to generate names (not necessary but makes it easier C for user to see what is happening) CHARACTER*12 CROW(*),CCOL(*),CTEMP INTEGER*4 MCOL(*),MROW(*) REAL*8 DELS(*) C We need NGRID for x and 2*NDEGREE*NGRID for y NCOL=NGRID+2*NDEGREE*NGRID+2 C We need two rows for x and 3*NDEGREE rows for y NROW=2+3*NDEGREE C Fill in all column bounds DO ICOL=1,NCOL-2 DCLO(ICOL)=0.0D0 DCUP(ICOL)=1.0D31 ENDDO C Make two result variables free DCLO(NCOL-1)=-1.0D31 DCUP(NCOL-1)=1.0D31 DCOST(NCOL-1)=0.0D0 DCLO(NCOL)=-1.0D31 DCUP(NCOL)=1.0D31 DCOST(NCOL)=0.0D0 C two rows for X - convexity row and result DRLO(1)=1.0D0 DRUP(1)=1.0D0 CROW(1)='X_convexity' DRLO(2)=0.0D0 DRUP(2)=0.0D0 DREFX=-2.0D0 CROW(2)='X_value_row' DINCRX=(2.0D0-(-2.0D0))/DFLOAT(NGRID-1) NEL=0 DO I=1,NGRID C Convexity entry is 1 NEL=NEL+1 MROW(NEL)=1 MCOL(NEL)=I DELS(NEL)=1.0D0 C result NEL=NEL+1 MROW(NEL)=2 MCOL(NEL)=I DELS(NEL)=DREFX C Cost - X only term DTOT=DCOEFF(0,0) DVAL=1.0D0 DO J=1,NDEGREE DVAL=DVAL*DREFX DTOT=DTOT+DVAL*DCOEFF(J,0) ENDDO DCOST(I)=DTOT WRITE(CTEMP,1) DREFX 1 FORMAT('Xat',F5.2) IF(DREFX.GT.-1.0D-6) CTEMP(4:4)='+' CCOL(I)=CTEMP DREFX=DREFX+DINCRX ENDDO C and put in X NEL=NEL+1 MROW(NEL)=2 MCOL(NEL)=NCOL-1 CCOL(NCOL-1)='X_value' DELS(NEL)=-1.0D0 CCOL(NCOL)='Y_value' C Now Y sets - two sets per degree of X and three rows IROW=2 ICOL=-NGRID DO IDEG=1,NDEGREE ICONV=IROW+1 WRITE(CTEMP,2) IDEG 2 FORMAT('Convexity_',I2.2) CROW(ICONV)=CTEMP IYROW=IROW+2 WRITE(CTEMP,3) IDEG 3 FORMAT('Y_value_',I2.2) CROW(IYROW)=CTEMP IXYROW=IROW+3 WRITE(CTEMP,4) IDEG 4 FORMAT('X_power_',I2.2) CROW(IXYROW)=CTEMP DRLO(ICONV)=1.0D0 DRUP(ICONV)=1.0D0 IROW=IROW+3 C base for max ICOL=ICOL+2*NGRID ICOL2=ICOL+NGRID DREFY=-2.0D0 DINCRY=(2.0D0-(-2.0D0))/DFLOAT(NGRID-1) C lower and upper limits on x**ideg IF(IAND(IDEG,1).EQ.0) THEN DLOVAL=0.0D0 DHIVAL=2.0D0**IDEG ELSE DLOVAL=-2.0D0**IDEG DHIVAL=2.0D0**IDEG ENDIF C And point back to correct power of x DREFX=-2.0D0 DINCRX=(2.0D0-(-2.0D0))/DFLOAT(NGRID-1) DO I=1,NGRID NEL=NEL+1 MROW(NEL)=IXYROW MCOL(NEL)=I DELS(NEL)=-DREFX**IDEG DREFX=DREFX+DINCRX ENDDO DO I=1,NGRID C Convexity NEL=NEL+1 MROW(NEL)=ICONV MCOL(NEL)=I+ICOL WRITE(CTEMP,5) 'L_',IDEG,DREFY C note on naming conventions C X_nn gives power of X involved C L_ or U_ refers to the minimum or maximum value of X_nn C last term gives value of Y 5 FORMAT(A2,'X_',I2.2,'_',F5.2) IF(DREFY.GT.-1.0D-6) CTEMP(8:8)='+' CCOL(I+ICOL)=CTEMP DELS(NEL)=1.0D0 NEL=NEL+1 MROW(NEL)=ICONV MCOL(NEL)=I+ICOL2 CTEMP(1:1)='U' CCOL(I+ICOL2)=CTEMP DELS(NEL)=1.0D0 C result NEL=NEL+1 MROW(NEL)=IYROW MCOL(NEL)=I+ICOL DELS(NEL)=DREFY NEL=NEL+1 MROW(NEL)=IYROW MCOL(NEL)=I+ICOL2 DELS(NEL)=DREFY C link so that we have x to the correct power NEL=NEL+1 MROW(NEL)=IXYROW MCOL(NEL)=I+ICOL DELS(NEL)=DLOVAL NEL=NEL+1 MROW(NEL)=IXYROW MCOL(NEL)=I+ICOL2 DELS(NEL)=DHIVAL C Cost for X**IDEG and all Y costs DTOT=0.0D0 DVAL=1.0D0 DO J=1,NDEGREE DVAL=DVAL*DREFY DTOT=DTOT+DVAL*DCOEFF(IDEG,J) ENDDO DCOST(I+ICOL)=DLOVAL*DTOT DCOST(I+ICOL2)=DHIVAL*DTOT C add in Y only part of objective on X**1 set IF(IDEG.EQ.1) THEN DTOT=0.0D0 DVAL=1.0D0 DO J=1,NDEGREE DVAL=DVAL*DREFY DTOT=DTOT+DVAL*DCOEFF(0,J) ENDDO DCOST(I+ICOL)=DCOST(I+ICOL)+DTOT DCOST(I+ICOL2)=DCOST(I+ICOL2)+DTOT ENDIF DREFY=DREFY+DINCRY ENDDO C and put in Y NEL=NEL+1 MROW(NEL)=IYROW MCOL(NEL)=NCOL DELS(NEL)=-1.0D0 ENDDO C C Check space availability. IF(NEL.GT.NELMAX) THEN WRITE(6,*)'Increase size of DSPACE and recompile program' WRITE(6,*)'Stopping your application now' STOP ENDIF RETURN END C C*********************************************************************** C This subroutine prints the character string RTNAME and the return C code RTCOD and stops if RTCOD is large enough to indicate that an C OSL error or severe error has occurred. C*********************************************************************** C SUBROUTINE CHKRT(RTNAME,RTCOD) CHARACTER*7 RTNAME INTEGER*4 RTCOD C WRITE(6,9000) RTNAME,RTCOD IF (RTCOD.GE.200) STOP 16 RETURN 9000 FORMAT (1X,'********** ',A7,' return code of ',I4,' **********') END
C*********************************************************************** C C EXNODU C C This user exit routine allows for non standard branching. C Data will be saved and restored to allow very flexible C branch and bound. It is intended for use with EXMSLV3 and C EXBRNU2. C C DSPACE is the user work area. C MSPACE is the user work area. C IREASN is the reason for calling. C MARRAY is the short node list -- integer values. C DARRAY is the short node list -- real*8 values (equivalenced). C C MARRAY(1) - Branch type SOS=1-3, Integer=4,Null=99 C MARRAY(2) - 0 for down branch , 1 for up C MARRAY(3) - Variable number / Set number C DARRAY(1) - Integer value/Set reference value C C MUSER - User array C NUSER - Length of array C C*********************************************************************** C SUBROUTINE EKKNODU(DSPACE,MSPACE,IREASN,MARRAY,DARRAY,MUSER,NUSER) C C Include files with control variable definitions. IMPLICIT REAL*8(D) INCLUDE (OSLI) INCLUDE (OSLN) C REAL*8 DSPACE(*),DARRAY(1),DTEMP(10) INTEGER*4 MARRAY(3),MUSER(NUSER),MSPACE(*),MTEMP(20) EQUIVALENCE (MTEMP,DTEMP) SAVE DTEMP,MTEMP C C Communication with EXBRNU2(EKKBRNU) COMMON/USRCOM1/DBRANCH C C This example assumes Imiplength was 20 IF(NUSER.NE.20) THEN WRITE(6,*)'Error detected in user exit EKKNODU' WRITE(6,*)'Bad example coding - Imiplength must be = 20' WRITE(6,*)'Stopping your application now' STOP ENDIF C GOTO (1000,2000) ,IREASN C C Return if ordinary branch 1000 IF(MARRAY(3).EQ.1) RETURN C C Save any information which will be needed when node is evaluated. C this information can be of any kind. The user could be C continually refining the grid, in which case information C on mesh size and X,Y values of basic variables would be stored C then in EKKEVNU the matrix could be completely rewritten! C DTEMP(1)=DBRANCH DO I = 1,NUSER MUSER(I)=MTEMP(I) ENDDO C Switch off ***ANY*** action by OSL to fix variables MARRAY(1)=99 GOTO 3000 C C Return if ordinary branch 2000 IF(MARRAY(1).NE.99) RETURN C Now set all variables to zero DO I = 1,NUSER MTEMP(I)=MUSER(I) ENDDO DBRANCH=DTEMP(1) IF(MARRAY(2).EQ.1) THEN C Up branch - fix all low reference values to zero DLO=-1.0D20 DHI=DBRANCH ELSE C Down branch - fix all high reference values to zero DLO=DBRANCH DHI=1.0D20 ENDIF CALL EKKIGET(IRTCOD,DSPACE,OSLI,OSLILN) CALL EKKNGET(IRTCOD,DSPACE,OSLN,OSLNLN) C C Number in grid NGRID=INUMINTS/INUMSETS IUP=NCOLUPPER-1+NGRID C C Through other sets getting *linked* reference entry C (In real use we would use true reference entries for grid) NCENTER= 1+(NGRID-1)/2 DPERGRID=2.0D0/DFLOAT((NGRID-1)/2) DO ISET =2,INUMSETS DO IVAR=1,NGRID IUP=IUP+1 DBRANCH=DFLOAT(IVAR-NCENTER)*DPERGRID C zero bound if in range IF(DBRANCH.GE.DLO.AND.DBRANCH.LE.DHI) THEN DSPACE(IUP)=0.0D0 ENDIF ENDDO ENDDO 3000 RETURN END
C*********************************************************************** C C EXBRNU2 C C This user exit routine modifies the default choice of the branching C variable. It is intended for use with EXMSLV3 and EXNODU. C C DSPACE is the user work area. C MSPACE is the user work area. C IREASN is the reason for calling. C MARRAY is the short node list -- integer values. C DARRAY is the short node list -- real*8 values (equivalenced). C JRTCOD is the return code. C C MARRAY(1) is the set number. (5) (negative allowed) C MARRAY(2) is the set type. C MARRAY(3) is the priority. (2) C MARRAY(4) is the number of variables in the set. C MARRAY(5) is the column number of current variable. (4) C (negative allow) C MARRAY(6) is which way. (4) C DARRAY(1) is the current value of variable. C DARRAY(2) is the down pseudo-cost. (3) C DARRAY(3) is the up cost (3) or ref row entry (4). C DARRAY(4) is the estimated degradation for going down. (4) (5) C DARRAY(5) is the estimated degradation for going up. (4) (5) C C Note: C The (n) to the right of the definitions of the MARRAY and DARRAY C elements indicate that these data items are available and/or may C be set when the calling parameter REASON=n. C C*********************************************************************** C SUBROUTINE EKKBRNU(DSPACE,MSPACE,IREASN,MARRAY,DARRAY,JRTCOD) C C Include files with control variable definitions. IMPLICIT REAL*8(D) INCLUDE (OSLI) INCLUDE (OSLN) C C Communication with EXNODU(EKKNODU). COMMON/USRCOM1/DBRANCH C REAL*8 DSPACE(*),DARRAY(5) INTEGER*4 MSPACE(*),MARRAY(6) SAVE DREFSUM,DSUM,DLO,DHI C C Only interested on first and last call GOTO (1000,6000,6000,6000,5000) ,IREASN C C Initial entry - do all work here 1000 CALL EKKIGET(IRTCOD,DSPACE,OSLI,OSLILN) CALL EKKNGET(IRTCOD,DSPACE,OSLN,OSLNLN) WRITE (6,101) 'There are ',INUMSETS,' sets and',INUMINTS, + ' variables' 101 FORMAT(1X,A,I5,A,I5,A) C C Number in grid NGRID=INUMINTS/INUMSETS ISOL=NCOLSOL-1 C C First set is X by itself NNZERO=0 DO IVAR=1,NGRID ISOL=ISOL+1 DVAL=DSPACE(ISOL) IF(ABS(DVAL).GT.1.0D-5) NNZERO=NNZERO+1 ENDDO C C Skip if first set unsatisfied as can do ordinary branch IF(NNZERO.GT.1) GO TO 6000 C C Through other sets getting *linked* reference entry C (In real use we would use true reference entries for grid) NCENTER= 1+(NGRID-1)/2 DPERGRID=2.0D0/DFLOAT((NGRID-1)/2) DREFSUM=0.0D0 DSUM=0.0D0 DLO=1.0D20 DHI=-1.0D20 DO ISET =2,INUMSETS DO IVAR=1,NGRID ISOL=ISOL+1 DVAL=DSPACE(ISOL) C accumulate IF(ABS(DVAL).GT.1.0D-5) THEN DREFVAL=DFLOAT(IVAR-NCENTER)*DPERGRID DREFSUM=DREFSUM+DVAL*DREFVAL DSUM=DSUM+DVAL DLO=MIN(DLO,DREFVAL) DHI=MAX(DHI,DREFVAL) ENDIF ENDDO ENDDO GOTO 6000 C C At end 5000 IF(MARRAY(1).EQ.1) RETURN C X satisfied - check Y IF(DHI-DLO.LE.1.0D-5) RETURN C choose set 2 MARRAY(1)=-2 C choose where to branch DBRANCH=DREFSUM/DSUM 6000 CONTINUE RETURN END
As with all the examples in this chapter, these samples are not intended
to be used without alteration. However, unlike the branch-and-bound user
exit subroutines EKKBRNU, EKKCUTU, and EKKEVNU, these user exit subroutines
can easily extend into hundreds of lines of code. Since these user exit
subroutines are also extremely problem dependent, it is difficult to provide
useful samples. Therefore, the samples that are provided are intended for
illustrative purposes only.
The following is a sample EKKCUTU user exit subroutine.
C*********************************************************************** C C EXCUTU C C Sample user exit routine EKKCUTU. EKKCUTU allows you to generate C "cuts", or additional constraint rows, to a problem matrix. It is C called during EKKMPRE preprocessing, or during supernode C processing in EKKMSLV. C C On Input MSTATBIN tells the user which 0-1 variables have been C fixed so far, -1 says fixed to 0, +1 to 1, 0 still free. C It is of length NBIN. MTOBIN is an array of INUMCOLS C length which has 0 for a continuous variable, the type for C other "Integer" variables, and points into MSTATBIN for all C 0-1 variables. C C The user returns a series of cuts by filling in MCADD,MRADD, C DEADD and NADD. C C MCADD - 1 is the first column of matrix, INUMCOLS is last. C minus 1 denotes a lower bound on a row and minus 2 C denotes an upper bound. C C If on entry NADD were 1002 and NROW were 205 then to add a C single cut of X + 2* Y <= 2 might be: C C IADD MCADD MRADD DEADD C C 1003 5 206 1.0 C 1003 407 206 2.0 C 1003 -2 206 2.0 C C and NADD would be changed to 1005 on exit, and NROW to 206. C C DSPACE - Main Data Array C MSTATBIN- Column status: 0 - free, 1 - fixed to 1, -1 - fixed to 0 C MTOBIN - 0 or 0-1 sequence for each real variable C MCADD - Column indices of cuts C MRADD - Row indices of cuts C DEADD - Elements of cuts C NBIN - Number of 0-1 variables in MSTAT01 C NADD - Current number of entries in MCADD,MRADD and DEADD C NADDMAX - Maximum number of entries allowed C NROW - Current number of rows in matrix (including cuts) C NROWMAX - Maximum number of rows allowed C C*********************************************************************** C SUBROUTINE EKKCUTU(DSPACE,MSTATBIN,MTOBIN,MCADD,MRADD,DEADD,NBIN, + NADD,NADDMAX,NROW,NROWMAX) C IMPLICIT REAL*8 (D,Z) REAL*8 DSPACE(*) INTEGER*4 MSTATBIN(NBIN),MTOBIN(NBIN) INTEGER*4 MCADD(NADDMAX),MRADD(NADDMAX) REAL*8 DEADD(NADDMAX) C C In actual use one or more cuts would be computed. C For explanations we assume that variables are X1,X2.... C Assume first cut is magically in a Data statement REAL*8 DELEMENT(5),DRHS INTEGER*4 MCOLUMN1(5),MCOLUMN2(3) C C First cut will be 5 X4 - 7 X6 + X8 +10 X30 - X2 <= 8 DATA DELEMENT/5.0D0,-7.0D0,1.0D0,10.0D0,-1.0D0/,DRHS/8.0D0/ DATA MCOLUMN1/4,6,8,30,2/ DATA MCOLUMN2/1,4,7/ C C Add in first cut if enough room IF(NROW.LT.NROWMAX.AND.NADD.LE.NROWMAX-6) THEN NROW=NROW+1 C elements DO I=1,4 NADD=NADD+1 MRADD(NADD)=NROW MCADD(NADD)=MCOLUMN1(I) DEADD(NADD)=DELEMENT(I) ENDDO C upper bound on row activity NADD=NADD+1 MRADD(NADD)=NROW MCADD(NADD)=-2 DEADD(NADD)=DRHS ENDIF C C Second cut will be X1 + X4 + X7 = 1 IF(NROW.LT.NROWMAX.AND.NADD.LE.NROWMAX-5) THEN NROW=NROW+1 C elements DO I=1,3 NADD=NADD+1 MRADD(NADD)=NROW MCADD(NADD)=MCOLUMN2(I) DEADD(NADD)=1.0D0 ENDDO C upper bound on row activity NADD=NADD+1 MRADD(NADD)=NROW MCADD(NADD)=-2 DEADD(NADD)=1.0D0 C and lower bound NADD=NADD+1 MRADD(NADD)=NROW MCADD(NADD)=-1 DEADD(NADD)=1.0D0 ENDIF RETURN END
The following is a sample EKKHEUU user exit subroutine.
C*********************************************************************** C C EXHEUU C C Sample user exit routine EKKHEUU. EKKHEUU gives you control over C the fixing of variables to their upper or lower bounds. C C On Input MSTATBIN tells the user which 0-1 variables have been C fixed so far, -1 says fixed to 0, +1 to 1, 0 still free. C It is of length NBIN. MTOBIN is an array of INUMCOLS C length which has 0 for a continuous variable, the type for C other "Integer" variables, and points into MSTATBIN for all C 0-1 variables. C C The user returns a valid branch. OSL will continue by C setting 0-1 variables to their bounds using the first C NFIX1 entries in MFIX. The postponed branch (which is only C used if JTYPE is 2) is given by the next NFIX2 entries of MFIX. C C In this example the user knows that going up to 1 is always C feasible, so the variable closest to one will be chosen C C Calling sequence - CALL EKKHEUU C C DSPACE - Main Data Array C MSTATBIN - Column status: 0 -free, 1 -fixed to 1, -1 -fixed to 0 C MTOBIN - 0 or 0-1 sequence for each real variable C MFIX - Stack of 0-1s - negative fixed to, positive to 1 C (maximum length 2*NBIN) C NBIN - Number of 0-1 variables in MSTATBIN C NFIX1 - Number of variables fixed on this way C NFIX2 - Number of variables fixed on postponed way C JTYPE - 1 for heuristic guess, 2 for heuristic branch C C********************************************************************* C SUBROUTINE EKKHEUU(DSPACE,MSTATBIN,MTOBIN,MFIX,NBIN,NFIX1, + NFIX2,JTYPE) C IMPLICIT REAL*8 (D,Z) REAL*8 DSPACE(*) INTEGER*4 MSTATBIN(NBIN),MFIX(*),MTOBIN(*) C C Bring in include files with control variable definitions. INCLUDE (OSLN) INCLUDE (OSLI) C C Do nothing unless in Branch and bound phase IF(JTYPE.NE.2) RETURN C C Get all values (GETs need only be done once) CALL EKKIGET(IRET,DSPACE,OSLI,OSLILN) CALL EKKNGET(IRET,DSPACE,OSLN,OSLNLN) C C Find largest 0-1 value (not at 1) DLARGE=0.0D0 ICHOSEN=0 DO ICOL=1,INUMCOLS IF(MTOBIN(ICOL).GT.0) THEN DVAL=DSPACE(NCOLSOL+ICOL-1) IF(DVAL.LT.0.9999D0.AND.DVAL.GT.DLARGE) THEN DLARGE=DVAL ICHOSEN=ICOL ENDIF ENDIF ENDDO IF(ICHOSEN.EQ.0) STOP 'This should never happen' C C Find 0-1 sequence ICHOSEN=MTOBIN(ICHOSEN) C C Carry on with this one to one NFIX1=1 MFIX(1)=ICHOSEN C C And postpone other way NFIX2=1 MFIX(2)=-ICHOSEN RETURN END