Chris@22: Chris@22: double precision FUNCTION DIFF(X,Y) Chris@22: c Chris@22: c Function used in tests that depend on machine precision. Chris@22: c Chris@22: c The original version of this code was developed by Chris@22: c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory Chris@22: c 1973 JUN 7, and published in the book Chris@22: c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Chris@22: c Revised FEB 1995 to accompany reprinting of the book by SIAM. Chris@22: C Chris@22: double precision X, Y Chris@22: DIFF=X-Y Chris@22: RETURN Chris@22: END Chris@22: Chris@22: Chris@22: SUBROUTINE G1 (A,B,CTERM,STERM,SIG) Chris@22: c Chris@22: C COMPUTE ORTHOGONAL ROTATION MATRIX.. Chris@22: c Chris@22: c The original version of this code was developed by Chris@22: c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory Chris@22: c 1973 JUN 12, and published in the book Chris@22: c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Chris@22: c Revised FEB 1995 to accompany reprinting of the book by SIAM. Chris@22: C Chris@22: C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) Chris@22: C (-S,C) (-S,C)(B) ( 0 ) Chris@22: C COMPUTE SIG = SQRT(A**2+B**2) Chris@22: C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT Chris@22: C SIG MAY BE IN THE SAME LOCATION AS A OR B . Chris@22: C ------------------------------------------------------------------ Chris@22: double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO Chris@22: parameter(ONE = 1.0d0, ZERO = 0.0d0) Chris@22: C ------------------------------------------------------------------ Chris@22: if (abs(A) .gt. abs(B)) then Chris@22: XR=B/A Chris@22: YR=sqrt(ONE+XR**2) Chris@22: CTERM=sign(ONE/YR,A) Chris@22: STERM=CTERM*XR Chris@22: SIG=abs(A)*YR Chris@22: RETURN Chris@22: endif Chris@22: Chris@22: if (B .ne. ZERO) then Chris@22: XR=A/B Chris@22: YR=sqrt(ONE+XR**2) Chris@22: STERM=sign(ONE/YR,B) Chris@22: CTERM=STERM*XR Chris@22: SIG=abs(B)*YR Chris@22: RETURN Chris@22: endif Chris@22: Chris@22: SIG=ZERO Chris@22: CTERM=ZERO Chris@22: STERM=ONE Chris@22: RETURN Chris@22: END Chris@22: Chris@22: Chris@22: C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) Chris@22: C Chris@22: C CONSTRUCTION AND/OR APPLICATION OF A SINGLE Chris@22: C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B Chris@22: C Chris@22: c The original version of this code was developed by Chris@22: c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory Chris@22: c 1973 JUN 12, and published in the book Chris@22: c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Chris@22: c Revised FEB 1995 to accompany reprinting of the book by SIAM. Chris@22: C ------------------------------------------------------------------ Chris@22: c Subroutine Arguments Chris@22: c Chris@22: C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a Chris@22: c Householder transformation, or Algorithm H2 to apply a Chris@22: c previously constructed transformation. Chris@22: C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. Chris@22: C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO Chris@22: C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M Chris@22: C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. Chris@22: C U(),IUE,UP On entry with MODE = 1, U() contains the pivot Chris@22: c vector. IUE is the storage increment between elements. Chris@22: c On exit when MODE = 1, U() and UP contain quantities Chris@22: c defining the vector U of the Householder transformation. Chris@22: c on entry with MODE = 2, U() and UP should contain Chris@22: c quantities previously computed with MODE = 1. These will Chris@22: c not be modified during the entry with MODE = 2. Chris@22: C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH Chris@22: c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE Chris@22: c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. Chris@22: c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. Chris@22: C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). Chris@22: C ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). Chris@22: C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 Chris@22: C NO OPERATIONS WILL BE DONE ON C(). Chris@22: C ------------------------------------------------------------------ Chris@22: SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) Chris@22: C ------------------------------------------------------------------ Chris@22: integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J Chris@22: integer L1, LPIVOT, M, MODE, NCV Chris@22: double precision B, C(*), CL, CLINV, ONE, SM Chris@22: c double precision U(IUE,M) Chris@22: double precision U(IUE,*) Chris@22: double precision UP Chris@22: parameter(ONE = 1.0d0) Chris@22: C ------------------------------------------------------------------ Chris@22: IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN Chris@22: CL=abs(U(1,LPIVOT)) Chris@22: IF (MODE.EQ.2) GO TO 60 Chris@22: C ****** CONSTRUCT THE TRANSFORMATION. ****** Chris@22: DO 10 J=L1,M Chris@22: 10 CL=MAX(abs(U(1,J)),CL) Chris@22: IF (CL) 130,130,20 Chris@22: 20 CLINV=ONE/CL Chris@22: SM=(U(1,LPIVOT)*CLINV)**2 Chris@22: DO 30 J=L1,M Chris@22: 30 SM=SM+(U(1,J)*CLINV)**2 Chris@22: CL=CL*SQRT(SM) Chris@22: IF (U(1,LPIVOT)) 50,50,40 Chris@22: 40 CL=-CL Chris@22: 50 UP=U(1,LPIVOT)-CL Chris@22: U(1,LPIVOT)=CL Chris@22: GO TO 70 Chris@22: C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** Chris@22: C Chris@22: 60 IF (CL) 130,130,70 Chris@22: 70 IF (NCV.LE.0) RETURN Chris@22: B= UP*U(1,LPIVOT) Chris@22: C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. Chris@22: C Chris@22: IF (B) 80,130,130 Chris@22: 80 B=ONE/B Chris@22: I2=1-ICV+ICE*(LPIVOT-1) Chris@22: INCR=ICE*(L1-LPIVOT) Chris@22: DO 120 J=1,NCV Chris@22: I2=I2+ICV Chris@22: I3=I2+INCR Chris@22: I4=I3 Chris@22: SM=C(I2)*UP Chris@22: DO 90 I=L1,M Chris@22: SM=SM+C(I3)*U(1,I) Chris@22: 90 I3=I3+ICE Chris@22: IF (SM) 100,120,100 Chris@22: 100 SM=SM*B Chris@22: C(I2)=C(I2)+SM*UP Chris@22: DO 110 I=L1,M Chris@22: C(I4)=C(I4)+SM*U(1,I) Chris@22: 110 I4=I4+ICE Chris@22: 120 CONTINUE Chris@22: 130 RETURN Chris@22: END Chris@22: Chris@22: Chris@22: Chris@22: C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) Chris@22: C Chris@22: C Algorithm NNLS: NONNEGATIVE LEAST SQUARES Chris@22: C Chris@22: c The original version of this code was developed by Chris@22: c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory Chris@22: c 1973 JUN 15, and published in the book Chris@22: c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Chris@22: c Revised FEB 1995 to accompany reprinting of the book by SIAM. Chris@22: c Chris@22: C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN Chris@22: C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM Chris@22: C Chris@22: C A * X = B SUBJECT TO X .GE. 0 Chris@22: C ------------------------------------------------------------------ Chris@22: c Subroutine Arguments Chris@22: c Chris@22: C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE Chris@22: C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N Chris@22: C MATRIX, A. ON EXIT A() CONTAINS Chris@22: C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN Chris@22: C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY Chris@22: C THIS SUBROUTINE. Chris@22: C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- Chris@22: C TAINS Q*B. Chris@22: C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL Chris@22: C CONTAIN THE SOLUTION VECTOR. Chris@22: C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE Chris@22: C RESIDUAL VECTOR. Chris@22: C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN Chris@22: C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. Chris@22: C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z Chris@22: C ZZ() AN M-ARRAY OF WORKING SPACE. Chris@22: C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. Chris@22: C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS Chris@22: C P AND Z AS FOLLOWS.. Chris@22: C Chris@22: C INDEX(1) THRU INDEX(NSETP) = SET P. Chris@22: C INDEX(IZ1) THRU INDEX(IZ2) = SET Z. Chris@22: C IZ1 = NSETP + 1 = NPP1 Chris@22: C IZ2 = N Chris@22: C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING Chris@22: C MEANINGS. Chris@22: C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. Chris@22: C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. Chris@22: C EITHER M .LE. 0 OR N .LE. 0. Chris@22: C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. Chris@22: C Chris@22: C ------------------------------------------------------------------ Chris@22: SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) Chris@22: C ------------------------------------------------------------------ Chris@22: integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L Chris@22: integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY Chris@22: c integer INDEX(N) Chris@22: c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) Chris@22: integer INDEX(*) Chris@22: double precision A(MDA,*), B(*), W(*), X(*), ZZ(*) Chris@22: double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM Chris@22: double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX Chris@22: double precision ZERO, ZTEST Chris@22: parameter(FACTOR = 0.01d0) Chris@22: parameter(TWO = 2.0d0, ZERO = 0.0d0) Chris@22: C ------------------------------------------------------------------ Chris@22: MODE=1 Chris@22: IF (M .le. 0 .or. N .le. 0) then Chris@22: MODE=2 Chris@22: RETURN Chris@22: endif Chris@22: ITER=0 Chris@22: ITMAX=3*N Chris@22: C Chris@22: C INITIALIZE THE ARRAYS INDEX() AND X(). Chris@22: C Chris@22: DO 20 I=1,N Chris@22: X(I)=ZERO Chris@22: 20 INDEX(I)=I Chris@22: C Chris@22: IZ2=N Chris@22: IZ1=1 Chris@22: NSETP=0 Chris@22: NPP1=1 Chris@22: C ****** MAIN LOOP BEGINS HERE ****** Chris@22: 30 CONTINUE Chris@22: C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. Chris@22: C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. Chris@22: C Chris@22: IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350 Chris@22: C Chris@22: C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). Chris@22: C Chris@22: DO 50 IZ=IZ1,IZ2 Chris@22: J=INDEX(IZ) Chris@22: SM=ZERO Chris@22: DO 40 L=NPP1,M Chris@22: 40 SM=SM+A(L,J)*B(L) Chris@22: W(J)=SM Chris@22: 50 continue Chris@22: C FIND LARGEST POSITIVE W(J). Chris@22: 60 continue Chris@22: WMAX=ZERO Chris@22: DO 70 IZ=IZ1,IZ2 Chris@22: J=INDEX(IZ) Chris@22: IF (W(J) .gt. WMAX) then Chris@22: WMAX=W(J) Chris@22: IZMAX=IZ Chris@22: endif Chris@22: 70 CONTINUE Chris@22: C Chris@22: C IF WMAX .LE. 0. GO TO TERMINATION. Chris@22: C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. Chris@22: C Chris@22: IF (WMAX .le. ZERO) go to 350 Chris@22: IZ=IZMAX Chris@22: J=INDEX(IZ) Chris@22: C Chris@22: C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. Chris@22: C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID Chris@22: C NEAR LINEAR DEPENDENCE. Chris@22: C Chris@22: ASAVE=A(NPP1,J) Chris@22: CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0) Chris@22: UNORM=ZERO Chris@22: IF (NSETP .ne. 0) then Chris@22: DO 90 L=1,NSETP Chris@22: 90 UNORM=UNORM+A(L,J)**2 Chris@22: endif Chris@22: UNORM=sqrt(UNORM) Chris@22: IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then Chris@22: C Chris@22: C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ Chris@22: C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). Chris@22: C Chris@22: DO 120 L=1,M Chris@22: 120 ZZ(L)=B(L) Chris@22: CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1) Chris@22: ZTEST=ZZ(NPP1)/A(NPP1,J) Chris@22: C Chris@22: C SEE IF ZTEST IS POSITIVE Chris@22: C Chris@22: IF (ZTEST .gt. ZERO) go to 140 Chris@22: endif Chris@22: C Chris@22: C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. Chris@22: C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL Chris@22: C COEFFS AGAIN. Chris@22: C Chris@22: A(NPP1,J)=ASAVE Chris@22: W(J)=ZERO Chris@22: GO TO 60 Chris@22: C Chris@22: C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM Chris@22: C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER Chris@22: C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN Chris@22: C COL J, SET W(J)=0. Chris@22: C Chris@22: 140 continue Chris@22: DO 150 L=1,M Chris@22: 150 B(L)=ZZ(L) Chris@22: C Chris@22: INDEX(IZ)=INDEX(IZ1) Chris@22: INDEX(IZ1)=J Chris@22: IZ1=IZ1+1 Chris@22: NSETP=NPP1 Chris@22: NPP1=NPP1+1 Chris@22: C Chris@22: IF (IZ1 .le. IZ2) then Chris@22: DO 160 JZ=IZ1,IZ2 Chris@22: JJ=INDEX(JZ) Chris@22: CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1) Chris@22: 160 continue Chris@22: endif Chris@22: C Chris@22: IF (NSETP .ne. M) then Chris@22: DO 180 L=NPP1,M Chris@22: 180 A(L,J)=ZERO Chris@22: endif Chris@22: C Chris@22: W(J)=ZERO Chris@22: C SOLVE THE TRIANGULAR SYSTEM. Chris@22: C STORE THE SOLUTION TEMPORARILY IN ZZ(). Chris@22: RTNKEY = 1 Chris@22: GO TO 400 Chris@22: 200 CONTINUE Chris@22: C Chris@22: C ****** SECONDARY LOOP BEGINS HERE ****** Chris@22: C Chris@22: C ITERATION COUNTER. Chris@22: C Chris@22: 210 continue Chris@22: ITER=ITER+1 Chris@22: IF (ITER .gt. ITMAX) then Chris@22: MODE=3 Chris@25: c write (*,'(/a)') ' NNLS quitting on iteration count.' Chris@22: GO TO 350 Chris@22: endif Chris@22: C Chris@22: C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. Chris@22: C IF NOT COMPUTE ALPHA. Chris@22: C Chris@22: ALPHA=TWO Chris@22: DO 240 IP=1,NSETP Chris@22: L=INDEX(IP) Chris@22: IF (ZZ(IP) .le. ZERO) then Chris@22: T=-X(L)/(ZZ(IP)-X(L)) Chris@22: IF (ALPHA .gt. T) then Chris@22: ALPHA=T Chris@22: JJ=IP Chris@22: endif Chris@22: endif Chris@22: 240 CONTINUE Chris@22: C Chris@22: C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL Chris@22: C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. Chris@22: C Chris@22: IF (ALPHA.EQ.TWO) GO TO 330 Chris@22: C Chris@22: C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO Chris@22: C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. Chris@22: C Chris@22: DO 250 IP=1,NSETP Chris@22: L=INDEX(IP) Chris@22: X(L)=X(L)+ALPHA*(ZZ(IP)-X(L)) Chris@22: 250 continue Chris@22: C Chris@22: C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I Chris@22: C FROM SET P TO SET Z. Chris@22: C Chris@22: I=INDEX(JJ) Chris@22: 260 continue Chris@22: X(I)=ZERO Chris@22: C Chris@22: IF (JJ .ne. NSETP) then Chris@22: JJ=JJ+1 Chris@22: DO 280 J=JJ,NSETP Chris@22: II=INDEX(J) Chris@22: INDEX(J-1)=II Chris@22: CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II)) Chris@22: A(J,II)=ZERO Chris@22: DO 270 L=1,N Chris@22: IF (L.NE.II) then Chris@22: c Chris@22: c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) Chris@22: c Chris@22: TEMP = A(J-1,L) Chris@22: A(J-1,L) = CC*TEMP + SS*A(J,L) Chris@22: A(J,L) =-SS*TEMP + CC*A(J,L) Chris@22: endif Chris@22: 270 CONTINUE Chris@22: c Chris@22: c Apply procedure G2 (CC,SS,B(J-1),B(J)) Chris@22: c Chris@22: TEMP = B(J-1) Chris@22: B(J-1) = CC*TEMP + SS*B(J) Chris@22: B(J) =-SS*TEMP + CC*B(J) Chris@22: 280 continue Chris@22: endif Chris@22: c Chris@22: NPP1=NSETP Chris@22: NSETP=NSETP-1 Chris@22: IZ1=IZ1-1 Chris@22: INDEX(IZ1)=I Chris@22: C Chris@22: C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD Chris@22: C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. Chris@22: C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY Chris@22: C THAT ARE NONPOSITIVE WILL BE SET TO ZERO Chris@22: C AND MOVED FROM SET P TO SET Z. Chris@22: C Chris@22: DO 300 JJ=1,NSETP Chris@22: I=INDEX(JJ) Chris@22: IF (X(I) .le. ZERO) go to 260 Chris@22: 300 CONTINUE Chris@22: C Chris@22: C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. Chris@22: C Chris@22: DO 310 I=1,M Chris@22: 310 ZZ(I)=B(I) Chris@22: RTNKEY = 2 Chris@22: GO TO 400 Chris@22: 320 CONTINUE Chris@22: GO TO 210 Chris@22: C ****** END OF SECONDARY LOOP ****** Chris@22: C Chris@22: 330 continue Chris@22: DO 340 IP=1,NSETP Chris@22: I=INDEX(IP) Chris@22: 340 X(I)=ZZ(IP) Chris@22: C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. Chris@22: GO TO 30 Chris@22: C Chris@22: C ****** END OF MAIN LOOP ****** Chris@22: C Chris@22: C COME TO HERE FOR TERMINATION. Chris@22: C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. Chris@22: C Chris@22: 350 continue Chris@22: SM=ZERO Chris@22: IF (NPP1 .le. M) then Chris@22: DO 360 I=NPP1,M Chris@22: 360 SM=SM+B(I)**2 Chris@22: else Chris@22: DO 380 J=1,N Chris@22: 380 W(J)=ZERO Chris@22: endif Chris@22: RNORM=sqrt(SM) Chris@22: RETURN Chris@22: C Chris@22: C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE Chris@22: C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). Chris@22: C Chris@22: 400 continue Chris@22: DO 430 L=1,NSETP Chris@22: IP=NSETP+1-L Chris@22: IF (L .ne. 1) then Chris@22: DO 410 II=1,IP Chris@22: ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1) Chris@22: 410 continue Chris@22: endif Chris@22: JJ=INDEX(IP) Chris@22: ZZ(IP)=ZZ(IP)/A(IP,JJ) Chris@22: 430 continue Chris@22: go to (200, 320), RTNKEY Chris@22: END Chris@22: