C C BASIS C C **** VERSION OF AUGUST 15, 1983 C C BASIS IS A COLLECTION OF FIVE SUBROUTINES THAT DEAL C WITH THE BASIS MATRIX BAS C C BAS IS A SQUARE MATRIX OF ORDER NP1 WHOSE COLUMNS C ARE OF THE FORM OF EITHER C C (F(X(.,I)),SCALE) WHERE X(.,I) CORRESPONDS TO THE ITH C VERTEX OF THE CURRENT PIECE (SEE SUBDIV) C C (ADERIV(.,K)*DELTA(K),0), A SCALED COLUMN OF THE APPROXIMATE C DERIVATIVE MATRIX ADERIV(N,N), OR C C (ADERIV*(XSPECIAL - XORIG),SCALE), WHERE XSPECIAL IS CLOSE TO C THE CURRENT PIECE AND XORIG C IS THE STARTING POINT OF THE C MAJOR CYCLE. C C THE LAST KIND OF COLUMN IS CALLED SPECIAL. THERE IS AT MOST ONE C OF THESE COLUMNS IN BAS, AND IF PRESENT IT IS IN THE NP1'ST C POSITION. C C SCALE IS CHOSEN TO MAKE THE ENTRIES OF BAS ROUGHLY EQUAL FOR C NUMERICAL REASONS. C C LPERM IS A PERMUTATION OF (1,2,...,NP1) WITH NP1 FIXED C SUPPOSE P IS A PERMUTATION MATRIX WITH P(L,J) = 1 IFF C LPERM(J) = L. THEN THE BASIS MATRIX IS USED VIA THE FACTORIZATION C C BAS * P = Q * R C C WHERE Q IS AN ORTHOGONAL AND R AN UPPER TRIANGULAR MATRIX OF C ORDER NP1. Q IS NOT STORED. C C THE POINTERS IPOINT AND KPOINT DESIGNATE THE ORIGIN OF THE C COLUMNS OF R. C C IF IPOINT(J) = I.GT.0, THE JTH COLUMN OF R AND THE LPERM(J)'TH C COLUMN OF BAS CORRESPOND TO THE ITH VERTEX C OF THE CURRENT PIECE, AND KPOINT(J) = 0 C C IF KPOINT(J) = K.GT.0, THESE COLUMNS CORRESPOND TO THE SCALED C KTH COLUMN OF ADERIV, AND IPOINT(J) = 0 C C KPOINT(NP1) IS ALWAYS ZERO AND IPOINT(NP1) IS ONLY NONZERO C WHEN THE SPECIAL COLUMN LEAVES THE BASIS AND A MAJOR CYCLE C IS ENDED. C C THE BASIS MATRIX IS ALWAYS NONSINGULAR AND WEIGHT IS AN NP1 - C VECTOR SOLVING C C BAS * P * WEIGHT = (0,0,...,0,SCALE) C C WEIGHT ALSO SATISFIES THE FOLLOWING INEQUALITIES C C WEIGHT(J).GE.0 IF (IPOINT(J).GT.0).OR.(J.EQ.NP1) C C WEIGHT(NP1).GE.WEIGHT(J).GE.( - WEIGHT(NP1)) OTHERWISE C C THE SUBROUTINES IN BASIS ARE FSTBAS,NEWBAS,PRTBAS,LABEL C AND MINRAT C C FSTBAS INITIALIZES THE MATRICES BAS AND R, THE POINTERS LPERM, C IPOINT AND KPOINT, AND THE VECTOR WEIGHT C C NEWBAS UPDATES THESE WHEN A NEW COLUMN IS BROUGHT INTO THE BASIS C IT ALSO DETERMINES THE LEAVING COLUMN C C PRTBAS PRINTS ALL THESE ARGUMENTS WHEN CALLED C C LABEL DETERMINES THE COLUMN THAT IS TO ENTER THE BASIS C C MINRAT IS CALLED BY NEWBAS TO DETERMINE THE LEAVING COLUMN C C FSTBAS AND NEWBAS REQUIRE THE LINALG SUBROUTINES QRGEN, QRHESS C AND SOLVE TO MANIPULATE THE FACTORIZATION C C ****** C SUBROUTINE FSTBAS(N,NP1,NEFFP1,PRTLEV,BAS,R,FR,WEIGHT,ADERIV, * DELTA,XORIG,SCALE,LPERM,FLPERM,IPOINT,KPOINT, * TYPJAC,JACFAC,STPDWN,SIGMA,WA,FLAG,VDLT) INTEGER N,NP1,NEFFP1,PRTLEV,LPERM(NP1),FLPERM(NP1) INTEGER IPOINT(NP1),KPOINT(NP1),TYPJAC,JACFAC,FLAG,VDLT DOUBLE PRECISION BAS(NP1,NP1),R(NP1,NP1),FR(*),WEIGHT(NP1) DOUBLE PRECISION ADERIV(N,N),DELTA(N),SCALE,SIGMA(NP1),WA(NP1) DOUBLE PRECISION XORIG(N) DOUBLE PRECISION STPDWN C C FSTBAS SETS UP THE INITIAL BASIS MATRIX BAS, THE PERMUTATION C LPERM, THE UPPER TRIANGULAR FACTOR R, THE POINTERS IPOINT AND C KPOINT AND THE VECTOR WEIGHT C C TO AVOID REPETITION OF THE FACTORIZATION, FR AND FLPERM ARE C SET TO THE INITIAL VALUES OF R AND LPERM C C TYPJAC IS 1 IF THE APPROXIMATE DERIVATIVE IS THE DIAGONAL C MATRIX WITH ENTRIES 1/DELTA(J), J = 1,...,N, C FOR THE INITIAL GRID SIZES DELTA(.) C 2 IF IT IS THE NEGATIVE OF THE ABOVE C 3 IF IT IS THE SAME AS IN THE LAST MAJOR CYCLE C 4 IF IT IS GENERAL C C IF TYPJAC.EQ.4, JACFAC IS ZERO IF ADERIV IS NOT FACTORED. C OTHERWISE, ADERIV(.,LPERM(.)) IS SOME ORTHOGONAL MATRIX C TIMES THE MATRIX R, WHERE C C JACFAC IS 1 IF R HAS SEVERAL SUBDIAGONAL ENTRIES ZERO C 2 IF R IS UPPER HESSENBERG C 3 IF R IS UPPER TRIANGULAR C C SIGMA AND WA ARE SCRATCH ARRAYS USED IN QRGEN C C SUBROUTINES AND FUNCTIONS USED C QRGEN,QRHESS, DBLE C INTEGER I,J,L,IIN,KIN,KINSG,IOUT,KOUT,KOUTSG DOUBLE PRECISION TEMP,ZERO,ONE, DBLE ,DNEFF1 DATA ZERO,ONE /0.0D0,1.0D0/ C C BEGIN PROGRAM C IF (TYPJAC.GT.2) GO TO 200 C C ADERIV IS A SCALED IDENTITY MATRIX. SET UP BAS AND R. C TEMP = ADERIV(1,1)*DELTA(1) DO 120 J = 1,N DO 110 I = 1,NP1 BAS(I,J) = ZERO R(I,J) = ZERO 110 CONTINUE BAS(J,J) = TEMP R(J,J) = TEMP LPERM(J) = J 120 CONTINUE GO TO 500 C C SET UP BAS. C 200 DO 220 L = 1,N TEMP = DELTA(L) DO 210 I = 1,N BAS(I,L) = ADERIV(I,L)*TEMP 210 CONTINUE BAS(NP1,L) = ZERO 220 CONTINUE IF (TYPJAC.GT.3) GO TO 300 C C IF ADERIV UNCHANGED SINCE LAST CYCLE, SET UP R. C L = 1 DO 240 J = 1,N DO 230 I = 1,J R(I,J) = FR(L)*STPDWN L = L + 1 230 CONTINUE LPERM(J) = FLPERM(J) 240 CONTINUE GO TO 500 300 IF (JACFAC.GT.0) GO TO 400 C C FACTORIZE BAS TO GET R. C DO 340 J = 1,N DO 330 I = 1,NP1 R(I,J) = BAS(I,J) 330 CONTINUE LPERM(J) = J 340 CONTINUE CALL QRGEN(N,N,R,NP1,.TRUE.,LPERM,SIGMA,WA) GO TO 500 400 DO 420 J = 1,N L = LPERM(J) TEMP = DELTA(L) DO 410 I = 1,N R(I,J) = R(I,J)*TEMP 410 CONTINUE 420 CONTINUE IF (JACFAC.EQ.1) CALL QRGEN(N,N,R,NP1,.FALSE.,LPERM, * SIGMA,WA) IF (JACFAC.EQ.2) CALL QRHESS(R,NP1,N,1) C C SET UP LAST COLUMNS OF BAS AND R, WEIGHT, FR AND FLPERM, C AND POINTERS. C 500 DO 510 I = 1,N BAS(I,NP1) = ZERO R(I,NP1) = ZERO 510 CONTINUE IF (NEFFP1.EQ.0) NEFFP1 = NP1 DNEFF1 = DBLE(NEFFP1) DO 530 J = 1,N L = LPERM(J) TEMP = DBLE(NEFFP1 - L)/DNEFF1 IF (TEMP.LT.ZERO) TEMP = ZERO XORIG(L) = XORIG(L) - TEMP*DELTA(L) IF (VDLT.GE.1) TEMP = ZERO WEIGHT(J) = TEMP DO 520 I = 1,N BAS(I,NP1) = BAS(I,NP1) - TEMP*BAS(I,L) R(I,NP1) = R(I,NP1) - TEMP*R(I,J) 520 CONTINUE KPOINT(J) = L IPOINT(J) = 0 FLPERM(J) = L 530 CONTINUE BAS(NP1,NP1) = SCALE R(NP1,NP1) = SCALE LPERM(NP1) = NP1 IPOINT(NP1) = 0 KPOINT(NP1) = 0 WEIGHT(NP1) = ONE L = 1 DO 550 J = 1,N DO 540 I = 1,J FR(L) = R(I,J) L = L + 1 540 CONTINUE IF (R(J,J).EQ.ZERO) FLAG = 2 550 CONTINUE IF (R(NP1,NP1).EQ.ZERO) FLAG = 2 IF(PRTLEV.GE.5) CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) RETURN C C LAST CARD OF FSTBAS FOLLOWS C END SUBROUTINE NEWBAS(N,NP1,PRTLEV,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * NEWCOL,UPDCOL,RCOL,LPERM,IPOINT,KPOINT,MOSTV, * IIN,KIN,KINSG,IOUT,KOUT,KOUTSG,OUTTOP, * NBOUND,NTRIV,TYPFUN,SIGMA,WA,IWA,LIWA) INTEGER N,NP1,PRTLEV,LIWA,LPERM(NP1),IPOINT(NP1),KPOINT(NP1) INTEGER IIN,KIN,KINSG,IOUT,KOUT,KOUTSG INTEGER NBOUND,NTRIV,TYPFUN,MOSTV(NP1),IWA(LIWA) DOUBLE PRECISION BAS(NP1,NP1),R(NP1,NP1),WEIGHT(NP1) DOUBLE PRECISION ADERIV(N,N),DELTA(N),SCALE DOUBLE PRECISION NEWCOL(NP1),UPDCOL(NP1),RCOL(NP1) DOUBLE PRECISION SIGMA(NP1),WA(NP1) LOGICAL OUTTOP C C NEWBAS BRINGS THE NEW COLUMN NEWCOL INTO THE CURRENT BASIS C BAS AND DETERMINES THE LEAVING COLUMN C C ON ENTRY, C C BAS, R, LPERM, IPOINT, KPOINT DETERMINE THE CURRENT BASIS, ITS C FACTORIZATION, AND ITS POINTERS C C WEIGHT IS THE CURRENT RIGHT HAND SIDE VECTOR C C NEWCOL IS THE NP1 - VECTOR ENTERING THE BASIS C C IF IIN.NE.0, NEWCOL CORRESPONDS TO THE IIN'TH VERTEX OF THE C CURRENT PIECE AND KIN.EQ.0 C C IF KIN.NE.0, NEWCOL CORRESPONDS TO THE KIN'TH COLUMN OF ADERIV C AND IIN.EQ.0. KINSG IS 1 ( - 1) IF THE COLUMN ENTERS AT A POSITIVE C (NEGATIVE) LEVEL C C ON RETURN, C C BAS, R, LPERM, IPOINT, KPOINT, AND WEIGHT ARE SUITABLY UPDATED C C IF (OUTTOP) THE SPECIAL COLUMN HAS LEFT THE BASIS AND A MAJOR C CYCLE IS COMPLETED C C IF (.NOT.OUTTOP) IOUT (IF .NE.0) IS THE INDEX OF THE LEAVING C VERTEX, AND KOUT (IF .NE.0) AND KOUTSG DETERMINE THE LEAVING C VERTICES (SEE NEWPCE) C C SUBROUTINES AND FUNCTIONS USED C QRHESS,SOLVE,MINRAT, DBLE ,DABS,DMAX1 C DOUBLE PRECISION LEV,DKINSG,DKOUTS, DBLE ,DABS,DMAX1,ZERO,ONE DOUBLE PRECISION TWO,DIFF,TEMP,WJ,WN,WNP1 INTEGER I,J,JOUT,JOUTSG,L,LPNP1,LOUT,NM1,JP1,JOUTM1,PRT2,FLAG INTEGER INFO,MO DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ C C BEGIN PROGRAM C C CHECK FOR TRIVIAL PIVOT, WITH INCOMING COLUMN IDENTICAL TO C BASIS COLUMN C IF ((IIN.EQ.0).OR.(TYPFUN.LT.2)) GO TO 30 MO = MOSTV(IIN) IF (MO.LE.0) GO TO 30 DO 10 J = 1,NP1 I = IPOINT(J) IF (I.EQ.0) GO TO 10 IF (MOSTV(I).EQ.MO) GO TO 20 10 CONTINUE GO TO 30 C C TRIVIAL PIVOT C 20 CONTINUE IOUT = I KOUT = 0 IPOINT(J) = IIN OUTTOP = .FALSE. NTRIV = NTRIV + 1 GO TO 600 C C SOLVE BAS*P*UPDCOL = NEWCOL AND SET RCOL = R*UPDCOL C 30 CALL SOLVE(BAS,NP1,NP1,LPERM,R,NP1,NEWCOL,UPDCOL,RCOL,NEWCOL, * .FALSE.) IF (KIN.EQ.0) GO TO 50 IF (KINSG.EQ.1) GO TO 50 DO 40 I = 1,NP1 UPDCOL(I) = - UPDCOL(I) 40 CONTINUE 50 CALL MINRAT(N,NP1,IPOINT,KPOINT,KIN,WEIGHT,UPDCOL,JOUT,JOUTSG, * LEV,BAS,LPERM,R,SIGMA,WA,IWA,LIWA,INFO,SCALE) PRT2 = PRTLEV IF (INFO.EQ.0) GO TO 60 PRT2 = PRTLEV + 2 IF (PRT2.GE.5) CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) 60 CONTINUE C C COLUMN JOUT LEAVES C IF (JOUT.EQ.NP1) GO TO 200 OUTTOP = .FALSE. IOUT = IPOINT(JOUT) KOUT = KPOINT(JOUT) IF (KOUT.NE.0) KOUTSG = - JOUTSG C C UPDATE LAST COLUMN OF BAS, R C LOUT = LPERM(JOUT) IF (KIN.EQ.0) GO TO 80 DKINSG = DBLE(KINSG) DO 70 I = 1,N BAS(I,NP1) = BAS(I,NP1) + DKINSG*NEWCOL(I) R(I,NP1) = R(I,NP1) + DKINSG*RCOL(I) 70 CONTINUE R(NP1,NP1) = R(NP1,NP1) + DKINSG*RCOL(NP1) 80 IF (KOUT.EQ.0) GO TO 100 DKOUTS = DBLE(KOUTSG) DO 90 I = 1,N BAS(I,NP1) = BAS(I,NP1) - DKOUTS*BAS(I,LOUT) R(I,NP1) = R(I,NP1) - DKOUTS*R(I,JOUT) 90 CONTINUE 100 IF (JOUT.EQ.N) GO TO 130 C C ROTATE COLUMNS JOUT + 1,...,N OF R AND THEIR POINTERS C TO THE LEFT C NM1 = N - 1 DO 120 J = JOUT,NM1 JP1 = J + 1 LPERM(J) = LPERM(JP1) IPOINT(J) = IPOINT(JP1) KPOINT(J) = KPOINT(JP1) DO 110 I = 1,JP1 R(I,J) = R(I,JP1) 110 CONTINUE 120 CONTINUE LPERM(N) = LOUT 130 IPOINT(N) = IIN KPOINT(N) = KIN C C INSERT RCOL INTO R AND NEWCOL INTO BAS C DO 140 I = 1,NP1 R(I,N) = RCOL(I) BAS(I,LOUT) = NEWCOL(I) 140 CONTINUE C C UPDATE R C CALL QRHESS(R,NP1,NP1,JOUT) C C UPDATE WEIGHTS C LEV = DMAX1(LEV,ZERO) WNP1 = WEIGHT(NP1) - LEV*UPDCOL(NP1) JOUTM1 = JOUT - 1 IF (JOUTM1.EQ.0) GO TO 160 DO 150 J = 1,JOUTM1 WJ = WEIGHT(J) - LEV*UPDCOL(J) IF (IPOINT(J).GT.0) WJ = DMAX1(WJ,ZERO) WEIGHT(J) = WJ IF (IPOINT(J).GT.0) GO TO 150 WNP1 = DMAX1(WNP1,DABS(WJ)) 150 CONTINUE 160 IF (JOUT.EQ.N) GO TO 180 NM1 = N - 1 DO 170 J = JOUT,NM1 JP1 = J + 1 WJ = WEIGHT(JP1) - LEV*UPDCOL(JP1) IF (IPOINT(J).GT.0) WJ = DMAX1(WJ,ZERO) WEIGHT(J) = WJ IF (IPOINT(J).EQ.0) WNP1 = DMAX1(WNP1,DABS(WJ)) 170 CONTINUE 180 IF (KIN.EQ.0) WN = LEV IF (KIN.NE.0) WN = DKINSG*(LEV - WNP1) IF (KIN.NE.0) WNP1 = DMAX1(WNP1,DABS(WN)) WEIGHT(N) = WN WEIGHT(NP1) = WNP1 GO TO 400 C 200 IF (KIN.EQ.0) GO TO 300 C C ENTERING COLUMN HAS HIT BOUND C UPDATE LAST COLUMN OF BAS AND R C OUTTOP = .FALSE. NBOUND = NBOUND + 1 TEMP = TWO* DBLE(KINSG) DO 210 I = 1,N BAS(I,NP1) = BAS(I,NP1) + TEMP*NEWCOL(I) R(I,NP1) = R(I,NP1) + TEMP*RCOL(I) 210 CONTINUE R(NP1,NP1) = R(NP1,NP1) + TEMP*RCOL(NP1) C C UPDATE WEIGHTS C LEV = DMAX1(LEV,ZERO) WNP1 = WEIGHT(NP1) - LEV*UPDCOL(NP1) DO 220 J = 1,N WJ = WEIGHT(J) - LEV*UPDCOL(J) IF (IPOINT(J).GT.0) WJ = DMAX1(WJ,ZERO) WEIGHT(J) = WJ IF (IPOINT(J).GT.0) GO TO 220 WNP1 = DMAX1(WNP1,DABS(WJ)) 220 CONTINUE WEIGHT(NP1) = WNP1 IOUT = 0 KOUT = KIN KOUTSG = - KINSG GO TO 400 C C MAJOR CYCLE ENDED C 300 OUTTOP = .TRUE. IOUT = 0 KOUT = 0 IPOINT(NP1) = IIN DO 310 I = 1,NP1 BAS(I,NP1) = NEWCOL(I) R(I,NP1) = RCOL(I) 310 CONTINUE LEV = DMAX1(LEV,ZERO) WEIGHT(NP1) = LEV DO 320 J = 1,N WEIGHT(J) = DMAX1(ZERO,WEIGHT(J) - LEV*UPDCOL(J)) 320 CONTINUE C C CHECK WHETHER TO RECOMPUTE WEIGHTS C 400 TEMP = ZERO DO 410 J = 1,NP1 IF (KPOINT(J).EQ.0) TEMP = TEMP + WEIGHT(J) 410 CONTINUE DIFF = ONE - TEMP IF (ONE + DIFF*DIFF.EQ.ONE) GO TO 600 C C RECOMPUTE WEIGHTS C FLAG = 0 420 DO 430 I = 1,N SIGMA(I) = ZERO 430 CONTINUE SIGMA(NP1) = SCALE CALL SOLVE(BAS,NP1,NP1,LPERM,R,NP1,SIGMA,WEIGHT,WEIGHT,WA,.TRUE.) TEMP = ZERO WNP1 = WEIGHT(NP1) DO 450 J = 1,N WJ = WEIGHT(J) IF (IPOINT(J).EQ.0) GO TO 440 WJ = DMAX1(WJ,ZERO) WEIGHT(J) = WJ TEMP = TEMP + WJ GO TO 450 440 WNP1 = DMAX1(WNP1,DABS(WJ)) 450 CONTINUE TEMP = TEMP + WNP1 WEIGHT(NP1) = WNP1 DIFF = ONE - TEMP IF (ONE + DIFF*DIFF.EQ.ONE) GO TO 600 C C NUMERICAL DIFFICULTIES. TRY REFACTORIZING BASIS C IF (FLAG.EQ.0) GO TO 500 PRT2 = PRTLEV + 2 OUTTOP = .TRUE. WRITE (6,460) 460 FORMAT(53H0***** NUMERICAL PROBLEMS UNSOLVABLE BY REFACTORIZING ) GO TO 600 500 FLAG = 1 WRITE (6,510) 510 FORMAT (35H0***** REFACTORIZING ) IF (PRT2.GE.5) CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) PRT2 = PRTLEV + 2 DO 530 J = 1,NP1 L = LPERM(J) LPNP1 = L + NP1 IWA(L) = IPOINT(J) IWA(LPNP1) = KPOINT(J) DO 520 I = 1,NP1 R(I,J) = BAS(I,J) 520 CONTINUE 530 CONTINUE CALL QRGEN(NP1,NP1,R,NP1,.TRUE.,LPERM,SIGMA,WA) DO 540 J = 1,NP1 L = LPERM(J) LPNP1 = L + NP1 IPOINT(J) = IWA(L) KPOINT(J) = IWA(LPNP1) 540 CONTINUE GO TO 420 600 IF (PRT2.GE.5)CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) RETURN C C LAST CARD OF NEWBAS FOLLOWS C END SUBROUTINE LABEL(N,NP1,NEFFP1,IIN,KIN,XNEW,MOSTV,ADERIV,DELTA, * SCALE,NEWCOL,BIG,TOLFUN,TYPFUN,FLAG,ITOP,VDLT, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT, * IFACT) INTEGER N,NP1,NEFFP1,IIN,KIN,MOSTV(NP1),TYPFUN,FLAG,ITOP,VDLT DOUBLE PRECISION XNEW(N),ADERIV(N,N),DELTA(N),NEWCOL(NP1) DOUBLE PRECISION SCALE,BIG,TOLFUN INTEGER IDAT(1000),IFACT(*) DOUBLE PRECISION XDAT(1000),FX(*),FXORIG(*),FF(*),FM(*) DOUBLE PRECISION FDIAG(*),FDELTA,FB(*),FA(*) C C LABEL DETERMINES THE NEW COLUMN TO BE INTRODUCED INTO THE BASIS C C ON ENTRY, C C IIN IF NONZERO IS THE INDEX OF THE ENTERING VERTEX C C KIN IF NONZERO IS THE INDEX OF THE COORDINATE OF THE C ENTERING VERTICES (SEE SUBDIV) C C XNEW IF IIN.NE.ZERO IS THE NEW POINT TO BE EVALUATED C C MOSTV(I) IS 0 IF VERTEX I SPECIAL, 1 OTHERWISE C C ADERIV IS THE APPROXIMATE DERIVATIVE MATRIX C C DELTA(J) IS THE SCALE FACTOR IN THE JTH COORDINATE C C SCALE IS THE EXPECTED SIZE OF ENTRIES IN THE BASIS C C ON RETURN, C C NEWCOL IS THE COLUMN TO ENTER THE BASIS C C MOSTV(IIN) IS 0 IF IIN'TH VERTEX SPECIAL, 1 OTHERWISE C C FLAG IS 1 IF NORM(F(XNEW)).LT.TOLFUN C 2 IF NORM(XNEW).GT.BIG AND C 0 OTHERWISE. C C SUBROUTINES USED -- FNEVL1, VSNORM C INTEGER INFO,J DOUBLE PRECISION ZERO,TEMP,VSNORM DATA ZERO/0.0D0/ C C BEGIN PROGRAM C IF (KIN.NE.0) GO TO 30 INFO = 0 CALL FNEVL1(N,XNEW,NEWCOL,0,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) NEWCOL(NP1) = SCALE MOSTV(IIN) = INFO FLAG = 0 TEMP = VSNORM(N,NEWCOL) IF (VSNORM(N,XNEW).GT.BIG) FLAG = 2 IF (TEMP.LT.TOLFUN) FLAG = 1 IF (FLAG.EQ.1) WRITE (6,100) IF (FLAG.EQ.2) WRITE (6,200) IF ((FLAG.GT.0).OR.(TYPFUN.GT.0)) GO TO 50 DO 20 J = 1,N NEWCOL(J) = NEWCOL(J)/TEMP 20 CONTINUE GO TO 50 30 CONTINUE DO 40 J = 1,N NEWCOL(J) = DELTA(KIN)*ADERIV(J,KIN) 40 CONTINUE NEWCOL(NP1) = ZERO 50 RETURN C C FORMATS C 100 FORMAT (35H0***** STOPPED WITH GOOD VERTEX ) 200 FORMAT (45H0***** GENERATED VERTEX TOO FAR AWAY ) C C LAST CARD OF LABEL FOLLOWS C END SUBROUTINE MINRAT(N,NP1,IPOINT,KPOINT,KIN,W,U,JOUT,JOUTSG,LEV, * BAS,LPERM,R,SIGMA,WA,EXIT,LEXIT,INFO,SCALE) INTEGER N,NP1,IPOINT(NP1),KPOINT(NP1),LPERM(NP1),KIN,JOUT,JOUTSG INTEGER LEXIT,EXIT(LEXIT),INFO DOUBLE PRECISION W(NP1),U(NP1),LEV,SIGMA(NP1),WA(NP1),SCALE DOUBLE PRECISION BAS(NP1,NP1),R(NP1,NP1) C C MINRAT PERFORMS A TASK SIMILAR TO THE MINIMUM RATIO TEST OF C LINEAR PROGRAMMING C C ON ENTRY, C C W IS AN NP1 - VECTOR SATISFYING THE FOLLOWING INEQUALITIES C C W(J).GE.0 IF IPOINT(J).GT.0.OR.J.EQ.NP1 AND C C W(NP1).GE.W(J).GE.( - W(NP1)) OTHERWISE C C THE SUBROUTINE DETERMINES THE MAXIMUM LAMBDA SUCH THAT THESE C INEQUALITIES (AND LAMBDA.LE.2*(NEW W)(NP1) IF KIN.GT.0) HOLD C ALSO FOR NEW W = OLD W + LAMBDA * U. C C LEXICOGRAPHIC RULES ARE USED IN CASE OF DEGENERACY. C C C ON RETURN, C C LEV IS THIS MAXIMUM LAMBDA C C JOUT IS THE INDEX PROHIBITING FURTHER INCREASE. JOUT IS C SET TO NP1 IF KIN.GT.0 AND LAMBDA HITS ITS UPPER BOUND C C IF IPOINT(JOUT).EQ.0 JOUTSG IS 1 ( - 1) IF W(JOUT) HITS ITS C UPPER (LOWER) BOUND C C W IS NOT CHANGED C C FUNCTIONS USED -- MOD, DMAX1, DBLE , DABS C INTEGER FLAG,I,ICOL,ICOLM1,J,J1,K,NEXIT,NEXIT1,MOD DOUBLE PRECISION LAMBDA,DIFF,LEV1,SIGN,ZERO,ONE,TWO,HALF,HUGE DOUBLE PRECISION TEMP,TEMP2,WJ,WNP1,DMAX1,DABS,DBLE DATA ZERO,ONE,TWO,HALF,HUGE/0.0D0,1.0D0,2.0D0,0.5D0,1.0D30/ C C BEGIN PROGRAM C INFO = 0 5 LEV = TWO NEXIT = 0 DO 100 J = 1,N IF (IPOINT(J).EQ.0) GO TO 10 TEMP = U(J) TEMP2 = W(J) IF (TEMP.LE.TEMP2/TWO.OR.ONE + TEMP*TEMP.EQ.ONE) GO TO 100 LAMBDA = TEMP2/TEMP FLAG = 0 GO TO 30 10 TEMP = U(NP1) - U(J) TEMP2 = W(NP1) - W(J) IF (TEMP.LE.TEMP2/TWO.OR.ONE + TEMP*TEMP.EQ.ONE) GO TO 20 LAMBDA = TEMP2/TEMP FLAG = 1 GO TO 30 20 TEMP = U(NP1) + U(J) TEMP2 = W(NP1) + W(J) IF (TEMP.LE.TEMP2/TWO.OR.ONE + TEMP*TEMP.EQ.ONE) GO TO 100 LAMBDA = TEMP2/TEMP FLAG = 2 30 IF (DMAX1(LAMBDA,LEV).EQ.ZERO) GO TO 35 DIFF = (LAMBDA - LEV)/DMAX1(LAMBDA,LEV) IF ((ONE + DIFF*DIFF).GT.ONE) GO TO 40 35 NEXIT = NEXIT + 1 GO TO 50 40 IF (DIFF.GT.ZERO) GO TO 60 NEXIT = 1 50 EXIT(NEXIT) = J + NP1*FLAG IF (DIFF.LT.ZERO) LEV = LAMBDA 60 IF (FLAG.EQ.1) GO TO 20 100 CONTINUE IF (KIN.NE.0) GO TO 120 TEMP = U(NP1) TEMP2 = W(NP1) IF (TEMP.LE.TEMP2/TWO.OR.ONE + TEMP*TEMP.EQ.ONE) GO TO 200 LAMBDA = TEMP2/TEMP FLAG = 1 GO TO 150 120 TEMP = U(NP1) + HALF TEMP2 = W(NP1) IF (TEMP.LE.TEMP2/TWO.OR.ONE + TEMP*TEMP.EQ.ONE) GO TO 200 LAMBDA = TEMP2/TEMP FLAG = 2 150 IF (DMAX1(LAMBDA,LEV).EQ.ZERO) GO TO 155 DIFF = (LAMBDA - LEV)/DMAX1(LAMBDA,LEV) IF ((ONE + DIFF*DIFF).GT.ONE) GO TO 160 155 IF (KIN.EQ.0) GO TO 165 NEXIT = NEXIT + 1 GO TO 170 160 IF (DIFF.GT.ZERO) GO TO 200 165 NEXIT = 1 170 EXIT(NEXIT) = NP1*FLAG IF (DIFF.LT.ZERO) LEV = LAMBDA 200 CONTINUE IF (NEXIT.EQ.1) GO TO 500 C C LEXICOGRAPHIC RESOLUTION OF DEGENERACY. C FIRST RESOLVE FOR THE WEIGHTS W AND TRY AGAIN. C IF (INFO.EQ.1) GO TO 300 INFO = 1 DO 210 I = 1,NP1 SIGMA(I) = ZERO 210 CONTINUE SIGMA(NP1) = SCALE CALL SOLVE(BAS,NP1,NP1,LPERM,R,NP1,SIGMA,W,W,WA,.TRUE.) WNP1 = W(NP1) DO 220 J = 1,NP1 WJ = W(J) IF (IPOINT(J).GT.0) W(J) = DMAX1(ZERO,WJ) IF (IPOINT(J).GT.0) GO TO 220 WNP1 = DMAX1(WNP1,DABS(WJ)) 220 CONTINUE W(NP1) = WNP1 GO TO 5 300 DO 400 ICOL = 1,N IF (ICOL.EQ.1) GO TO 310 ICOLM1 = ICOL - 1 DO 305 I = 1,ICOLM1 SIGMA(I) = ZERO 305 CONTINUE 310 CONTINUE DO 315 I = ICOL,NP1 SIGMA(I) = DBLE(I) 315 CONTINUE CALL SOLVE(BAS,NP1,NP1,LPERM,R,NP1,SIGMA,WA,WA,SIGMA,.TRUE.) LEV1 = HUGE NEXIT1 = NEXIT NEXIT = 0 DO 390 K = 1,NEXIT1 J = EXIT(K) IF (J.GT.NP1) GO TO 320 TEMP = U(J) TEMP2 = WA(J) IF (DABS(TEMP).LE.DABS(TEMP2)/HUGE.OR.ONE + TEMP*TEMP.EQ.ONE) * GO TO 390 LAMBDA = TEMP2/TEMP GO TO 350 320 IF (J.NE.2*NP1) GO TO 330 TEMP = U(NP1) + HALF TEMP2 = WA(J) IF (DABS(TEMP).LE.DABS(TEMP2)/HUGE.OR.ONE + TEMP*TEMP.EQ.ONE) * GO TO 390 LAMBDA = TEMP2/TEMP GO TO 350 330 SIGN = ONE IF (J.GT.2*NP1) SIGN = - ONE J1 = MOD(J,NP1) TEMP = U(NP1) - SIGN*U(J1) TEMP2 = WA(NP1) - SIGN*WA(J1) IF (DABS(TEMP).LE.DABS(TEMP2)/HUGE.OR.ONE + TEMP*TEMP.EQ.ONE) * GO TO 390 LAMBDA = TEMP2/TEMP 350 IF (DMAX1(DABS(LAMBDA),DABS(LEV1)).EQ.ZERO) GO TO 355 DIFF = (LAMBDA - LEV1)/DMAX1(DABS(LAMBDA),DABS(LEV1)) IF ((ONE + DIFF*DIFF).GT.ONE) GO TO 360 355 NEXIT = NEXIT + 1 GO TO 370 360 IF (DIFF.GT.ZERO) GO TO 390 NEXIT = 1 370 EXIT(NEXIT) = J IF (DIFF.LT.ZERO) LEV1 = LAMBDA 390 CONTINUE IF (NEXIT.LE.1) GO TO 450 400 CONTINUE C C ONLY ONE CANDIDATE LEFT TO LEAVE THE BASIS. C 450 CONTINUE I = ICOL + 1 WRITE (6,460) I 460 FORMAT (35H0***** LEXICOGRAPHIC RULES USED -- ,I3,8H COLUMNS, * 15H EXAMINED ) IF (NEXIT.EQ.1) GO TO 500 WRITE (6,470) 470 FORMAT (50H0***** NUMERICAL PROBLEMS IN LEXICOGRAPHIC PIVOT ) 500 JOUT = EXIT(1) IF (JOUT.LE.NP1) GO TO 520 JOUTSG = - 1 IF (JOUT.LE.2*NP1) JOUTSG = 1 520 CONTINUE JOUT = MOD(JOUT,NP1) IF (JOUT.EQ.0) JOUT = NP1 RETURN C C LAST CARD OF MINRAT FOLLOWS C END SUBROUTINE PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) INTEGER N,NP1,LPERM(NP1),IPOINT(NP1),KPOINT(NP1) INTEGER IIN,KIN,KINSG,IOUT,KOUT,KOUTSG DOUBLE PRECISION BAS(NP1,NP1),R(NP1,NP1),WEIGHT(NP1) DOUBLE PRECISION ADERIV(N,N),DELTA(N),SCALE DOUBLE PRECISION WA(NP1) C C PRTBAS PRINTS THE APPROXIMATE DERIVATIVE MATRIX AND THE C COORDINATE SCALING DELTA. THEN THE BASIS MATRIX BAS AND THE C BASIS SCALING FACTOR SCALE. ALSO THE UPPER TRIANGULAR C FACTOR R, THE PERMUTATION LPERM AND THE POINTERS IPOINT C AND KPOINT. FINALLY THE ENTERING VERTEX (VERTICES) GIVEN C BY IIN, KIN, KINSG AND THE LEAVING ONE (ONES) GIVEN BY C IOUT,KOUT, KOUTSG. C INTEGER I,J,L DOUBLE PRECISION ZERO DATA ZERO/0.0D0/ C BEGIN PROGRAM C WRITE (6,100) DO 10 I = 1,N WRITE (6,200) (ADERIV(I,J),J = 1,N) 10 CONTINUE WRITE (6,300) WRITE (6,200) (DELTA(J),J = 1,N) WRITE (6,400) DO 20 I = 1,NP1 WRITE (6,200) (BAS(I,J),J = 1,NP1) 20 CONTINUE WRITE (6,500) SCALE WRITE (6,600) DO 30 I = 1,NP1 WRITE (6,200) (R(I,J),J = 1,NP1) 30 CONTINUE WRITE (6,700) WRITE (6,800) (LPERM(J),J = 1,NP1) WRITE (6,850) WRITE (6,200) (WEIGHT(J),J = 1,NP1) WRITE (6,900) WRITE (6,800) (IPOINT(J),J = 1,NP1) WRITE (6,1000) WRITE (6,800) (KPOINT(J),J = 1,NP1) WRITE (6,1100) IIN,KIN,KINSG,IOUT,KOUT,KOUTSG DO 35 I = 1,NP1 WA(I) = ZERO 35 CONTINUE DO 50 I = 1,NP1 DO 40 J = 1,NP1 L = LPERM(J) WA(I) = WA(I) + BAS(I,L)*WEIGHT(J) 40 CONTINUE 50 CONTINUE WRITE (6,1200) WRITE (6,200) (WA(I),I = 1,NP1) C C FORMATS C 100 FORMAT (45H0 THE APPROXIMATE DERIVATIVE MATRIX IS /) 200 FORMAT (6X,7D10.3) 300 FORMAT (45H0 THE COORDINATE SCALINGS ARE /) 400 FORMAT (45H0 THE BASIS MATRIX IS /) 500 FORMAT (35H0 THE BASIS SCALING FACTOR IS ,D10.3) 600 FORMAT (45H0 THE UPPER TRIANGULAR MATRIX R IS /) 700 FORMAT (45H0 THE PERMUTATION LPERM IS ) 800 FORMAT (1H0,5X,7(I5,5X)) 850 FORMAT (25H0 THE WEIGHTS ARE /) 900 FORMAT (30H0 THE POINTER IPOINT IS ) 1000 FORMAT (30H0 THE POINTER KPOINT IS ) 1100 FORMAT (10H0 IIN= ,I3,5H KIN= ,I3,7H KINSG= ,I3, * 6H IOUT= ,I3,6H KOUT= ,I3,8H KOUTSG= ,I3) 1200 FORMAT (30H0 BAS * P * WEIGHT = ) RETURN C C LAST CARD OF PRTBAS FOLLOWS C END