C ALGORITHM 666, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 14, NO. 4, PP. 330-336.
PROGRAM MAINSP
***********************************************************************
* *
* This is a sample driver program for our mathematical software *
* package CHABIS (CHAracteristic BISection) which solves systems *
* of nonlinear (algebraic and/or transcendental) equations of the *
* form *
* - - - *
* F ( X ) = 0 , *
* where *
* - T *
* F = ( f1, f2, ... , fn ) is a continuous mapping of a *
* bounded region in the real N-space into the real N-space. *
* *
* First, CHABIS locates at least one solution of the given system *
* within an N-dimensional Polyhedron. Then it obtains an approxi- *
* mate solution using a new generalized method of bisection. *
* *
* This driver program runs all the examples of the accompanying *
* paper. It echoes the example number and the input values (ini- *
* tial guess and stepsizes), and prints out the final approximate *
* root. Also, it prints out performance information such as that *
* given in the accompanying paper. It takes the input values from *
* the file SINPUT and provides the print out in the file OUTPUT. *
* *
* *
* CHABIS. Version of June 1988 *
* Michael N. Vrahatis *
* *
***********************************************************************
INTRINSIC INT, FLOAT
INTEGER NEXT, N, NFCALL, LIU, LOU, I, ICON, LWA, J, INF1, INF2
REAL DELTA, EPSILO, XO(9), H(9), AS(9), VAS(9)
REAL WA(28178)
CHARACTER*9 EXAMP
CHARACTER*70 F1, F2
EXTERNAL FNC
COMMON / BLK1 / NEXT, N, NFCALL
*
* The Logical Input Unit is assumed to be the number 7. *
*
DATA LIU / 7 /
*
* The Logical Output Unit is assumed to be the number 8. *
*
DATA LOU / 8 /
*
* Open the Input and Output files. *
*
OPEN ( UNIT = LIU, FILE = 'SINPUT', STATUS = 'OLD' )
OPEN ( UNIT = LOU, FILE = 'OUTPUT', STATUS = 'NEW' )
DO 10 I = 1, 24
*
* Set the value of DELTA. DELTA is a positive variable which *
* determines the accuracy of the computation of the roots, of *
* the components of the given function, which are located on *
* the edges of the Initial N-Polyhedron, and are used for the *
* construction of a Characteristic N-Polyhedron. Note that if *
* DELTA is less than the machine precision EPSMCH, or if DELTA *
* is not defined, the default value of DELTA becomes equal to *
* 0.0625 . The EPSMCH is computed within CHABIS. *
*
DELTA = 0.625E-1
*
* Set the value of EPSILO , which is a nonnegative variable. *
* Termination occurs when the algorithm estimates that the *
* Infinity Norm, of the function values at an approximate *
* solution, is at most EPSILO. If EPSILO is less than the *
* machine precision EPSMCH, or if EPSILO is not defined, the *
* default value of EPSILO becomes equal to EPSMCH. *
*
EPSILO = 1.0E-8
*
* Set the condition variable ICON. Note that if ICON is equal *
* to 1, then CHABIS attempts to find an approximate solution *
* even if the default Characteristic N-Polyhedron construction *
* fails. Otherwise, set ICON equal to another integer value. *
*
ICON = 1
*
* Set the starting values. *
*
READ (LIU, *) EXAMP
READ (LIU, *) F1, F2
READ (LIU, *) N
READ (LIU, *) ( XO(J), J = 1, N )
READ (LIU, *) ( H(J), J = 1, N )
READ (LIU, *) NEXT
*
* Set the length LWA of the workspace array WA. Note that LWA *
* must be at least equal to the value ( 2*N+(6*N+1)*2**N ). *
*
LWA = 2*N + (6*N+1) * 2**N
*
WRITE (LOU, 9999) EXAMP, F1, F2, N
WRITE (LOU, 9998) ( XO(J), H(J), J = 1, N )
*
* Call the interface subroutine INTSUB. *
*
NFCALL = 0
CALL INTSUB( FNC, N, XO, H, DELTA, EPSILO, ICON, INF1,
+ AS, VAS, INF2, WA, LWA)
*
IF ( INF1 .EQ. 0 ) THEN
WRITE (LOU, 9997)
ELSEIF ( INF1 .EQ. 2 .AND. ICON .NE. 1 ) THEN
WRITE (LOU, 9996) N
ELSE
WRITE (LOU, 9995) DELTA, EPSILO, (AS(J), VAS(J), J = 1, N)
ENDIF
NFCALL = INT( FLOAT( NFCALL ) / FLOAT( N ) )
WRITE (LOU, 9994) INF1, INF2, NFCALL
10 CONTINUE
STOP
*
* Format statements. *
*
9999 FORMAT (//3X, 14('====='),//14X,' EXAMPLE ', A /15X,9('--')//
+ 3X, A / 3X, A //2X,' n =', I2 )
9998 FORMAT (//2X,' INITIAL GUESS :', 17X, ' STEPSIZES :'//
+ 9( F16.7, 16X, F16.7 / ) )
9997 FORMAT (//5X,' * * * IMPROPER INPUT PARAMETERS * * *'// )
9996 FORMAT (//5X,'* * * THE CHARACTERISTIC', I2, '-POLYHEDRON HAS',
+ ' NOT BEEN COMPLETED * * *'// )
9995 FORMAT (//2X,' DELTA = ',F20.18 //2X,' EPSILO = ',F20.18
+ ////2X,' FINAL APPROXIMATE SOLUTION :', 5X,
+ 'VERIFICATION OF THE SOLUTION :'//9(F24.18, 9X,F24.18/))
9994 FORMAT (//2X,' EXIT PARAMETERS : INF1 = ',I1,1X,', INF2 = '
+ ,I1//2X,' NUMBER OF FUNCTION CALLS : NFCALL =',I4 )
*
* Last statement of the main program. *
*
END
*---------------------------------------------------------------------*
*
REAL FUNCTION FNC( X, IFLAG )
***********************************************************************
* *
* FUNCTION FNC *
* *
* The purpose of this subprogram is to evaluate the IFLAG-th *
* component of the given function. *
* *
* *
* The function statement is : *
* *
* REAL FUNCTION FNC( X, IFLAG ) *
* *
* where *
* *
* X is an array of length N and contains the corresponding com- *
* ponents of the independent variable. *
* *
* IFLAG determines the component of the function to be evalu- *
* ated. *
* *
***********************************************************************
INTEGER IFLAG, NEXT, N, NFCALL, I
REAL X(9), ZERO, ONE, TWO, FOUR, TEN, ONETEN
COMMON / BLK1 / NEXT, N, NFCALL
PARAMETER ( ZERO = 0.0, ONE = 1.0, TWO = 2.0,
+ FOUR = 4.0, TEN = 10.0, ONETEN = 0.1 )
NFCALL = NFCALL + 1
GO TO ( 100, 200, 300, 400, 500, 600 ), NEXT
100 GO TO ( 1, 2 ), IFLAG
1 FNC = X(1)**2 - FOUR * X(2)
RETURN
2 FNC = X(2)**2 - TWO * X(1) + FOUR * X(2)
RETURN
200 GO TO ( 3, 4 ), IFLAG
3 FNC = ONE - X(1)
RETURN
4 FNC = TEN * ( X(2) - X(1) ** 2 )
RETURN
300 IF ( X(1) .EQ. ZERO ) THEN
IF ( X(2) .EQ. ZERO ) THEN
FNC = ZERO
RETURN
ENDIF
ENDIF
GO TO ( 5, 6 ), IFLAG
5 FNC = ( X(1)**3 - X(2)**3 ) / ( X(1)**2 + X(2)**2 )
RETURN
6 FNC = ( X(1)**3 + X(2)**3 ) / ( X(1)**2 + X(2)**2 )
RETURN
400 FNC = X(IFLAG)
RETURN
500 I = IFLAG + 1
IF ( IFLAG .EQ. N ) I = 1
FNC = ( X(IFLAG) - ONETEN ) ** 2 + ( X(I) - ONETEN )
RETURN
600 I = IFLAG + 1
IF ( IFLAG .EQ. N ) I = 1
FNC = X(IFLAG) ** 2 - X(I)
RETURN
*
* Last statement of the function FNC. *
*
END
*---------------------------------------------------------------------*
*
SUBROUTINE INTSUB( FNC, N, XO, H, DELTA, EPSILO, ICON, INF1,
+ AS, VAS, INF2, WA, LWA )
INTEGER N, ICON, INF1, INF2, LWA
REAL XO(N), H(N), DELTA, EPSILO
REAL AS(N), VAS(N), WA(LWA)
EXTERNAL FNC
***********************************************************************
* *
* SUBROUTINE INTSUB *
* *
* This is an interface subroutine between the main program and *
* the subroutines CHAPOL and GENBIS, through which a single *
* workspace array WA is passed. INTSUB evokes CHAPOL and GENBIS *
* in such a way that the addresses of the corresponding arrays *
* within the subroutines are computed in WA. In this way, their *
* arrays can appear with variable dimensions. *
* *
* *
* The subroutine statement is : *
* *
* SUBROUTINE INTSUB( FNC, N, XO, H, DELTA, EPSILO, ICON, INF1, *
* + AS, VAS, INF2, WA, LWA ) *
* *
* where *
* *
* FNC is the name of the user-supplied function which evaluates *
* components of the given function. FNC should be declared in *
* an external statement in the user - calling program, and *
* should be written as follows : *
* *
* REAL FUNCTION FNC( X, IFLAG ) *
* INTEGER IFLAG *
* REAL X(N) *
* ------------------------------------------------------ *
* Calculate the IFLAG-th component of the function at X. *
* ------------------------------------------------------ *
* RETURN *
* END *
* *
* N is a positive integer input variable which defines the num- *
* ber of equations and variables. *
* *
* XO is an input array of length N which defines the initial *
* guess of the solution. *
* *
* H is an input array of length N, with positive entries, which *
* determines the stepsizes in each coordinate direction. *
* *
* DELTA is a positive input variable which determines the accu- *
* racy of the computation of the roots, of the components of *
* the given function, which are located on the edges of the *
* Initial N-Polyhedron and are used for the construction of a *
* Characteristic N-Polyhedron. If DELTA is less than the *
* machine precision EPSMCH, DELTA takes the value 0.0625 . *
* The value of EPSMCH is computed within INTSUB. *
* *
* EPSILO is a nonnegative input variable. Termination occurs *
* when the algorithm estimates that the Infinity Norm, of the *
* function values at an approximate solution, is at most *
* EPSILO. If EPSILO is less than the machine precision EPSMCH,*
* EPSILO becomes equal to EPSMCH. *
* *
* ICON is an integer input variable. If ICON is equal to 1, *
* then GENBIS is evoked even if a Characteristic N-Polyhedron *
* has not been constructed. *
* *
* INF1 is an integer output variable set as follows : *
* *
* INF1 = 0 Improper input parameters. *
* CHAPOL and GENBIS have not been evoked. *
* *
* INF1 = 1 A Characteristic N-Polyhedron has been con- *
* structed. *
* *
* INF1 = 2 A Characteristic N-Polyhedron has not been con- *
* structed. *
* *
* INF1 = 3 More than two vertices of the constructed Char- *
* acteristic N-Polyhedron are located on the same *
* edge of the Initial N-Polyhedron. *
* *
* INF1 = 4 An approximate solution according to the preci- *
* sion EPSILO has been found during the construc- *
* tion of the Characteristic N-Polyhedron. *
* GENBIS has not been evoked. *
* *
* AS is an output array of length N which determines the final *
* approximate solution. *
* *
* VAS is an output array of length N which specifies the func- *
* tion values at AS. *
* *
* INF2 is an integer output variable set as follows : *
* *
* INF2 = 0 GENBIS has not been evoked. *
* *
* INF2 = 1 The solution has been found within the required *
* accuracy of EPSILO. *
* *
* INF2 = 2 Iterations have reached their upper limit. *
* The answer may not be accurate. *
* *
* INF2 = 3 The length of the longest diagonal has been *
* found to be less than 2.0 * N * EPSILO. *
* *
* WA is a workspace array of length LWA. *
* *
* LWA is a positive integer input variable not less than the *
* value of ( 2*N + (6*N+1) * 2**N ). *
* *
* *
* Subprograms required : *
* *
* USER-Supplied ..... FNC *
* *
* CHABIS-Supplied ... CHAPOL, GENBIS *
* *
* *
* CHABIS. Version of June 1988 *
* Michael N. Vrahatis *
* *
***********************************************************************
INTEGER NV, J, NDG, NN, LD1, LD2, LD3, LD4, LD5, LD6, LD7, LD8
INTEGER NS, LD9, LD10, LD11
REAL ZERO, ONE, TWO, DEF, EPSMCH, TOL
PARAMETER ( ZERO = 0.0, ONE = 1.0,
+ TWO = 2.0, DEF = 0.0625 )
*
* Check the input parameters for errors. *
*
INF1 = 0
INF2 = 0
NV = 2 ** N
IF ( N .LE. 1 .OR. LWA .LT. (2*N+(6*N+1)*NV) ) GO TO 30
DO 10 J = 1, N
IF ( H(J) .LE. ZERO ) GO TO 30
10 CONTINUE
*
* Compute the machine precision. *
*
EPSMCH = ONE
20 EPSMCH = EPSMCH / TWO
TOL = ONE + EPSMCH
IF ( TOL .GT. ONE ) GO TO 20
EPSMCH = TWO * EPSMCH
*
* Determine if the input values of EPSILO and DELTA are less than *
* the machine precision. *
*
IF ( EPSILO .LT. EPSMCH ) EPSILO = EPSMCH
IF ( DELTA .LT. EPSMCH ) DELTA = DEF
*
* Call CHAPOL. *
*
NDG = NV / 2
NN = N * NV
LD1 = 1 + NN
LD2 = LD1 + NN
LD3 = LD2 + NN
LD4 = LD3 + NN
LD5 = LD4 + NN
LD6 = LD5 + NV
LD7 = LD6 + NN
LD8 = LD7 + N
CALL CHAPOL( FNC, N, NV, NDG, XO, H, EPSMCH, DELTA, EPSILO,
+ ICON, WA(1), WA(LD1), WA(LD3), WA(LD4), WA(LD2),
+ WA(LD5),INF1, AS, VAS, WA(LD6), WA(LD7), WA(LD8))
IF ( INF1 .EQ. 4 ) GO TO 30
IF ( INF1 .EQ. 2 .AND. ICON .NE. 1 ) GO TO 30
*
* Call GENBIS. *
*
NS = N * NDG
LD9 = LD6 + NS
LD10 = LD9 + NDG
LD11 = LD10 + N
CALL GENBIS( FNC, N, NV, NS, NDG, WA(1), WA(LD1), WA(LD2), EPSMCH,
+ EPSILO, WA(LD3), WA(LD6), WA(LD9), AS, VAS, INF2,
+ WA(LD5), WA(LD10), WA(LD4), WA(LD11), XO, H )
30 RETURN
*
* Last statement of the interface subroutine INTSUB. *
*
END
*---------------------------------------------------------------------*
*
SUBROUTINE CHAPOL( FNC, N, NV, NDG, XO, H, EPSMCH, DELTA, EPSILO,
+ ICON, B, C, PO, SPO, CP,
+ A, INF1, X, FX, WM, WA1, WA2 )
INTEGER N, NV, NDG, ICON, INF1
REAL XO(N), H(N), EPSMCH, DELTA, EPSILO
REAL B(NV,N), C(NV,N), PO(NV,N), SPO(NV,N), CP(NV,N)
REAL A(NV), X(N), FX(N), WM(NV,N), WA1(N), WA2(N)
***********************************************************************
* *
* SUBROUTINE CHAPOL *
* *
* The purpose of this subroutine is to construct a Character- *
* istic N-Polyhedron for the localization of at least one solu- *
* tion of a system of N non-linear equations in N variables. *
* *
* *
* The subroutine statement is : *
* *
* SUBROUTINE CHAPOL( FNC, N, NV, NDG, XO, H, EPSMCH, DELTA,EPSILO,*
* + ICON, B, C, PO, SPO, CP, *
* + A, INF1, X, FX, WM, WA1, WA2 ) *
* *
* where *
* *
* FNC is the name of the user-supplied function which evaluates *
* components of the given function. FNC should be declared in *
* an external statement in the user - calling program, and *
* should be written as follows : *
* *
* REAL FUNCTION FNC( X, IFLAG ) *
* INTEGER IFLAG *
* REAL X(N) *
* ------------------------------------------------------ *
* Calculate the IFLAG-th component of the function at X. *
* ------------------------------------------------------ *
* RETURN *
* END *
* *
* N is a positive integer input variable which defines the num- *
* ber of equations and variables. *
* *
* NV is a positive integer input variable equal to 2**N. NV *
* specifies the number of vertices of the Initial and Char- *
* acteristic N-Polyhedra. *
* *
* NDG is a positive integer input variable equal to NV/2. NDG *
* specifies the number of diagonals of the Characteristic N- *
* Polyhedron. *
* *
* XO is an input array of length N which defines the initial *
* guess of the solution. *
* *
* H is an input array of length N, with positive entries, which *
* determines the stepsizes in each coordinate direction. *
* *
* EPSMCH is a positive input variable which determines the ma- *
* chine precision. *
* *
* DELTA is a positive input variable which determines the accu- *
* racy of the computation of the roots, of the components of *
* the given function, which are located on the edges of the *
* Initial N-Polyhedron and are used for the construction of a *
* Characteristic N-Polyhedron. If DELTA is less than the *
* machine precision EPSMCH, DELTA takes the value 0.0625 . *
* *
* EPSILO is a nonnegative input variable. Termination occurs *
* when the algorithm estimates that the Infinity Norm, of the *
* function values at an approximate solution, is at most *
* EPSILO. If EPSILO is less than the machine precision EPSMCH,*
* EPSILO becomes equal to EPSMCH. *
* *
* ICON is an integer input variable. If ICON is equal to 1 and *
* the Characteristic N-Polyhedron construction fails, then *
* CHAPOL will complete the unconstructed Characteristic N- *
* Polyhedron vertices, using the Initial N-Polyhedron verti- *
* ces. Note that this action is taken when the Characteristic *
* N-Polyhedron is incomplete and the subroutine GENBIS is *
* utilized. *
* *
* B is an output NV by N matrix which defines the N-binary *
* matrix. *
* *
* C is an output NV by N matrix which defines the N-complete *
* matrix. *
* *
* PO is an output NV by N matrix,with entries in each row which *
* are the corresponding components of the vertices of the *
* Initial N-Polyhedron. *
* *
* SPO is an output NV by N matrix, with entries in each row *
* which are the corresponding components of the vectors of *
* signs of the function relative to the vertices of PO. *
* *
* CP is an output NV by N matrix,with entries in each row which *
* are the corresponding components of the vertices of the *
* Characteristic N-Polyhedron. *
* *
* A is an output array of length NV and determines which verti- *
* ces of the Characteristic N-Polyhedron have not been con- *
* structed. When all its entries are equal to zero, the Char- *
* acteristic N-Polyhedron has been constructed. *
* *
* INF1 is an integer output variable set as follows : *
* *
* INF1 = 0 Improper input parameters. *
* CHAPOL and GENBIS have not been evoked. *
* *
* INF1 = 1 A Characteristic N-Polyhedron has been con- *
* structed. *
* *
* INF1 = 2 A Characteristic N-Polyhedron has not been con- *
* structed. *
* *
* INF1 = 3 More than two vertices of the constructed Char- *
* acteristic N-Polyhedron are located on the same *
* edge of the Initial N-Polyhedron. *
* *
* INF1 = 4 An approximate solution, according to the preci- *
* sion EPSILO, has been found during the construc- *
* tion of the Characteristic N-Polyhedron. *
* GENBIS has not been evoked. *
* *
* X is an output array of length N which determines a point *
* which will be a vertex of the Characteristic N-Polyhedron. *
* In the case of INF1 = 4, X contains the coordinates of the *
* approximate solution. *
* *
* FX is an output array of length N which determines the func- *
* tion values at X . *
* *
* WM is an NV by N work matrix. *
* *
* WA1, WA2 are work arrays of length N. *
* *
* *
* Subprograms required : *
* *
* USER-Supplied ...... FNC *
* *
* BLAS-Supplied ...... SCOPY *
* *
* CHABIS-Supplied .... SGN, LSINNC, ISCOMC, SAVERC, SSCALC, *
* SMNLGC *
* *
* FORTRAN-Supplied ... FLOAT, LOG, INT, ABS, SIGN *
* *
* *
* CHABIS. Version of June 1988 *
* Michael N. Vrahatis *
* *
***********************************************************************
INTRINSIC FLOAT, LOG, INT, ABS, SIGN
INTEGER J, K, L, I, I1, J1, J2, J3, J4, J5, I2, ITR, I3, NNV, I4
REAL ZERO, ONE, TWO, DSTAR, T, EL, ER, RLEN
REAL RITR, RO, SRO, S, RN, RT, FE, SRN, SC1, SC2
LOGICAL TEST, LSINNC
PARAMETER ( ZERO = 0.0, ONE = 1.0, TWO = 2.0 )
*
* Set DSTAR, a small positive number, which determines a toler- *
* ance of the roots of the components of the given function, *
* which are located on the edges of the Initial N-Polyhedron. *
* DSTAR must be greater than, or equal to, DELTA. *
*
DSTAR = DELTA + TWO * EPSMCH
*
* Construct the N-Complete Matrix and the matrix PO. Also, set *
* all the entries of CP equal to the corresponding entries of PO. *
*
DO 20 J = 1, N
K = 2 ** ( N - J )
L = 2 * K
DO 10 I = 1, NV
I1 = I - 1
T = FLOAT( I1 / K - 2 * ( I1 / L ) )
B(I,J) = T
C(I,J) = TWO * T - ONE
PO(I,J) = XO(J) + T * H(J)
CP(I,J) = PO(I,J)
10 CONTINUE
20 CONTINUE
*
* Construct the indexing array, A. *
*
DO 30 I = 1, NV
A(I) = FLOAT( I )
30 CONTINUE
*
* Construct the matrix of signs of the given function, relative *
* to vertices of the Initial N-Polyhedron. *
*
DO 50 I = 1, NV
CALL SCOPY( N, PO(I,1), NV, X, 1 )
DO 40 J = 1, N
FX(J) = FNC( X, J )
SPO(I,J) = SGN( FX(J) )
40 CONTINUE
TEST = LSINNC( N, FX, EPSILO )
IF ( TEST ) GO TO 270
*
* Start the construction of a Characteristic N-Polyhedron. *
*
CALL SCOPY( N, SPO(I,1), NV, WA1, 1 )
K = ISCOMC( N, WA1, C, NV )
IF ( K .GT. 0 ) THEN
IF ( A(K) .GT. ZERO ) THEN
CALL SCOPY( N, X, 1, CP(K,1), NV )
A(K) = ZERO
ENDIF
ENDIF
50 CONTINUE
*
* Determine if a Characteristic N-Polyhedron has been completed. *
*
TEST = LSINNC( NV, A, EPSMCH )
IF ( TEST ) GO TO 250
*
* Compute the roots of the components of the function, which are *
* located on the edges of the Initial N-Polyhedron, and use them *
* to construct a Characteristic N-Polyhedron. *
*
DO 130 J1 = 1, N
J2 = 2 ** J1 - 2
J3 = 2 ** ( N - J1 )
DO 120 J4 = 0, J2, 2
DO 110 J5 = 1, J3
I1 = J4 * J3 + J5
I2 = I1 + J3
CALL SCOPY( N, PO(I1,1), NV, X, 1 )
EL = PO(I1,J1)
ER = PO(I2,J1)
RLEN = ER - EL
*
* Compute the number of iterations which are required *
* to obtain a root to the predetermined accuracy, DELTA. *
*
RITR = LOG( RLEN / DELTA ) / LOG( TWO )
K = INT( RITR )
ITR = K + 1
IF ( RITR .EQ. FLOAT( K ) ) ITR = ITR - 1
*
DO 100 J = 1, N
IF ( SPO(I1,J) * SPO(I2,J) .GT. ZERO ) GO TO 100
RO = EL
SRO = SPO(I1,J)
S = RLEN / TWO
RN = RO + S
RT = ER
*
* Begin the iterations. *
*
DO 60 I3 = 1, ITR
X(J1) = RN
FE = FNC( X, J )
SRN = SGN( FE )
RO = RN
S = S / TWO
*
* Compute the new approximate root. *
*
RN = RO + SRO * SRN * S
*
* Specify another stopping test using the differ- *
* ence between two successive iterations. *
*
IF ( ABS( RN - RO ) .LE. DELTA ) GO TO 70
60 CONTINUE
70 RT = RN
IF ( ( ER - RT ) .LE. DELTA ) GO TO 100
*
* Construct a Characteristic N-Polyhedron using the *
* computed roots. *
*
DO 90 K = -1, 1, 2
X(J1) = RT - FLOAT( K ) * DSTAR
DO 80 L = 1, N
FX(L) = FNC( X, L )
WA1(L) = SGN( FX(L) )
80 CONTINUE
TEST = LSINNC( N, FX, EPSILO )
IF ( TEST ) GO TO 270
I = ISCOMC( N, WA1, C, NV )
IF ( I .NE. 0 ) THEN
CALL SCOPY( N, X, 1, CP(I,1), NV )
A(I) = ZERO
*
* Determine if a Characteristic N-Polyhedron *
* has been completed. *
*
TEST = LSINNC( NV, A, EPSMCH )
IF ( TEST ) GO TO 150
ENDIF
90 CONTINUE
100 CONTINUE
110 CONTINUE
120 CONTINUE
130 CONTINUE
*
* Complete the Characteristic N-Polyhedron vertices, using the *
* Initial N-Polyhedron vertices. Note that this action is taken *
* when the Characteristic N-Polyhedron is incomplete, and the *
* the subroutine GENBIS is utilized. *
*
IF ( ICON .NE. 1 ) GO TO 260
CALL SAVERC( N, PO, NV, X )
DO 140 I = 1, NV
L = 1 - I + NV
IF ( A(I) .GT. ZERO ) THEN
IF ( A(L) .EQ. ZERO ) THEN
SC1 = TWO
SC2 = - ONE
CALL SSCALC( N, SC1, X, 1, 1, SC2, CP, L, NV,
+ CP, I, NV )
ENDIF
ENDIF
140 CONTINUE
GO TO 260
*
* Reconstruct the Characteristic N-Polyhedron in such a way that *
* its vertices are the vertices of a scaled translation of the *
* unit N-cube. *
*
150 INF1 = 1
NNV = N * NV
CALL SCOPY( NNV, CP(1,1), 1, WM(1,1), 1 )
CALL SMNLGC( N, CP, NV, WA1, WA2 )
DO 190 I1 = 1, NV
DO 160 J = 1, N
X(J) = WA1(J) + B(I1,J) * WA2(J)
160 CONTINUE
K = ISCOMC( N, X, WM, NV )
IF ( K .NE. 0 ) GO TO 190
K = ISCOMC( N, X, PO, NV )
IF ( K .NE. 0 ) GO TO 190
DO 170 J = 1, N
FX(J) = FNC( X, J )
170 CONTINUE
TEST = LSINNC( N, FX, EPSILO )
IF ( TEST ) GO TO 270
*
* Change vertices. *
*
DO 180 J = 1, N
FX(J) = SIGN( ONE, FX(J) )
180 CONTINUE
I = ISCOMC( N, FX, C, NV )
IF ( I .EQ. 0 ) GO TO 190
CALL SCOPY( N, X, 1, CP(I,1), NV )
190 CONTINUE
*
* Reconstruct the Characteristic N-Polyhedron in such a way that *
* not more than two of its vertices are located on the same edge *
* of the Initial N-Polyhedron. *
*
DO 240 I1 = 1, NDG
I2 = 1 - I1 + NV
DO 200 K = 1, N
IF ( ABS( CP(I1,K) - CP(I2,K) ) .LT. EPSMCH ) GO TO 210
200 CONTINUE
GO TO 240
210 INF1 = 3
DO 230 I3 = 1, 2
I4 = I1
IF ( I3 .EQ. 2 ) I4 = I2
CALL SCOPY( N, CP(I4,1), NV, X, 1 )
IF ( ABS( X(K) - XO(K) ) .LT. EPSMCH ) THEN
X(K) = XO(K) + H(K)
ELSE
X(K) = X(K) - H(K)
ENDIF
DO 220 J = 1, N
FX(J) = FNC( X, J )
WA1(J) = SIGN( ONE, FX(J) )
220 CONTINUE
TEST = LSINNC( N, FX, EPSILO )
IF ( TEST ) GO TO 270
*
* Change vertices. *
*
I = ISCOMC( N, WA1, C, NV )
IF ( I .EQ. 0 ) GO TO 240
IF ( I4 .EQ. I ) THEN
INF1 = 1
CALL SCOPY( N, X, 1, CP(I,1), NV )
GO TO 240
ENDIF
230 CONTINUE
240 CONTINUE
*
RETURN
250 INF1 = 1
RETURN
260 INF1 = 2
RETURN
270 INF1 = 4
RETURN
*
* Last statement of the subroutine CHAPOL. *
*
END
*---------------------------------------------------------------------*
*
SUBROUTINE GENBIS( FNC, N, NV, NS, NDG, B, C, CP, EPSMCH,
+ EPSILO, CPR, DM, DG, BP, FVB, INF2,
+ A, ARP, WM, WA1, WA2, WA3 )
INTEGER N, NV, NS, NDG, INF2
REAL B(NV,N), C(NV,N), CP(NV,N), EPSMCH, EPSILO
REAL CPR(NV,N), DM(NS), DG(NDG), BP(N), FVB(N)
REAL A(NV), ARP(N), WM(NV,N), WA1(N), WA2(N), WA3(N)
***********************************************************************
* *
* SUBROUTINE GENBIS *
* *
* The purpose of this subroutine is to bisect the Character- *
* istic N-Polyhedron in order to find an approximate solution, *
* according to the predetermined accuracy of EPSILO. *
* *
* *
* The subroutine statement is : *
* *
* SUBROUTINE GENBIS( FNC, N, NV, NS, NDG, B, C, CP, EPSMCH, *
* + EPSILO, CPR, DM, DG, BP, FVB, INF2, *
* + A, ARP, WM, WA1, WA2, WA3 ) *
* *
* where *
* *
* FNC is the name of the user-supplied function which evaluates *
* components of the given function. FNC should be declared in *
* an external statement in the user - calling program, and *
* should be written as follows : *
* *
* REAL FUNCTION FNC( X, IFLAG ) *
* INTEGER IFLAG *
* REAL X(N) *
* ------------------------------------------------------ *
* Calculate the IFLAG-th component of the function at X. *
* ------------------------------------------------------ *
* RETURN *
* END *
* *
* N is a positive integer input variable which defines the num- *
* ber of equations and variables. *
* *
* NV is a positive integer input variable equal to 2**N. NV *
* specifies the number of vertices of the Initial and Char- *
* acteristic N-Polyhedra. *
* *
* NS is a positive integer input variable equal to N*NV/2. *
* NS specifies the number of proper 1-simplexes of the Char- *
* acteristic N-Polyhedron. *
* *
* NDG is a positive integer input variable equal to NV/2. NDG *
* specifies the number of diagonals of the Characteristic N- *
* Polyhedron. *
* *
* B is an input NV by N matrix which defines the N-binary *
* matrix. *
* *
* C is an input NV by N matrix which defines the N-complete *
* matrix. *
* *
* CP is an input NV by N matrix, with entries in each row which *
* are the corresponding components of the vertices of the *
* Characteristic N-Polyhedron. In the output, CP remains un- *
* changed. *
* *
* EPSMCH is a positive input variable which determines the ma- *
* chine precision. *
* *
* EPSILO is a nonnegative input variable. Termination occurs *
* when the algorithm estimates that the Infinity Norm, of the *
* function values at an approximate solution, is at most *
* EPSILO. If EPSILO is less than the machine precision EPSMCH,*
* EPSILO becomes equal to EPSMCH. *
* *
* CPR is an output NV by N matrix, with entries in each row *
* which are the corresponding elements of CP, and determines *
* the Characteristic N-Polyhedron used for the iterative re- *
* finement. *
* *
* DM is an output array of length NS which contains the lengths *
* of all the proper 1-simplexes of the Characteristic N-Poly- *
* hedron. *
* *
* DG is an output array of length NDG which contains the lengths*
* of all the diagonals of the Characteristic N-Polyhedron. *
* *
* BP is an output array of length N which determines an approxi-*
* mate solution. *
* *
* FVB is an output array of length N which determines the func- *
* tion values at BP. *
* *
* INF2 is an integer output variable set as follows : *
* *
* INF2 = 0 GENBIS has not been evoked. *
* *
* INF2 = 1 The solution has been found within the required *
* accuracy of EPSILO. *
* *
* INF2 = 2 Iterations have reached their upper limit. *
* The answer may not be accurate. *
* *
* INF2 = 3 The length of the longest diagonal has been *
* found to be less than 2.0 * N * EPSILO. *
* *
* A is an output array of length NV which determines which *
* vertices of the Characteristic N-Polyhedron are substituted *
* during the bisection of its proper 1-simplexes. *
* *
* ARP is an output array of length N which is used during the *
* relaxation procedure. *
* *
* WM is an NV by N work matrix. *
* *
* WA1, WA2, WA3 are work arrays of length N. *
* *
* *
* Subprograms required : *
* *
* USER-Supplied ...... FNC *
* *
* BLAS-Supplied ...... SCOPY, SNRM2, ISAMAX *
* *
* CHABIS-Supplied .... SSCALC, LSINNC, ISCOMC, SMNLGC *
* *
* FORTRAN-Supplied ... FLOAT, LOG, INT, SIGN *
* *
* *
* CHABIS. Version of June 1988 *
* Michael N. Vrahatis *
* *
***********************************************************************
INTRINSIC FLOAT, LOG, INT, SIGN
INTEGER NIRP, NNV, I3, J1, J2, J3, J4, J5, I1, I2, I, K, ITR
INTEGER IRP1, J, IRP2, M
REAL ZERO, ONE, TWO, HALF, ZETA, SC1, SC2, GDM, RITR
REAL GDG
LOGICAL TEST, LSINNC
PARAMETER ( ZERO = 0.0, ONE = 1.0,
+ TWO = 2.0, HALF = 0.5 )
*
* Set ZETA, a small positive real variable, which determines the *
* accuracy of a stopping criterion during the bisection process. *
* The algorithm terminates if it is determined that the longest *
* diagonal of a refined Characteristic N-Polyhehedron has a *
* length less than ZETA. *
*
ZETA = TWO * FLOAT( N ) * EPSILO
*
* Set the maximum number of iterations of the relaxation process. *
*
NIRP = 2
*
* Set all the entries of the matrix CPR equal to the correspond- *
* ing entries of CP. *
*
NNV = N * NV
CALL SCOPY( NNV, CP(1,1), 1, CPR(1,1), 1 )
*
* Compute the lengths of all the proper 1-simplexes of the Char- *
* acteristic N-Polyhedron and store them in the array, DM. *
*
I3 = 0
DO 30 J1 = 1, N
J2 = 2 ** J1 - 2
J3 = 2 ** ( N - J1 )
DO 20 J4 = 0, J2, 2
DO 10 J5 = 1, J3
I3 = I3 + 1
I1 = J4 * J3 + J5
I2 = I1 + J3
SC1 = ONE
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I1, NV, SC2, CPR, I2, NV,
+ WA1, 1, 1 )
DM(I3) = SNRM2( N, WA1, 1 )
10 CONTINUE
20 CONTINUE
30 CONTINUE
*
* Compute the length of the longest proper 1-simplexes. *
*
I = ISAMAX( NS, DM, 1 )
GDM = DM( I )
*
* Compute the maximum number of iterations. *
*
RITR = LOG( GDM * TWO / (FLOAT(N)*EPSILO) ) / LOG( TWO )
K = INT( RITR )
ITR = K + 1
IF ( RITR .EQ. FLOAT( K ) ) ITR = ITR - 1
*
* Begin the iterations. *
*
IRP1 = 0
DO 230 I3 = 1, ITR
*
* Bisect the diagonals. *
*
IF ( IRP1 .NE. 0 ) GO TO 80
DO 60 I1 = 1, NDG
I2 = 1 - I1 + NV
40 SC1 = HALF
CALL SSCALC( N, SC1, CPR, I1, NV, SC1, CPR, I2, NV,
+ BP, 1, 1 )
DO 50 J = 1, N
FVB(J) = FNC( BP, J )
WA1(J) = SIGN( ONE, FVB(J) )
50 CONTINUE
*
* Stopping test. *
* The solution has been found within the required accuracy. *
*
TEST = LSINNC( N, FVB, EPSILO )
IF ( TEST ) GO TO 270
*
* Change vertices. *
*
I = ISCOMC( N, WA1, C, NV )
IF ( I .EQ. 0 ) GO TO 60
CALL SCOPY( N, BP, 1, CPR(I,1), NV )
IF (I .EQ. I1 .OR. I .EQ. I2) GO TO 40
60 CONTINUE
*
* Compute the length of the longest diagonal, GDG. *
*
DO 70 I1 = 1, NDG
I2 = 1 - I1 + NV
SC1 = ONE
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I1, NV, SC2, CPR, I2, NV,
+ WA1, 1, 1 )
DG(I1) = SNRM2( N, WA1, 1 )
70 CONTINUE
I = ISAMAX( NDG, DG, 1 )
GDG = DG( I )
*
* Determine if GDG is less than ZETA. *
*
IF ( GDG .LT. ZETA ) THEN
INF2 = 3
GO TO 250
ENDIF
*
80 IRP1 = 0
*
* Construct the indexing array, A. *
*
DO 90 I = 1, NV
A(I) = FLOAT( I )
90 CONTINUE
*
* Bisect the proper 1-simplexes. *
*
DO 170 J1 = 1, N
J2 = 2 ** J1 - 2
J3 = 2 ** ( N - J1 )
DO 160 J4 = 0, J2, 2
DO 150 J5 = 1, J3
I1 = J4 * J3 + J5
I2 = I1 + J3
IRP2 = 0
SC1 = HALF
CALL SSCALC( N, SC1, CPR, I1, NV, SC1, CPR, I2, NV,
+ BP, 1, 1 )
GO TO 110
*
* Execute the relaxation procedure. *
*
100 IRP2 = IRP2 + 1
IRP1 = 1
SC1 = TWO
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I, NV, SC2, ARP, 1, 1,
+ BP, 1, 1 )
110 DO 120 J = 1, N
WA3(J) = ZERO
FVB(J) = FNC( BP, J )
IF ( FVB(J) .EQ. ZERO ) WA3(J) = FLOAT( J )
WA1(J) = SIGN( ONE, FVB(J) )
120 CONTINUE
*
* Stopping test. *
* The solution has been found within the required *
* accuracy. *
*
TEST = LSINNC( N, FVB, EPSILO )
IF ( TEST ) GO TO 270
*
* Change vertices. *
*
I = ISCOMC ( N, WA1, C, NV )
IF ( I .EQ. 0 ) GO TO 150
DO 130 K = 1, N
J = INT( WA3(K) )
IF ( WA3(K) .EQ. ZERO ) GO TO 130
M = I - INT( C(I,J) * TWO ** ( N - J ) )
SC1 = ONE
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I, NV,
+ SC2, CPR, M, NV, WA2, 1, 1 )
TEST = LSINNC( N, WA2, EPSMCH )
IF ( TEST ) GO TO 140
130 CONTINUE
CALL SCOPY( N, CPR(I,1), NV, ARP, 1 )
CALL SCOPY( N, BP, 1, CPR(I,1), NV )
A(I) = ZERO
140 IF ( I .NE. I1 .AND. I .NE. I2
+ .AND. IRP2 .LT. NIRP ) GO TO 100
150 CONTINUE
160 CONTINUE
170 CONTINUE
*
* Compute the length of the longest diagonal, GDG. *
*
DO 180 I1 = 1, NDG
I2 = 1 - I1 + NV
SC1 = ONE
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I1, NV, SC2, CPR, I2, NV,
+ WA1, 1, 1 )
DG(I1) = SNRM2( N, WA1, 1 )
180 CONTINUE
I = ISAMAX( NDG, DG, 1 )
GDG = DG( I )
*
* Determine if GDG is less than ZETA. *
*
IF ( GDG .LT. ZETA ) THEN
INF2 = 3
GO TO 250
ENDIF
*
K = 0
DO 190 I = 1, NV
IF ( A(I) .NE. ZERO ) K = K + 1
190 CONTINUE
IF ( K .EQ. 0 .OR. IRP1 .EQ. 0 ) THEN
IRP1 = 0
GO TO 230
ENDIF
*
* Reconstruct the refined Characteristic N-Polyhedron in such *
* a way that its vertices are the vertices of a scaled trans- *
* lation of the unit N-cube. *
*
CALL SCOPY( NNV, CPR(1,1), 1, WM(1,1), 1 )
CALL SMNLGC( N, CPR, NV, WA1, WA2 )
DO 220 I1 = 1, NV
DO 200 J = 1, N
BP(J) = WA1(J) + B(I1,J) * WA2(J)
200 CONTINUE
I = ISCOMC( N, BP, WM, NV )
IF ( I .NE. 0 ) GO TO 220
DO 210 J = 1, N
FVB(J) = FNC( BP, J )
WA3(J) = SIGN( ONE, FVB(J) )
210 CONTINUE
*
* Stopping test. *
* The solution has been found within the required accuracy. *
*
TEST = LSINNC( N, FVB, EPSILO )
IF ( TEST ) GO TO 270
*
* Change vertices. *
*
I = ISCOMC( N, WA3, C, NV )
IF ( I .EQ. 0 ) GO TO 220
CALL SCOPY( N, BP, 1, CPR(I,1), NV )
220 CONTINUE
230 CONTINUE
*
INF2 = 2
*
* Find the longest diagonal of the refined Characteristic N-Poly- *
* hedron and take its midpoint as the final approximate solution. *
*
DO 240 I1 = 1, NDG
I2 = 1 - I1 + NV
SC1 = ONE
SC2 = - ONE
CALL SSCALC( N, SC1, CPR, I1, NV, SC2, CPR, I2, NV, WA1, 1, 1)
DG(I1) = SNRM2( N, WA1, 1 )
240 CONTINUE
I = ISAMAX( NDG, DG, 1 )
250 I1 = 1 - I + NV
SC1 = HALF
CALL SSCALC( N, SC1, CPR, I, NV, SC1, CPR, I1, NV, BP, 1, 1 )
DO 260 J = 1, N
FVB(J) = FNC( BP, J )
260 CONTINUE
*
* Stopping test. *
* The solution has been found within the required accuracy. *
*
TEST = LSINNC( N, FVB, EPSILO )
IF ( TEST ) GO TO 270
*
RETURN
270 INF2 = 1
RETURN
*
* Last statement of the subroutine GENBIS. *
*
END
*---------------------------------------------------------------------*
*---------------------------------------------------------------------*
* *
* B L A S - L I K E U T I L I T I E S *
* *
* CHABIS. Version of June 1988 *
* *
* Michael N. Vrahatis *
* *
*---------------------------------------------------------------------*
*---------------------------------------------------------------------*
*
LOGICAL FUNCTION LSINNC( N, SX, EPSILO )
*
* Determine if the Infinity Norm of a single precision N-vector *
* stored in SX( ) is less than, or equal to, the single precision *
* positive constant EPSILO; if so, set LSINNC equal to " TRUE ";*
* otherwise, set LSINNC equal to " FALSE ". *
*
INTRINSIC ABS
INTEGER N, I
REAL SX(1), EPSILO
*
LSINNC = .FALSE.
DO 10 I = 1, N
IF ( ABS( SX(I) ) .GT. EPSILO ) RETURN
10 CONTINUE
LSINNC = .TRUE.
*
RETURN
END
*---------------------------------------------------------------------*
*
INTEGER FUNCTION ISCOMC( N, SX, SY, INCY )
*
* Find the value of I, where I can be I = 1,2,...,INCY, such that *
* the entries of a single precision N-vector, stored in SX( J ), *
* J = 1,2,...,N, coincide with the corresponding single precision *
* entries SY( I+(J-1)*INCY ) for all J = 1,2,...,N. If there is *
* no such value, then return the value of zero. *
*
INTEGER N, INCY, I, K, J
REAL SX(1), SY(1)
*
DO 20 I = 1, INCY
K = I
DO 10 J = 1, N
IF ( SX(J) .NE. SY(K) ) GO TO 20
K = K + INCY
10 CONTINUE
ISCOMC = I
RETURN
20 CONTINUE
ISCOMC = 0
*
RETURN
END
*---------------------------------------------------------------------*
*
REAL FUNCTION SGN( SV )
*
* Transfer the sign of the single precision variable SV in such a *
* way that if SV is equal to zero, the sign of SV will be equal *
* to zero. *
*
INTRINSIC SIGN
REAL SV, ZERO, ONE
PARAMETER ( ZERO = 0.0, ONE = 1.0 )
*
SGN = SIGN( ONE, SV )
IF ( SV .EQ. ZERO ) SGN = ZERO
*
RETURN
END
*---------------------------------------------------------------------*
*
SUBROUTINE SAVERC( N, SX, INCX, SY )
*
* Find the average value for each of the N vector segments (each *
* of which contains INCX elements) of the single precision N- *
* vector stored in SX( ). Store these values in the single preci- *
* sion N-vector SY( ). *
*
INTRINSIC MOD, FLOAT
INTEGER N, INCX, M, K, J, K1, K2, I
REAL SX(1), SY(1), ZERO, SUM
PARAMETER ( ZERO = 0.0 )
*
M = MOD( INCX, 4 )
K = - INCX
DO 50 J = 1, N
K = K + INCX
SUM = ZERO
IF ( M .EQ. 0 ) GO TO 20
K1 = K + 1
K2 = K + M
DO 10 I = K1, K2
SUM = SUM + SX( I )
10 CONTINUE
IF ( INCX .LT. 4 ) GO TO 40
20 K1 = K + M + 1
K2 = K + INCX
DO 30 I = K1, K2, 4
SUM = SUM + SX( I ) + SX( I+1 ) + SX( I+2 ) + SX( I+3 )
30 CONTINUE
40 SY(J) = SUM / FLOAT( INCX )
50 CONTINUE
*
RETURN
END
*---------------------------------------------------------------------*
*
SUBROUTINE SSCALC( N, V1, SX, LX, INCX, V2, SY, LY, INCY,
+ SZ, LZ, INCZ )
*
* Multiply the single precision variable V1 by the single preci- *
* sion N-vector stored in SX( ), with leading dimension LX and *
* storage increment INCX ; to the product, add the single preci- *
* sion variable V2 multiplied by the single precision N-vector *
* stored in SY( ),with leading dimension LY and storage increment *
* INCY. Store the result in the single precision N-vector SZ( ), *
* with leading dimension LZ and storage increment INCZ. *
*
INTRINSIC MOD
INTEGER N, LX, INCX, LY, INCY, LZ, INCZ, K, IX, IY, IZ, I
REAL V1, SX(1), V2, SY(1), SZ(1)
*
K = MOD( N, 4 )
IX = LX
IY = LY
IZ = LZ
IF ( K .EQ. 0 ) GO TO 20
DO 10 I = 1, K
SZ( IZ ) = V1 * SX( IX ) + V2 * SY( IY )
IX = IX + INCX
IY = IY + INCY
IZ = IZ + INCZ
10 CONTINUE
IF ( N .LT. 4 ) RETURN
20 DO 30 I = K, N-1, 4
SZ( IZ ) = V1 * SX( IX ) + V2 * SY( IY )
IX = IX + INCX
IY = IY + INCY
IZ = IZ + INCZ
SZ( IZ ) = V1 * SX( IX ) + V2 * SY( IY )
IX = IX + INCX
IY = IY + INCY
IZ = IZ + INCZ
SZ( IZ ) = V1 * SX( IX ) + V2 * SY( IY )
IX = IX + INCX
IY = IY + INCY
IZ = IZ + INCZ
SZ( IZ ) = V1 * SX( IX ) + V2 * SY( IY )
IX = IX + INCX
IY = IY + INCY
IZ = IZ + INCZ
30 CONTINUE
*
RETURN
END
*---------------------------------------------------------------------*
*
SUBROUTINE SMNLGC( N, SX, INCX, SMN, SLG )
* *
* Find the maximum value, and the difference between the maximum *
* and minimum values, within each of the N vector segments (each *
* of which contains INCX elements) of the single precision N- *
* vector stored in SX( ). Store these values in the single pre- *
* cision N-vectors SMN( ) and SLG( ), respectively. *
*
INTRINSIC MIN, MAX
INTEGER N, INCX, K, J, I
REAL SX(1), SMN(1), SLG(1), SMIN, SMAX, X
*
K = - INCX
DO 20 J = 1, N
K = K + INCX
SMIN = SX( K + 1 )
SMAX = SMIN
DO 10 I = K+2, K+INCX
X = SX( I )
SMIN = MIN( SMIN, X )
SMAX = MAX( SMAX, X )
10 CONTINUE
SMN(J) = SMIN
SLG(J) = SMAX - SMIN
20 CONTINUE
*
RETURN
END