C ALGORITHM 587, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL. 8, NO. 3, C SEP., 1982, P.323. SUBROUTINE LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, LSEI 10 * MODE, WS, IP) C C DIMENSION W(MDW,N+1),PRGOPT(*),X(N), C WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2) C ABOVE, K=MAX(MA+MG,N). C C ABSTRACT C C THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES C PROBLEM WITH BOTH EQUALITY AND INEQUALITY CONSTRAINTS, AND, IF THE C USER REQUESTS, OBTAINS A COVARIANCE MATRIX OF THE SOLUTION C PARAMETERS. C C SUPPOSE THERE ARE GIVEN MATRICES E, A AND G OF RESPECTIVE C DIMENSIONS ME BY N, MA BY N AND MG BY N, AND VECTORS F, B AND H OF C RESPECTIVE LENGTHS ME, MA AND MG. THIS SUBROUTINE SOLVES THE C LINEARLY CONSTRAINED LEAST SQUARES PROBLEM C C EX = F, (E ME BY N) (EQUATIONS TO BE EXACTLY C SATISFIED) C AX = B, (A MA BY N) (EQUATIONS TO BE C APPROXIMATELY SATISFIED, C LEAST SQUARES SENSE) C GX.GE.H,(G MG BY N) (INEQUALITY CONSTRAINTS) C C THE INEQUALITIES GX.GE.H MEAN THAT EVERY COMPONENT OF THE PRODUCT C GX MUST BE .GE. THE CORRESPONDING COMPONENT OF H. C C IN CASE THE EQUALITY CONSTRAINTS CANNOT BE SATISFIED, A C GENERALIZED INVERSE SOLUTION RESIDUAL VECTOR LENGTH IS OBTAINED C FOR F-EX. THIS IS THE MINIMAL LENGTH POSSIBLE FOR F-EX. C C C ANY VALUES ME.GE.0, MA.GE.0, OR MG.GE.0 ARE PERMITTED. THE C RANK OF THE MATRIX E IS ESTIMATED DURING THE COMPUTATION. WE CALL C THIS VALUE KRANKE. IT IS AN OUTPUT PARAMETER IN IP(1) DEFINED C BELOW. USING A GENERALIZED INVERSE SOLUTION OF EX=F, A REDUCED C LEAST SQUARES PROBLEM WITH INEQUALITY CONSTRAINTS IS OBTAINED. C THE TOLERANCES USED IN THESE TESTS FOR DETERMINING THE RANK C OF E AND THE RANK OF THE REDUCED LEAST SQUARES PROBLEM ARE C GIVEN IN SANDIA TECH. REPT. SAND 78-1290. THEY CAN BE C MODIFIED BY THE USER IF NEW VALUES ARE PROVIDED IN C THE OPTION LIST OF THE ARRAY PRGOPT(*). C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START EDITING AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SASUM/DASUM/,/SDOT/DDOT/, C /SNRM2/DNRM2/,/ SQRT/ DSQRT/,/ ABS/ DABS/,/AMAX1/DMAX1/, C /SCOPY/DCOPY/,/SSCAL/DSCAL/,/SAXPY/DAXPY/,/SSWAP/DSWAP/,/E0/D0/, C /, DUMMY/,SNGL(DUMMY)/,/SRELPR/DRELPR/ C C WRITTEN BY R. J. HANSON AND K. H. HASKELL. FOR FURTHER MATH. C AND ALGORITHMIC DETAILS SEE SANDIA LABORATORIES TECH. REPTS. C SAND 77-0552, (1978), SAND 78-1290, (1979), AND C MATH. PROGRAMMING, VOL. 21, (1981), P.98-118. C C THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST.. C W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2) C WHERE K=MAX(MA+MG,N). THIS ALLOWS FOR A SOLUTION OF A RANGE OF C PROBLEMS IN THE GIVEN WORKING SPACE. THE DIMENSION OF WS(*) C GIVEN IS A NECESSARY OVERESTIMATE. ONCE A PARTICULAR PROBLEM C HAS BEEN RUN, THE OUTPUT PARAMETER IP(3) GIVES THE ACTUAL C DIMENSION REQUIRED FOR THAT PROBLEM. C C THE PARAMETERS FOR LSEI( ) ARE C C INPUT.. C C W(*,*),MDW, THE ARRAY W(*,*) IS DOUBLY SUBSCRIPTED WITH C ME,MA,MG,N FIRST DIMENSIONING PARAMETER EQUAL TO MDW. C FOR THIS DISCUSSION LET US CALL M = ME+MA+MG. THEN C MDW MUST SATISFY MDW.GE.M. THE CONDITION C MDW.LT.M IS AN ERROR. C C THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS C C (E F) C (A B) C (G H) C C IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1 C RESPECTIVELY. C C THE INTEGERS ME, MA, AND MG ARE THE C RESPECTIVE MATRIX ROW DIMENSIONS C OF E, A AND G. EACH MATRIX HAS N COLUMNS. C C PRGOPT(*) THIS ARRAY IS THE OPTION VECTOR. C IF THE USER IS SATISFIED WITH THE NOMINAL C SUBPROGRAM FEATURES SET C C PRGOPT(1)=1 (OR PRGOPT(1)=1.0) C C OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF C GROUPS OF DATA OF THE FOLLOWING FORM C C LINK C KEY C DATA SET C C THE PARAMETERS LINK AND KEY ARE EACH ONE WORD. C THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS. C THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY. C THE VALUE OF LINK POINTS TO THE FIRST C ENTRY OF THE NEXT GROUP OF DATA WITHIN C PRGOPT(*). THE EXCEPTION IS WHEN THERE ARE C NO MORE OPTIONS TO CHANGE. IN THAT C CASE LINK=1 AND THE VALUES KEY AND DATA SET C ARE NOT REFERENCED. THE GENERAL LAYOUT OF C PRGOPT(*) IS AS FOLLOWS. C C ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP) C . PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE) C . PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE) C . . C . . C . . C ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF C . NEXT GROUP) C . PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE) C . PRGOPT(LINK1+2)=DATA VALUE C ... . C . . C . . C ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE) C C VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS. C A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR. C THIS HELPS PREVENT USING INVALID BUT POSITIVE C VALUES OF LINK THAT WILL PROBABLY EXTEND C BEYOND THE PROGRAM LIMITS OF PRGOPT(*). C UNRECOGNIZED VALUES OF KEY ARE IGNORED. THE C ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER C OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING C RESTRICTION. TO PREVENT CYCLING IN THE C PROCESSING OF THE OPTION ARRAY A COUNT OF THE C NUMBER OF OPTIONS CHANGED IS MAINTAINED. C WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR C MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS. C C OPTIONS.. C C KEY=1 C COMPUTE IN W(*,*) THE N BY N C COVARIANCE MATRIX OF THE SOLUTION VARIABLES C AS AN OUTPUT PARAMETER. NOMINALLY THE C COVARIANCE MATRIX WILL NOT BE COMPUTED. C (THIS REQUIRES NO USER INPUT.) C THE DATA SET FOR THIS OPTION IS A SINGLE VALUE. C IT MUST BE NONZERO WHEN THE COVARIANCE MATRIX C IS DESIRED. IF IT IS ZERO, THE COVARIANCE C MATRIX IS NOT COMPUTED. WHEN THE COVARIANCE MATRIX C IS COMPUTED, THE FIRST DIMENSIONING PARAMETER C OF THE ARRAY W(*,*) MUST SATISFY MDW.GE.MAX0(M,N). C C KEY=2 C SCALE THE NONZERO COLUMNS OF THE C ENTIRE DATA MATRIX. C (E) C (A) C (G) C C TO HAVE LENGTH ONE. THE DATA SET FOR THIS C OPTION IS A SINGLE VALUE. IT MUST BE C NONZERO IF UNIT LENGTH COLUMN SCALING C IS DESIRED. C C KEY=3 C SCALE COLUMNS OF THE ENTIRE DATA MATRIX C (E) C (A) C (G) C C WITH A USER-PROVIDED DIAGONAL MATRIX. C THE DATA SET FOR THIS OPTION CONSISTS C OF THE N DIAGONAL SCALING FACTORS, ONE FOR C EACH MATRIX COLUMN. C C KEY=4 C CHANGE THE RANK DETERMINATION TOLERANCE FOR C THE EQUALITY CONSTRAINT EQUATIONS FROM C THE NOMINAL VALUE OF SQRT(SRELPR). THIS QUANTITY CAN C BE NO SMALLER THAN SRELPR, THE ARITHMETIC- C STORAGE PRECISION. THE QUANTITY SRELPR IS THE C LARGEST POSITIVE NUMBER SUCH THAT T=1.+SRELPR C SATISFIES T.EQ.1. THE QUANTITY USED C HERE IS INTERNALLY RESTRICTED TO BE AT C LEAST SRELPR. THE DATA SET FOR THIS OPTION C IS THE NEW TOLERANCE. C C KEY=5 C CHANGE THE RANK DETERMINATION TOLERANCE FOR C THE REDUCED LEAST SQUARES EQUATIONS FROM C THE NOMINAL VALUE OF SQRT(SRELPR). THIS QUANTITY CAN C BE NO SMALLER THAN SRELPR, THE ARITHMETIC- C STORAGE PRECISION. THE QUANTITY USED C HERE IS INTERNALLY RESTRICTED TO BE AT C LEAST SRELPR. THE DATA SET FOR THIS OPTION C IS THE NEW TOLERANCE. C C FOR EXAMPLE, SUPPOSE WE WANT TO CHANGE C THE TOLERANCE FOR THE REDUCED LEAST SQUARES C PROBLEM, COMPUTE THE COVARIANCE MATRIX OF C THE SOLUTION PARAMETERS, AND PROVIDE C COLUMN SCALING FOR THE DATA MATRIX. FOR C THESE OPTIONS THE DIMENSION OF PRGOPT(*) C MUST BE AT LEAST N+9. THE FORTRAN STATEMENTS C DEFINING THESE OPTIONS WOULD BE AS FOLLOWS. C C PRGOPT(1)=4 (LINK TO ENTRY 4 IN PRGOPT(*)) C PRGOPT(2)=1 (COVARIANCE MATRIX KEY) C PRGOPT(3)=1 (COVARIANCE MATRIX WANTED) C C PRGOPT(4)=7 (LINK TO ENTRY 7 IN PRGOPT(*)) C PRGOPT(5)=5 (LEAST SQUARES EQUAS. TOLERANCE KEY) C PRGOPT(6)=... (NEW VALUE OF THE TOLERANCE) C C PRGOPT(7)=N+9 (LINK TO ENTRY N+9 IN PRGOPT(*)) C PRGOPT(8)=3 (USER-PROVIDED COLUMN SCALING KEY) C C CALL SCOPY(N,D,1,PRGOPT(9),1) (COPY THE N C SCALING FACTORS FROM THE USER ARRAY D(*) C TO PRGOPT(9)-PRGOPT(N+8)) C C PRGOPT(N+9)=1 (NO MORE OPTIONS TO CHANGE) C C THE CONTENTS OF PRGOPT(*) ARE NOT MODIFIED C BY THE SUBPROGRAM. C THE OPTIONS FOR WNNLS( ) CAN ALSO BE INCLUDED C IN THIS ARRAY. THE VALUES OF KEY RECOGNIZED C BY WNNLS( ) ARE 6, 7 AND 8. THEIR FUNCTIONS C ARE DOCUMENTED IN THE USAGE INSTRUCTIONS FOR C SUBROUTINE WNNLS( ). NORMALLY THESE OPTIONS C DO NOT NEED TO BE MODIFIED WHEN USING LSEI( ). C C IP(1), THE AMOUNTS OF WORKING STORAGE ACTUALLY C IP(2) ALLOCATED FOR THE WORKING ARRAYS WS(*) AND C IP(*), RESPECTIVELY. THESE QUANTITIES ARE C COMPARED WITH THE ACTUAL AMOUNTS OF STORAGE C NEEDED BY LSEI( ). INSUFFICIENT STORAGE C ALLOCATED FOR EITHER WS(*) OR IP(*) IS AN C ERROR. THIS FEATURE WAS INCLUDED IN LSEI( ) C BECAUSE MISCALCULATING THE STORAGE FORMULAS C FOR WS(*) AND IP(*) MIGHT VERY WELL LEAD TO C SUBTLE AND HARD-TO-FIND EXECUTION ERRORS. C C THE LENGTH OF WS(*) MUST BE AT LEAST C C LW = 2*(ME+N)+K+(MG+2)*(N+7) C C C WHERE K = MAX(MA+MG,N) C THIS TEST WILL NOT BE MADE IF IP(1).LE.0. C C THE LENGTH OF IP(*) MUST BE AT LEAST C C LIP = MG+2*N+2 C THIS TEST WILL NOT BE MADE IF IP(2).LE.0. C C OUTPUT.. C C X(*),RNORME, THE ARRAY X(*) CONTAINS THE SOLUTION PARAMETERS C RNORML IF THE INTEGER OUTPUT FLAG MODE = 0 OR 1. C THE DEFINITION OF MODE IS GIVEN DIRECTLY BELOW. C WHEN MODE = 0 OR 1, RNORME AND RNORML C RESPECTIVELY CONTAIN THE RESIDUAL VECTOR C EUCLIDEAN LENGTHS OF F - EX AND B - AX. WHEN C MODE=1 THE EQUALITY CONSTRAINT EQUATIONS EX=F C ARE CONTRADICTORY, SO RNORME.NE.0. THE RESIDUAL C VECTOR F-EX HAS MINIMAL EUCLIDEAN LENGTH. FOR C MODE.GE.2, NONE OF THESE PARAMETERS ARE C DEFINED. C C MODE INTEGER FLAG THAT INDICATES THE SUBPROGRAM C STATUS AFTER COMPLETION. IF MODE.GE.2, NO C SOLUTION HAS BEEN COMPUTED. C C MODE = C C 0 BOTH EQUALITY AND INEQUALITY CONSTRAINTS C ARE COMPATIBLE AND HAVE BEEN SATISFIED. C C 1 EQUALITY CONSTRAINTS ARE CONTRADICTORY. C A GENERALIZED INVERSE SOLUTION OF EX=F WAS USED C TO MINIMIZE THE RESIDUAL VECTOR LENGTH F-EX. C IN THIS SENSE, THE SOLUTION IS STILL MEANINGFUL. C C 2 INEQUALITY CONSTRAINTS ARE CONTRADICTORY. C C 3 BOTH EQUALITY AND INEQUALITY CONSTRAINTS C ARE CONTRADICTORY. C C THE FOLLOWING INTERPRETATION OF C MODE=1,2 OR 3 MUST BE MADE. THE C SETS CONSISTING OF ALL SOLUTIONS C OF THE EQUALITY CONSTRAINTS EX=F C AND ALL VECTORS SATISFYING GX.GE.H C HAVE NO POINTS ON COMMON. (IN C PARTICULAR THIS DOES NOT SAY THAT C EACH INDIVIDUAL SET HAS NO POINTS C AT ALL, ALTHOUGH THIS COULD BE THE C CASE.) C C 4 USAGE ERROR OCCURRED. THE VALUE C OF MDW IS .LT. ME+MA+MG, MDW IS C .LT. N AND A COVARIANCE MATRIX IS C REQUESTED, THE OPTION VECTOR C PRGOPT(*) IS NOT PROPERLY DEFINED, C OR THE LENGTHS OF THE WORKING ARRAYS C WS(*) AND IP(*), WHEN SPECIFIED IN C IP(1) AND IP(2) RESPECTIVELY, ARE NOT C LONG ENOUGH. C C W(*,*) THE ARRAY W(*,*) CONTAINS THE N BY N SYMMETRIC C COVARIANCE MATRIX OF THE SOLUTION PARAMETERS, C PROVIDED THIS WAS REQUESTED ON INPUT WITH C THE OPTION VECTOR PRGOPT(*) AND THE OUTPUT C FLAG IS RETURNED WITH MODE = 0 OR 1. C C IP(*) THE INTEGER WORKING ARRAY HAS THREE ENTRIES C THAT PROVIDE RANK AND WORKING ARRAY LENGTH C INFORMATION AFTER COMPLETION. C C IP(1) = RANK OF EQUALITY CONSTRAINT C MATRIX. DEFINE THIS QUANTITY C AS KRANKE. C C IP(2) = RANK OF REDUCED LEAST SQUARES C PROBLEM. C C IP(3) = THE AMOUNT OF STORAGE IN THE C WORKING ARRAY WS(*) THAT WAS C ACTUALLY USED BY THE SUBPROGRAM. C THE FORMULA GIVEN ABOVE FOR THE LENGTH C OF WS(*) IS A NECESSARY OVERESTIMATE. C USER DESIGNATED C WORKING ARRAYS.. C C WS(*),IP(*) THESE ARE RESP. TYPE FLOATING POINT C AND TYPE INTEGER WORKING ARRAYS. C THEIR REQUIRED MINIMAL LENGTHS ARE C GIVEN ABOVE. C C C SUBROUTINES CALLED C C LSI PART OF THIS PACKAGE. SOLVES A C CONSTRAINED LEAST SQUARES PROBLEM WITH C INEQUALITY CONSTRAINTS. C C++ C SDOT,SSCAL, SUBROUTINES FROM THE BLAS PACKAGE. C SAXPY,SASUM, SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308. C SCOPY,SNRM2, C SSWAP C C H12 SUBROUTINE TO CONSTRUCT AND APPLY A C HOUSEHOLDER TRANSFORMATION. C C XERROR FROM SLATEC ERROR PROCESSING PACKAGE. C THIS IS DOCUMENTED IN SANDIA TECH. REPT., C SAND78-1189. C C SUBROUTINE LSEI(W,MDW,ME,MA,MG,N,PRGOPT,X,RNORME,RNORML,MODE,WS, C 1 IP) C C REVISED OCT. 1, 1981. C REAL W(MDW,1), PRGOPT(1), X(1), WS(1), RNORME, RNORML INTEGER IP(3) REAL DUMMY, ENORM, SRELPR, FNORM, GAM, HALF, ONE, RB, *RN, RNMAX, SIZE, SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME, ZERO REAL SASUM, SDOT, SNRM2, SQRT, ABS, AMAX1 LOGICAL COV DATA ZERO /0.E0/, SRELPR /0.E0/, HALF /0.5E0/, ONE /1.E0/ C C CHECK THAT ENOUGH STORAGE WAS ALLOCATED IN WS(*) AND IP(*). IF (.NOT.(IP(1).GT.0)) GO TO 20 LCHK = 2*(ME+N) + MAX0(MA+MG,N) + (MG+2)*(N+7) IF (.NOT.(IP(1).LT.LCHK)) GO TO 10 MODE = 4 NERR = 2 IOPT = 1 CALL XERRWV(67HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR WS(*), *NEED LW=I1 BELOW, 67, NERR, IOPT, 1, LCHK, 0, * 0, DUMMY, DUMMY) RETURN 10 CONTINUE 20 IF (.NOT.(IP(2).GT.0)) GO TO 40 LCHK = MG + 2*N + 2 IF (.NOT.(IP(2).LT.LCHK)) GO TO 30 MODE = 4 NERR = 2 IOPT = 1 CALL XERRWV(68HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR IP(*), *NEED LIP=I1 BELOW, 68, NERR, IOPT, 1, LCHK, 0, * 0, DUMMY, DUMMY) RETURN 30 CONTINUE C C COMPUTE MACHINE PRECISION=SRELPR ONLY WHEN NECESSARY. 40 IF (.NOT.(SRELPR.EQ.ZERO)) GO TO 70 SRELPR = ONE 50 IF (ONE+SRELPR.EQ.ONE) GO TO 60 SRELPR = SRELPR*HALF GO TO 50 60 SRELPR = SRELPR + SRELPR C C COMPUTE NUMBER OF POSSIBLE RIGHT MULTIPLYING HOUSEHOLDER C TRANSFORMATIONS. 70 M = ME + MA + MG MODE = 0 IF (N.LE.0 .OR. M.LE.0) RETURN IF (.NOT.(MDW.LT.M)) GO TO 80 NERR = 1 IOPT = 1 CALL XERROR(36HLSEI( ), MDW.LT.ME+MA+MG IS AN ERROR, 36, NERR, * IOPT) MODE = 4 RETURN 80 NP1 = N + 1 KRANKE = MIN0(ME,N) N1 = 2*KRANKE + 1 N2 = N1 + N C C PROCESS-OPTION-VECTOR ASSIGN 90 TO IGO990 GO TO 480 90 IF (.NOT.(COV .AND. MDW.LT.N)) GO TO 100 NERR = 2 IOPT = 1 CALL XERROR( * 54HLSEI( ), MDW.LT.N, WHEN COV MATRIX NEEDED, IS AN ERROR, 54, * NERR, IOPT) MODE = 4 RETURN 100 L = KRANKE C C COMPUTE NORM OF EQUALITY CONSTRAINT MATRIX AND RT SIDE. ENORM = ZERO DO 110 J=1,N ENORM = AMAX1(ENORM,SASUM(ME,W(1,J),1)) 110 CONTINUE FNORM = SASUM(ME,W(1,NP1),1) IF (.NOT.(L.GT.0)) GO TO 190 SNMAX = ZERO RNMAX = ZERO DO 180 I=1,L C C COMPUTE MAXIMUM RATIO OF VECTOR LENGTHS. PARTITION C IS AT COL. I. DO 150 K=I,ME SN = SDOT(N-I+1,W(K,I),MDW,W(K,I),MDW) RN = SDOT(I-1,W(K,1),MDW,W(K,1),MDW) IF (.NOT.(RN.EQ.ZERO .AND. SN.GT.SNMAX)) GO TO 120 SNMAX = SN IMAX = K GO TO 140 120 IF (.NOT.(K.EQ.I .OR. (SN*RNMAX.GT.RN*SNMAX))) GO TO 130 SNMAX = SN RNMAX = RN IMAX = K 130 CONTINUE 140 CONTINUE 150 CONTINUE C C INTERCHANGE ROWS IF NECESSARY. IF (I.NE.IMAX) CALL SSWAP(NP1, W(I,1), MDW, W(IMAX,1), MDW) IF (.NOT.(SNMAX.GT.TAU**2*RNMAX)) GO TO 160 C C ELIMINATE ELEMS I+1,...,N IN ROW I. CALL H12(1, I, I+1, N, W(I,1), MDW, WS(I), W(I+1,1), MDW, 1, * M-I) GO TO 170 160 KRANKE = I - 1 GO TO 200 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE C C SAVE DIAG. TERMS OF LOWER TRAP. MATRIX. CALL SCOPY(KRANKE, W, MDW+1, WS(KRANKE+1), 1) C C USE HOUSEHOLDER TRANS FROM LEFT TO ACHIEVE KRANKE BY KRANKE UPPER C TRIANGULAR FORM. IF (.NOT.(KRANKE.GT.0 .AND. KRANKE.LT.ME)) GO TO 220 DO 210 KK=1,KRANKE K = KRANKE + 1 - KK C C APPLY TRANFORMATION TO MATRIX COLS. 1,...,K-1. CALL H12(1, K, KRANKE+1, ME, W(1,K), 1, UP, W, 1, MDW, K-1) C C APPLY TO RT SIDE VECTOR. CALL H12(2, K, KRANKE+1, ME, W(1,K), 1, UP, W(1,NP1), 1, 1, 1) 210 CONTINUE 220 IF (.NOT.(KRANKE.GT.0)) GO TO 240 C C SOLVE FOR VARIABLES 1,...,KRANKE IN NEW COORDINATES. CALL SCOPY(KRANKE, W(1,NP1), 1, X, 1) DO 230 I=1,KRANKE X(I) = (X(I)-SDOT(I-1,W(I,1),MDW,X,1))/W(I,I) 230 CONTINUE C C COMPUTE RESIDUALS FOR REDUCED PROBLEM. 240 MEP1 = ME + 1 RNORML = ZERO IF (.NOT.(ME.LT.M)) GO TO 270 DO 260 I=MEP1,M W(I,NP1) = W(I,NP1) - SDOT(KRANKE,W(I,1),MDW,X,1) SN = SDOT(KRANKE,W(I,1),MDW,W(I,1),MDW) RN = SDOT(N-KRANKE,W(I,KRANKE+1),MDW,W(I,KRANKE+1),MDW) IF (.NOT.(RN.LE.TAU**2*SN .AND. KRANKE.LT.N)) GO TO 250 W(I,KRANKE+1) = ZERO CALL SCOPY(N-KRANKE, W(I,KRANKE+1), 0, W(I,KRANKE+1), MDW) 250 CONTINUE 260 CONTINUE C C COMPUTE EQUAL. CONSTRAINT EQUAS. RESIDUAL LENGTH. 270 RNORME = SNRM2(ME-KRANKE,W(KRANKE+1,NP1),1) C C MOVE REDUCED PROBLEM DATA UPWARD IF KRANKE.LT.ME. IF (.NOT.(KRANKE.LT.ME)) GO TO 290 DO 280 J=1,NP1 CALL SCOPY(M-ME, W(ME+1,J), 1, W(KRANKE+1,J), 1) 280 CONTINUE C C COMPUTE SOLN OF REDUCED PROBLEM. 290 CALL LSI(W(KRANKE+1,KRANKE+1), MDW, MA, MG, N-KRANKE, PRGOPT, * X(KRANKE+1), RNORML, MODE, WS(N2), IP(2)) IF (.NOT.(ME.GT.0)) GO TO 330 C C TEST FOR CONSISTENCY OF EQUALITY CONSTRAINTS. MDEQC = 0 XNRME = SASUM(KRANKE,W(1,NP1),1) IF (RNORME.GT.TAU*(ENORM*XNRME+FNORM)) MDEQC = 1 MODE = MODE + MDEQC C C CHECK IF SOLN TO EQUAL. CONSTRAINTS SATISFIES INEQUAL. C CONSTRAINTS WHEN THERE ARE NO DEGREES OF FREEDOM LEFT. IF (.NOT.(KRANKE.EQ.N .AND. MG.GT.0)) GO TO 320 XNORM = SASUM(N,X,1) MAPKE1 = MA + KRANKE + 1 MEND = MA + KRANKE + MG DO 310 I=MAPKE1,MEND SIZE = SASUM(N,W(I,1),MDW)*XNORM + ABS(W(I,NP1)) IF (.NOT.(W(I,NP1).GT.TAU*SIZE)) GO TO 300 MODE = MODE + 2 GO TO 450 300 CONTINUE 310 CONTINUE 320 CONTINUE 330 IF (.NOT.(KRANKE.GT.0)) GO TO 420 C C REPLACE DIAG. TERMS OF LOWER TRAP. MATRIX. CALL SCOPY(KRANKE, WS(KRANKE+1), 1, W, MDW+1) C C REAPPLY TRANS TO PUT SOLN IN ORIGINAL COORDINATES. DO 340 II=1,KRANKE I = KRANKE + 1 - II CALL H12(2, I, I+1, N, W(I,1), MDW, WS(I), X, 1, 1, 1) 340 CONTINUE C C COMPUTE COV MATRIX OF EQUAL. CONSTRAINED PROBLEM. IF (.NOT.(COV)) GO TO 410 DO 400 JJ=1,KRANKE J = KRANKE + 1 - JJ IF (.NOT.(J.LT.N)) GO TO 390 RB = WS(J)*W(J,J) IF (RB.NE.ZERO) RB = ONE/RB JP1 = J + 1 DO 350 I=JP1,N W(I,J) = SDOT(N-J,W(I,JP1),MDW,W(J,JP1),MDW)*RB 350 CONTINUE GAM = SDOT(N-J,W(JP1,J),1,W(J,JP1),MDW)*RB GAM = HALF*GAM CALL SAXPY(N-J, GAM, W(J,JP1), MDW, W(JP1,J), 1) DO 370 I=JP1,N DO 360 K=I,N W(I,K) = W(I,K) + W(J,I)*W(K,J) + W(I,J)*W(J,K) W(K,I) = W(I,K) 360 CONTINUE 370 CONTINUE UJ = WS(J) VJ = GAM*UJ W(J,J) = UJ*VJ + UJ*VJ DO 380 I=JP1,N W(J,I) = UJ*W(I,J) + VJ*W(J,I) 380 CONTINUE CALL SCOPY(N-J, W(J,JP1), MDW, W(JP1,J), 1) 390 CONTINUE 400 CONTINUE 410 CONTINUE C C APPLY THE SCALING TO THE COVARIANCE MATRIX. 420 IF (.NOT.(COV)) GO TO 440 DO 430 I=1,N L = N1 + I CALL SSCAL(N, WS(L-1), W(I,1), MDW) CALL SSCAL(N, WS(L-1), W(1,I), 1) 430 CONTINUE 440 CONTINUE 450 CONTINUE C C RESCALE SOLN. VECTOR. IF (.NOT.(MODE.LE.1)) GO TO 470 DO 460 J=1,N L = N1 + J X(J) = X(J)*WS(L-1) 460 CONTINUE 470 IP(1) = KRANKE IP(3) = IP(3) + 2*KRANKE + N RETURN 480 CONTINUE C TO PROCESS-OPTION-VECTOR C C THE NOMINAL TOLERANCE USED IN THE CODE C FOR THE EQUALITY CONSTRAINT EQUATIONS. TAU = SQRT(SRELPR) C C THE NOMINAL COLUMN SCALING USED IN THE CODE IS C THE IDENTITY SCALING. WS(N1) = ONE CALL SCOPY(N, WS(N1), 0, WS(N1), 1) C C NO COVARIANCE MATRIX IS NOMINALLY COMPUTED. COV = .FALSE. C C DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE. NOPT = 1000 NTIMES = 0 C C DEFINE BOUND FOR POSITIVE VALUES OF LINK. NLINK = 100000 LAST = 1 LINK = PRGOPT(1) IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 490 NERR = 3 IOPT = 1 CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR, * IOPT) MODE = 4 RETURN 490 IF (.NOT.(LINK.GT.1)) GO TO 540 NTIMES = NTIMES + 1 IF (.NOT.(NTIMES.GT.NOPT)) GO TO 500 NERR = 3 IOPT = 1 CALL XERROR( * 52HLSEI( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 52, * NERR, IOPT) MODE = 4 RETURN 500 KEY = PRGOPT(LAST+1) IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO IF (.NOT.(KEY.EQ.2 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 520 DO 510 J=1,N T = SNRM2(M,W(1,J),1) IF (T.NE.ZERO) T = ONE/T L = N1 + J WS(L-1) = T 510 CONTINUE 520 IF (KEY.EQ.3) CALL SCOPY(N, PRGOPT(LAST+2), 1, WS(N1), 1) IF (KEY.EQ.4) TAU = AMAX1(SRELPR,PRGOPT(LAST+2)) NEXT = PRGOPT(LINK) IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 530 NERR = 3 IOPT = 1 CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR, * IOPT) MODE = 4 RETURN 530 LAST = LINK LINK = NEXT GO TO 490 540 DO 550 J=1,N L = N1 + J CALL SSCAL(M, WS(L-1), W(1,J), 1) 550 CONTINUE GO TO 560 560 GO TO IGO990, (90) END SUBROUTINE LSI(W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, IP) LSI 10 C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START EDITING AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SASUM/DASUM/,/SDOT/DDOT/, C / SQRT/ DSQRT/,/AMAX1/DMAX1/,/SSWAP/DSWAP/, C /SCOPY/DCOPY/,/SSCAL/DSCAL/,/SAXPY/DAXPY/,/E0/D0/,/SRELPR/DRELPR/ C C THIS IS A COMPANION SUBPROGRAM TO LSEI( ). C THE DOCUMENTATION FOR LSEI( ) HAS MORE COMPLETE C USAGE INSTRUCTIONS. C WRITTEN BY R. J. HANSON, SLA. C C SOLVE.. C AX = B, A MA BY N (LEAST SQUARES EQUATIONS) C SUBJECT TO.. C C GX.GE.H, G MG BY N (INEQUALITY CONSTRAINTS) C C INPUT.. C C W(*,*) CONTAINS (A B) IN ROWS 1,...,MA+MG, COLS 1,...,N+1. C (G H) C C MDW,MA,MG,N C CONTAIN (RESP) VAR. DIMENSION OF W(*,*), C AND MATRIX DIMENSIONS. C C PRGOPT(*), C PROGRAM OPTION VECTOR. C C OUTPUT.. C C X(*),RNORM C C SOLUTION VECTOR(UNLESS MODE=2), LENGTH OF AX-B. C C MODE C =0 INEQUALITY CONSTRAINTS ARE COMPATIBLE. C =2 INEQUALITY CONSTRAINTS CONTRADICTORY. C C WS(*), C WORKING STORAGE OF DIMENSION K+N+(MG+2)*(N+7), C WHERE K=MAX(MA+MG,N). C IP(MG+2*N+1) C INTEGER WORKING STORAGE C REVISED OCT. 1, 1981. C C SUBROUTINES CALLED C C LPDP THIS SUBPROGRAM MINIMIZES A SUM OF SQUARES C OF UNKNOWNS SUBJECT TO LINEAR INEQUALITY C CONSTRAINTS. PART OF THIS PACKAGE. C C++ C SDOT,SSCAL SUBROUTINES FROM THE BLAS PACKAGE. C SAXPY,SASUM, SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308. C SCOPY,SSWAP C C HFTI SOLVES AN UNCONSTRAINED LINEAR LEAST SQUARES C PROBLEM. PART OF THIS PACKAGE. C C H12 SUBROUTINE TO CONSTRUCT AND APPLY A HOUSEHOLDER C TRANSFORMATION. C C SUBROUTINE LSI(W,MDW,MA,MG,N,PRGOPT,X,RNORM,MODE,WS,IP) C REAL W(MDW,1), PRGOPT(1), RNORM, WS(1), X(1) INTEGER IP(1) REAL ANORM, SRELPR, FAC, GAM, HALF, ONE, RB, TAU, TOL, * XNORM, ZERO REAL SASUM, SDOT, SQRT, AMAX1 LOGICAL COV C DATA ZERO /0.E0/, SRELPR /0.E0/, ONE /1.E0/, HALF /.5E0/ C C COMPUTE MACHINE PRECISION=SRELPR ONLY WHEN NECESSARY. IF (.NOT.(SRELPR.EQ.ZERO)) GO TO 30 SRELPR = ONE 10 IF (ONE+SRELPR.EQ.ONE) GO TO 20 SRELPR = SRELPR*HALF GO TO 10 20 SRELPR = SRELPR + SRELPR 30 MODE = 0 RNORM = ZERO M = MA + MG NP1 = N + 1 KRANK = 0 IF (N.LE.0 .OR. M.LE.0) GO TO 70 ASSIGN 40 TO IGO994 GO TO 500 C C PROCESS-OPTION-VECTOR C C COMPUTE MATRIX NORM OF LEAST SQUARES EQUAS. 40 ANORM = ZERO DO 50 J=1,N ANORM = AMAX1(ANORM,SASUM(MA,W(1,J),1)) 50 CONTINUE C C SET TOL FOR HFTI( ) RANK TEST. TAU = TOL*ANORM C C COMPUTE HOUSEHOLDER ORTHOGONAL DECOMP OF MATRIX. IF (N.GT.0) WS(1) = ZERO CALL SCOPY(N, WS, 0, WS, 1) CALL SCOPY(MA, W(1,NP1), 1, WS, 1) K = MAX0(M,N) MINMAN = MIN0(MA,N) N1 = K + 1 N2 = N1 + N CALL HFTI(W, MDW, MA, N, WS, 1, 1, TAU, KRANK, RNORM, WS(N2), * WS(N1), IP) FAC = ONE GAM=MA-KRANK IF (KRANK.LT.MA) FAC = RNORM**2/GAM ASSIGN 60 TO IGO990 GO TO 80 C C REDUCE-TO-LPDP-AND-SOLVE 60 CONTINUE 70 IP(1) = KRANK IP(2) = N + MAX0(M,N) + (MG+2)*(N+7) RETURN 80 CONTINUE C C TO REDUCE-TO-LPDP-AND-SOLVE MAP1 = MA + 1 C C COMPUTE INEQ. RT-HAND SIDE FOR LPDP. IF (.NOT.(MA.LT.M)) GO TO 260 IF (.NOT.(MINMAN.GT.0)) GO TO 160 DO 90 I=MAP1,M W(I,NP1) = W(I,NP1) - SDOT(N,W(I,1),MDW,WS,1) 90 CONTINUE DO 100 I=1,MINMAN J = IP(I) C C APPLY PERMUTATIONS TO COLS OF INEQ. CONSTRAINT MATRIX. CALL SSWAP(MG, W(MAP1,I), 1, W(MAP1,J), 1) 100 CONTINUE C C APPLY HOUSEHOLDER TRANSFORMATIONS TO CONSTRAINT MATRIX. IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 120 DO 110 II=1,KRANK I = KRANK + 1 - II L = N1 + I CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), W(MAP1,1), * MDW, 1, MG) 110 CONTINUE C C COMPUTE PERMUTED INEQ. CONSTR. MATRIX TIMES R-INVERSE. 120 DO 150 I=MAP1,M IF (.NOT.(0.LT.KRANK)) GO TO 140 DO 130 J=1,KRANK W(I,J) = (W(I,J)-SDOT(J-1,W(1,J),1,W(I,1),MDW))/W(J,J) 130 CONTINUE 140 CONTINUE 150 CONTINUE C C SOLVE THE REDUCED PROBLEM WITH LPDP ALGORITHM, C THE LEAST PROJECTED DISTANCE PROBLEM. 160 CALL LPDP(W(MAP1,1), MDW, MG, KRANK, N-KRANK, PRGOPT, X, XNORM, * MDLPDP, WS(N2), IP(N+1)) IF (.NOT.(MDLPDP.EQ.1)) GO TO 240 IF (.NOT.(KRANK.GT.0)) GO TO 180 C C COMPUTE SOLN IN ORIGINAL COORDINATES. DO 170 II=1,KRANK I = KRANK + 1 - II X(I) = (X(I)-SDOT(II-1,W(I,I+1),MDW,X(I+1),1))/W(I,I) 170 CONTINUE C C APPLY HOUSEHOLDER TRANS. TO SOLN VECTOR. 180 IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 200 DO 190 I=1,KRANK L = N1 + I CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), X, 1, 1, 1) 190 CONTINUE 200 IF (.NOT.(MINMAN.GT.0)) GO TO 230 C C REPERMUTE VARIABLES TO THEIR INPUT ORDER. DO 210 II=1,MINMAN I = MINMAN + 1 - II J = IP(I) CALL SSWAP(1, X(I), 1, X(J), 1) 210 CONTINUE C C VARIABLES ARE NOW IN ORIG. COORDINATES. C ADD SOLN OF UNSCONSTRAINED PROB. DO 220 I=1,N X(I) = X(I) + WS(I) 220 CONTINUE C C COMPUTE THE RESIDUAL VECTOR NORM. RNORM = SQRT(RNORM**2+XNORM**2) 230 GO TO 250 240 MODE = 2 250 GO TO 270 260 CALL SCOPY(N, WS, 1, X, 1) 270 IF (.NOT.(COV .AND. KRANK.GT.0)) GO TO 490 C C COMPUTE COVARIANCE MATRIX BASED ON THE ORTHOGONAL DECOMP. C FROM HFTI( ). C KRM1 = KRANK - 1 KRP1 = KRANK + 1 C C COPY DIAG. TERMS TO WORKING ARRAY. CALL SCOPY(KRANK, W, MDW+1, WS(N2), 1) C C RECIPROCATE DIAG. TERMS. DO 280 J=1,KRANK W(J,J) = ONE/W(J,J) 280 CONTINUE IF (.NOT.(KRANK.GT.1)) GO TO 310 C C INVERT THE UPPER TRIANGULAR QR FACTOR ON ITSELF. DO 300 I=1,KRM1 IP1 = I + 1 DO 290 J=IP1,KRANK W(I,J) = -SDOT(J-I,W(I,I),MDW,W(I,J),1)*W(J,J) 290 CONTINUE 300 CONTINUE C C COMPUTE THE INVERTED FACTOR TIMES ITS TRANSPOSE. 310 DO 330 I=1,KRANK DO 320 J=I,KRANK W(I,J) = SDOT(KRANK+1-J,W(I,J),MDW,W(J,J),MDW) 320 CONTINUE 330 CONTINUE IF (.NOT.(KRANK.LT.N)) GO TO 450 C C ZERO OUT LOWER TRAPEZOIDAL PART. C COPY UPPER TRI. TO LOWER TRI. PART. DO 340 J=1,KRANK CALL SCOPY(J, W(1,J), 1, W(J,1), MDW) 340 CONTINUE DO 350 I=KRP1,N W(I,1) = ZERO CALL SCOPY(I, W(I,1), 0, W(I,1), MDW) 350 CONTINUE C C APPLY RIGHT SIDE TRANSFORMATIONS TO LOWER TRI. N3 = N2 + KRP1 DO 430 I=1,KRANK L = N1 + I K = N2 + I RB = WS(L-1)*WS(K-1) IF (.NOT.(RB.LT.ZERO)) GO TO 420 C C IF RB.GE.ZERO, TRANSFORMATION CAN BE REGARDED AS ZERO. RB = ONE/RB C C STORE UNSCALED RANK-ONE HOUSEHOLDER UPDATE IN WORK ARRAY. WS(N3) = ZERO CALL SCOPY(N, WS(N3), 0, WS(N3), 1) L = N1 + I K = N3 + I WS(K-1) = WS(L-1) DO 360 J=KRP1,N K = N3 + J WS(K-1) = W(I,J) 360 CONTINUE DO 370 J=1,N L = N3 + I K = N3 + J WS(J) = SDOT(J-I,W(J,I),MDW,WS(L-1),1) + SDOT(N-J+1,W(J,J),1, * WS(K-1),1) WS(J) = WS(J)*RB 370 CONTINUE L = N3 + I GAM = SDOT(N-I+1,WS(L-1),1,WS(I),1)*RB GAM = GAM*HALF CALL SAXPY(N-I+1, GAM, WS(L-1), 1, WS(I), 1) DO 410 J=I,N IF (.NOT.(I.GT.1)) GO TO 390 IM1 = I - 1 K = N3 + J DO 380 L=1,IM1 W(J,L) = W(J,L) + WS(K-1)*WS(L) 380 CONTINUE 390 K = N3 + J DO 400 L=I,J IL = N3 + L W(J,L) = W(J,L) + WS(J)*WS(IL-1) + WS(L)*WS(K-1) 400 CONTINUE 410 CONTINUE 420 CONTINUE 430 CONTINUE C C COPY LOWER TRI. TO UPPER TRI. TO SYMMETRIZE THE COVARIANCE MATRIX. DO 440 I=1,N CALL SCOPY(I, W(I,1), MDW, W(1,I), 1) 440 CONTINUE C C REPERMUTE ROWS AND COLS. 450 DO 470 II=1,MINMAN I = MINMAN + 1 - II K = IP(I) IF (.NOT.(I.NE.K)) GO TO 460 CALL SSWAP(1, W(I,I), 1, W(K,K), 1) CALL SSWAP(I-1, W(1,I), 1, W(1,K), 1) CALL SSWAP(K-I-1, W(I,I+1), MDW, W(I+1,K), 1) CALL SSWAP(N-K, W(I,K+1), MDW, W(K,K+1), MDW) 460 CONTINUE 470 CONTINUE C C PUT IN NORMALIZED RESIDUAL SUM OF SQUARES SCALE FACTOR C AND SYMMETRIZE THE RESULTING COVARIANCE MARIX. DO 480 J=1,N CALL SSCAL(J, FAC, W(1,J), 1) CALL SCOPY(J, W(1,J), 1, W(J,1), MDW) 480 CONTINUE 490 GO TO 540 500 CONTINUE C C TO PROCESS-OPTION-VECTOR C C THE NOMINAL TOLERANCE USED IN THE CODE, TOL = SQRT(SRELPR) COV = .FALSE. LAST = 1 LINK = PRGOPT(1) 510 IF (.NOT.(LINK.GT.1)) GO TO 520 KEY = PRGOPT(LAST+1) IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO IF (KEY.EQ.5) TOL = AMAX1(SRELPR,PRGOPT(LAST+2)) NEXT = PRGOPT(LINK) LAST = LINK LINK = NEXT GO TO 510 520 GO TO 530 530 GO TO IGO994, (40) 540 GO TO IGO990, (60) END SUBROUTINE LPDP(A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, IS)LPDP 10 C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START EDITING AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SNRM2/DNRM2/,/SDOT/DDOT/, C /SCOPY/DCOPY/,/SSCAL/DSCAL/,/ABS(/DABS(/, ABS/, DABS/,/E0/D0/ C C DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1), C WHERE N=N1+N2. THIS IS A SLIGHT OVERESTIMATE FOR WS(*). C C WRITTEN BY R. J. HANSON AND K. H. HASKELL, SANDIA LABS C REVISED OCT. 1, 1981. C C DETERMINE AN N1-VECTOR W, AND C AN N2-VECTOR Z C WHICH MINIMIZES THE EUCLIDEAN LENGTH OF W C SUBJECT TO G*W+H*Z .GE. Y. C THIS IS THE LEAST PROJECTED DISTANCE PROBLEM, LPDP. C THE MATRICES G AND H ARE OF RESPECTIVE C DIMENSIONS M BY N1 AND M BY N2. C C CALLED BY SUBPROGRAM LSI( ). C C THE MATRIX C (G H Y) C C OCCUPIES ROWS 1,...,M AND COLS 1,...,N1+N2+1 OF A(*,*). C C THE SOLUTION (W) IS RETURNED IN X(*). C (Z) C C THE VALUE OF MODE INDICATES THE STATUS OF C THE COMPUTATION AFTER RETURNING TO THE USER. C C MODE=1 THE SOLUTION WAS SUCCESSFULLY OBTAINED. C C MODE=2 THE INEQUALITIES ARE INCONSISTENT. C C SUBROUTINES CALLED C C WNNLS SOLVES A NONNEGATIVELY CONSTRAINED LINEAR LEAST C SQUARES PROBLEM WITH LINEAR EQUALITY CONSTRAINTS. C PART OF THIS PACKAGE. C C++ C SDOT, SUBROUTINES FROM THE BLAS PACKAGE. C SSCAL,SNRM2, SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308. C SCOPY C REAL A(MDA,1), PRGOPT(1), WS(1), WNORM, X(1) INTEGER IS(1) REAL FAC, ONE, RNORM, SC, YNORM, ZERO REAL SDOT, SNRM2, ABS DATA ZERO, ONE /0.E0,1.E0/, FAC /0.1E0/ N = N1 + N2 MODE = 1 IF (.NOT.(M.LE.0)) GO TO 20 IF (.NOT.(N.GT.0)) GO TO 10 X(1) = ZERO CALL SCOPY(N, X, 0, X, 1) 10 WNORM = ZERO RETURN 20 NP1 = N + 1 C C SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE. DO 40 I=1,M SC = SNRM2(N,A(I,1),MDA) IF (.NOT.(SC.NE.ZERO)) GO TO 30 SC = ONE/SC CALL SSCAL(NP1, SC, A(I,1), MDA) 30 CONTINUE 40 CONTINUE C C SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO). YNORM = SNRM2(M,A(1,NP1),1) IF (.NOT.(YNORM.NE.ZERO)) GO TO 50 SC = ONE/YNORM CALL SSCAL(M, SC, A(1,NP1), 1) C C SCALE COLS OF MATRIX H. 50 J = N1 + 1 60 IF (.NOT.(J.LE.N)) GO TO 70 SC = SNRM2(M,A(1,J),1) IF (SC.NE.ZERO) SC = ONE/SC CALL SSCAL(M, SC, A(1,J), 1) X(J) = SC J = J + 1 GO TO 60 70 IF (.NOT.(N1.GT.0)) GO TO 130 C C COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*). IW = 0 DO 80 I=1,M C C MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY. CALL SCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1) IW = IW + N2 C C MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY. CALL SCOPY(N1, A(I,1), MDA, WS(IW+1), 1) IW = IW + N1 C C MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY. WS(IW+1) = A(I,NP1) IW = IW + 1 80 CONTINUE WS(IW+1) = ZERO CALL SCOPY(N, WS(IW+1), 0, WS(IW+1), 1) IW = IW + N WS(IW+1) = ONE IW = IW + 1 C C SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0. THE C MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR C F = TRANSPOSE OF (0,...,0,1). IX = IW + 1 IW = IW + M C C DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ). IS(1) = 0 IS(2) = 0 CALL WNNLS(WS, NP1, N2, NP1-N2, M, 0, PRGOPT, WS(IX), RNORM, * MODEW, IS, WS(IW+1)) C C COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W. SC = ONE - SDOT(M,A(1,NP1),1,WS(IX),1) IF (.NOT.(ONE+FAC*ABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 110 SC = ONE/SC DO 90 J=1,N1 X(J) = SC*SDOT(M,A(1,J),1,WS(IX),1) 90 CONTINUE C C COMPUTE THE VECTOR Q=Y-GW. OVERWRITE Y WITH THIS VECTOR. DO 100 I=1,M A(I,NP1) = A(I,NP1) - SDOT(N1,A(I,1),MDA,X,1) 100 CONTINUE GO TO 120 110 MODE = 2 RETURN 120 CONTINUE 130 IF (.NOT.(N2.GT.0)) GO TO 180 C C COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*). IW = 0 DO 140 I=1,M CALL SCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1) IW = IW + N2 WS(IW+1) = A(I,NP1) IW = IW + 1 140 CONTINUE WS(IW+1) = ZERO CALL SCOPY(N2, WS(IW+1), 0, WS(IW+1), 1) IW = IW + N2 WS(IW+1) = ONE IW = IW + 1 IX = IW + 1 IW = IW + M C C SOLVE RV=S SUBJECT TO V.GE.0. THE MATRIX R =(TRANSPOSE C OF (H Q)), WHERE Q=Y-GW. THE (N2+1)-VECTOR S =(TRANSPOSE C OF (0,...,0,1)). C C DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ). IS(1) = 0 IS(2) = 0 CALL WNNLS(WS, N2+1, 0, N2+1, M, 0, PRGOPT, WS(IX), RNORM, MODEW, * IS, WS(IW+1)) C C COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z. SC = ONE - SDOT(M,A(1,NP1),1,WS(IX),1) IF (.NOT.(ONE+FAC*ABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 160 SC = ONE/SC DO 150 J=1,N2 L = N1 + J X(L) = SC*SDOT(M,A(1,L),1,WS(IX),1)*X(L) 150 CONTINUE GO TO 170 160 MODE = 2 RETURN 170 CONTINUE C C ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION. 180 CALL SSCAL(N, YNORM, X, 1) WNORM = SNRM2(N1,X,1) RETURN END SUBROUTINE WNNLS(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, WNN 10 * IWORK, WORK) C C DIMENSION W(MDW,N+1),PRGOPT(*),X(N),IWORK(M+N),WORK(M+5*N) C C ABSTRACT C C THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES C PROBLEM. SUPPOSE THERE ARE GIVEN MATRICES E AND A OF C RESPECTIVE DIMENSIONS ME BY N AND MA BY N, AND VECTORS F C AND B OF RESPECTIVE LENGTHS ME AND MA. THIS SUBROUTINE C SOLVES THE PROBLEM C C EX = F, (EQUATIONS TO BE EXACTLY SATISFIED) C C AX = B, (EQUATIONS TO BE APPROXIMATELY SATISFIED, C IN THE LEAST SQUARES SENSE) C C SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE C C ANY VALUES ME.GE.0, MA.GE.0 AND 0.LE. L .LE.N ARE PERMITTED. C C THE PROBLEM IS REPOSED AS PROBLEM WNNLS C C (WT*E)X = (WT*F) C ( A) ( B), (LEAST SQUARES) C SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE. C C THE SUBPROGRAM CHOOSES THE HEAVY WEIGHT (OR PENALTY PARAMETER) WT. C C THE PARAMETERS FOR WNNLS ARE C C INPUT.. C C W(*,*),MDW, THE ARRAY W(*,*) IS DOUBLE SUBSCRIPTED WITH FIRST C ME,MA,N,L DIMENSIONING PARAMETER EQUAL TO MDW. FOR THIS C DISCUSSION LET US CALL M = ME + MA. THEN MDW C MUST SATISFY MDW.GE.M. THE CONDITION MDW.LT.M C IS AN ERROR. C C THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS C C (E F) C (A B) C C IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1 C RESPECTIVELY. COLUMNS 1,...,L CORRESPOND TO C UNCONSTRAINED VARIABLES X(1),...,X(L). THE C REMAINING VARIABLES ARE CONSTRAINED TO BE C NONNEGATIVE. THE CONDITION L.LT.0 .OR. L.GT.N IS C AN ERROR. C C PRGOPT(*) THIS ARRAY IS THE OPTION VECTOR. C IF THE USER IS SATISFIED WITH THE NOMINAL C SUBPROGRAM FEATURES SET C C PRGOPT(1)=1 (OR PRGOPT(1)=1.0) C C OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF C GROUPS OF DATA OF THE FOLLOWING FORM C C LINK C KEY C DATA SET C C THE PARAMETERS LINK AND KEY ARE EACH ONE WORD. C THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS. C THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY. C THE VALUE OF LINK POINTS TO THE FIRST C ENTRY OF THE NEXT GROUP OF DATA WITHIN C PRGOPT(*). THE EXCEPTION IS WHEN THERE ARE C NO MORE OPTIONS TO CHANGE. IN THAT C CASE LINK=1 AND THE VALUES KEY AND DATA SET C ARE NOT REFERENCED. THE GENERAL LAYOUT OF C PRGOPT(*) IS AS FOLLOWS. C C ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP) C . PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE) C . PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE) C . . C . . C . . C ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF C . NEXT GROUP) C . PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE) C . PRGOPT(LINK1+2)=DATA VALUE C ... . C . . C . . C ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE) C C VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS. C A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR. C THIS HELPS PREVENT USING INVALID BUT POSITIVE C VALUES OF LINK THAT WILL PROBABLY EXTEND C BEYOND THE PROGRAM LIMITS OF PRGOPT(*). C UNRECOGNIZED VALUES OF KEY ARE IGNORED. THE C ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER C OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING C RESTRICTION. TO PREVENT CYCLING IN THE C PROCESSING OF THE OPTION ARRAY A COUNT OF THE C NUMBER OF OPTIONS CHANGED IS MAINTAINED. C WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR C MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS. C C OPTIONS.. C C KEY=6 C SCALE THE NONZERO COLUMNS OF THE C ENTIRE DATA MATRIX C (E) C (A) C TO HAVE LENGTH ONE. THE DATA SET FOR C THIS OPTION IS A SINGLE VALUE. IT MUST C BE NONZERO IF UNIT LENGTH COLUMN SCALING IS C DESIRED. C C KEY=7 C SCALE COLUMNS OF THE ENTIRE DATA MATRIX C (E) C (A) C WITH A USER-PROVIDED DIAGONAL MATRIX. C THE DATA SET FOR THIS OPTION CONSISTS C OF THE N DIAGONAL SCALING FACTORS, ONE FOR C EACH MATRIX COLUMN. C C KEY=8 C CHANGE THE RANK DETERMINATION TOLERANCE FROM C THE NOMINAL VALUE OF SQRT(EPS). THIS QUANTITY CAN C BE NO SMALLER THAN EPS, THE ARITHMETIC- C STORAGE PRECISION. THE QUANTITY USED C HERE IS INTERNALLY RESTRICTED TO BE AT C LEAST EPS. THE DATA SET FOR THIS OPTION C IS THE NEW TOLERANCE. C C KEY=9 C CHANGE THE BLOW-UP PARAMETER FROM THE C NOMINAL VALUE OF SQRT(EPS). THE RECIPROCAL OF C THIS PARAMETER IS USED IN REJECTING SOLUTION C COMPONENTS AS TOO LARGE WHEN A VARIABLE IS C FIRST BROUGHT INTO THE ACTIVE SET. TOO LARGE C MEANS THAT THE PROPOSED COMPONENT TIMES THE C RECIPROCAL OF THE PARAMETERIS NOT LESS THAN C THE RATIO OF THE NORMS OF THE RIGHT-SIDE C VECTOR AND THE DATA MATRIX. C THIS PARAMETER CAN BE NO SMALLER THAN EPS, C THE ARITHMETIC-STORAGE PRECISION. C C FOR EXAMPLE, SUPPOSE WE WANT TO PROVIDE C A DIAGONAL MATRIX TO SCALE THE PROBLEM C MATRIX AND CHANGE THE TOLERANCE USED FOR C DETERMINING LINEAR DEPENDENCE OF DROPPED COL C VECTORS. FOR THESE OPTIONS THE DIMENSIONS OF C PRGOPT(*) MUST BE AT LEAST N+6. THE FORTRAN C STATEMENTS DEFINING THESE OPTIONS WOULD C BE AS FOLLOWS. C C PRGOPT(1)=N+3 (LINK TO ENTRY N+3 IN PRGOPT(*)) C PRGOPT(2)=7 (USER-PROVIDED SCALING KEY) C C CALL SCOPY(N,D,1,PRGOPT(3),1) (COPY THE N C SCALING FACTORS FROM A USER ARRAY CALLED D(*) C INTO PRGOPT(3)-PRGOPT(N+2)) C C PRGOPT(N+3)=N+6 (LINK TO ENTRY N+6 OF PRGOPT(*)) C PRGOPT(N+4)=8 (LINEAR DEPENDENCE TOLERANCE KEY) C PRGOPT(N+5)=... (NEW VALUE OF THE TOLERANCE) C C PRGOPT(N+6)=1 (NO MORE OPTIONS TO CHANGE) C C IWORK(1), THE AMOUNTS OF WORKING STORAGE ACTUALLY ALLOCATED C IWORK(2) FOR THE WORKING ARRAYS WORK(*) AND IWORK(*), C RESPECTIVELY. THESE QUANTITIES ARE COMPARED WITH C THE ACTUAL AMOUNTS OF STORAGE NEEDED FOR WNNLS( ). C INSUFFICIENT STORAGE ALLOCATED FOR EITHER WORK(*) C OR IWORK(*) IS CONSIDERED AN ERROR. THIS FEATURE C WAS INCLUDED IN WNNLS( ) BECAUSE MISCALCULATING C THE STORAGE FORMULAS FOR WORK(*) AND IWORK(*) C MIGHT VERY WELL LEAD TO SUBTLE AND HARD-TO-FIND C EXECUTION ERRORS. C C THE LENGTH OF WORK(*) MUST BE AT LEAST C C LW = ME+MA+5*N C THIS TEST WILL NOT BE MADE IF IWORK(1).LE.0. C C THE LENGTH OF IWORK(*) MUST BE AT LEAST C C LIW = ME+MA+N C THIS TEST WILL NOT BE MADE IF IWORK(2).LE.0. C C OUTPUT.. C C X(*) AN ARRAY DIMENSIONED AT LEAST N, WHICH WILL C CONTAIN THE N COMPONENTS OF THE SOLUTION VECTOR C ON OUTPUT. C C RNORM THE RESIDUAL NORM OF THE SOLUTION. THE VALUE OF C RNORM CONTAINS THE RESIDUAL VECTOR LENGTH OF THE C EQUALITY CONSTRAINTS AND LEAST SQUARES EQUATIONS. C C MODE THE VALUE OF MODE INDICATES THE SUCCESS OR FAILURE C OF THE SUBPROGRAM. C C MODE = 0 SUBPROGRAM COMPLETED SUCCESSFULLY. C C = 1 MAX. NUMBER OF ITERATIONS (EQUAL TO C 3*(N-L)) EXCEEDED. NEARLY ALL PROBLEMS C SHOULD COMPLETE IN FEWER THAN THIS C NUMBER OF ITERATIONS. AN APPROXIMATE C SOLUTION AND ITS CORRESPONDING RESIDUAL C VECTOR LENGTH ARE IN X(*) AND RNORM. C C = 2 USAGE ERROR OCCURRED. THE OFFENDING C CONDITION IS NOTED WITH THE ERROR C PROCESSING SUBPROGRAM, XERROR( ). C C USER-DESIGNATED C WORKING ARRAYS.. C C WORK(*) A WORKING ARRAY OF LENGTH AT LEAST C M + 5*N. C C IWORK(*) AN INTEGER-VALUED WORKING ARRAY OF LENGTH AT LEAST C M+N. C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/, DUMMY/,SNGL(DUMMY)/ C C WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES, C AND R.J. HANSON, SANDIA LABORATORIES. C REVISED FEB.25, 1982. C C SUBROUTINES CALLED BY WNNLS( ) C C++ C WNLSM COMPANION SUBROUTINE TO WNNLS( ), WHERE C MOST OF THE COMPUTATION TAKES PLACE. C C XERROR,XERRWV FROM SLATEC ERROR PROCESSING PACKAGE. C THIS IS DOCUMENTED IN SANDIA TECH. REPT., C SAND78-1189. C C REFERENCES C C 1. SOLVING LEAST SQUARES PROBLEMS, BY C.L. LAWSON C AND R.J. HANSON. PRENTICE-HALL, INC. (1974). C C 2. BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, BY C C.L. LAWSON, R.J. HANSON, D.R. KINCAID, AND F.T. KROGH. C TOMS, V. 5, NO. 3, P. 308. ALSO AVAILABLE AS C SANDIA TECHNICAL REPORT NO. SAND77-0898. C C 3. AN ALGORITHM FOR LINEAR LEAST SQUARES WITH EQUALITY C AND NONNEGATIVITY CONSTRAINTS, BY K.H. HASKELL AND C R.J. HANSON. AVAILABLE AS SANDIA TECHNICAL REPORT NO. C SAND77-0552, AND MATH. PROGRAMMING, VOL. 21, (1981), P. 98-118. C C 4. SLATEC COMMON MATH. LIBRARY ERROR HANDLING C PACKAGE. BY R. E. JONES. AVAILABLE AS SANDIA C TECHNICAL REPORT SAND78-1189. C REAL DUMMY, W(MDW,1), PRGOPT(1), X(1), WORK(1), RNORM INTEGER IWORK(1) C C MODE = 0 IF (MA+ME.LE.0 .OR. N.LE.0) RETURN IF (.NOT.(IWORK(1).GT.0)) GO TO 20 LW = ME + MA + 5*N IF (.NOT.(IWORK(1).LT.LW)) GO TO 10 NERR = 2 IOPT = 1 CALL XERRWV(70HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR WORK(* *), NEED LW=I1 BELOW, 70, NERR, IOPT, 1, LW, 0, 0, DUMMY, DUMMY) MODE = 2 RETURN 10 CONTINUE 20 IF (.NOT.(IWORK(2).GT.0)) GO TO 40 LIW = ME + MA + N IF (.NOT.(IWORK(2).LT.LIW)) GO TO 30 NERR = 2 IOPT = 1 CALL XERRWV(72HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR IWORK( **), NEED LIW=I1 BELOW, 72, NERR, IOPT, 1, LIW, 0, 0, DUMMY, DUMMY) MODE = 2 RETURN 30 CONTINUE 40 IF (.NOT.(MDW.LT.ME+MA)) GO TO 50 NERR = 1 IOPT = 1 CALL XERROR(44HWNNLS( ), THE VALUE MDW.LT.ME+MA IS AN ERROR, 44, * NERR, IOPT) MODE = 2 RETURN 50 IF (0.LE.L .AND. L.LE.N) GO TO 60 NERR = 2 IOPT = 1 CALL XERROR(39HWNNLS( ), L.LE.0.AND.L.LE.N IS REQUIRED, 39, NERR, * IOPT) MODE = 2 RETURN C C THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS C WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS C REQUIRED BY THE MAIN SUBROUTINE WNLSM( ). C 60 L1 = N + 1 L2 = L1 + N L3 = L2 + ME + MA L4 = L3 + N L5 = L4 + N C CALL WNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK, * IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3), WORK(L4), * WORK(L5)) RETURN END SUBROUTINE WNLSM(W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, WNL 10 * IPIVOT, ITYPE, WD, H, SCALE, Z, TEMP, D) C C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START CHANGES AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SASUM/DASUM/,/SROTMG/DROTMG/, C /SNRM2/DNRM2/,/ SQRT/ DSQRT/,/SROTM/DROTM/,/AMAX1/DMAX1/, C /SCOPY/DCOPY/,/SSCAL/DSCAL/,/SAXPY/DAXPY/,/E0/D0/,/SSWAP/DSWAP/, C /ISAMAX/IDAMAX/,/SRELPR/DRELPR/,/.E-/.D-/ REMK C C THIS IS A COMPANION SUBPROGRAM TO WNNLS( ). C THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE C USAGE INSTRUCTIONS. C C WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES, C WITH THE HELP OF R.J. HANSON, SANDIA LABORATORIES, C DECEMBER 1976 - JANUARY 1978. C REVISED MAR. 4, 1982. C C IN ADDITION TO THE PARAMETERS DISCUSSED IN THE PROLOGUE TO C SUBROUTINE WNNLS, THE FOLLOWING WORK ARRAYS ARE USED IN C SUBROUTINE WNLSM (THEY ARE PASSED THROUGH THE CALLING C SEQUENCE FROM WNNLS FOR PURPOSES OF VARIABLE DIMENSIONING). C THEIR CONTENTS WILL IN GENERAL BE OF NO INTEREST TO THE USER. C C IPIVOT(*) C AN ARRAY OF LENGTH N. UPON COMPLETION IT CONTAINS THE C PIVOTING INFORMATION FOR THE COLS OF W(*,*). C C ITYPE(*) C AN ARRAY OF LENGTH M WHICH IS USED TO KEEP TRACK C OF THE CLASSIFICATION OF THE EQUATIONS. ITYPE(I)=0 C DENOTES EQUATION I AS AN EQUALITY CONSTRAINT. C ITYPE(I)=1 DENOTES EQUATION I AS A LEAST SQUARES C EQUATION. C C WD(*) C AN ARRAY OF LENGTH N. UPON COMPLETION IT CONTAINS THE C DUAL SOLUTION VECTOR. C C H(*) C AN ARRAY OF LENGTH N. UPON COMPLETION IT CONTAINS THE C PIVOT SCALARS OF THE HOUSEHOLDER TRANSFORMATIONS PERFORMED C IN THE CASE KRANK.LT.L. C C SCALE(*) C AN ARRAY OF LENGTH M WHICH IS USED BY THE SUBROUTINE C TO STORE THE DIAGONAL MATRIX OF WEIGHTS. C THESE ARE USED TO APPLY THE MODIFIED GIVENS C TRANSFORMATIONS. C C Z(*),TEMP(*) C WORKING ARRAYS OF LENGTH N. C C D(*) C AN ARRAY OF LENGTH N THAT CONTAINS THE C COLUMN SCALING FOR THE MATRIX (E). C (A) C C SUBROUTINE WNLSM (W,MDW,MME,MA,N,L,PRGOPT,X,RNORM,MODE, C 1 IPIVOT,ITYPE,WD,H,SCALE,Z,TEMP,D) C++ REAL W(MDW,1), X(1), WD(1), H(1), SCALE(1), DOPE(4) REAL Z(1), TEMP(1), PRGOPT(1), D(1), SPARAM(5) REAL ALAMDA, ALPHA, ALSQ, AMAX, BNORM, EANORM REAL SRELPR, FAC, ONE, BLOWUP REAL RNORM, SM, T, TAU, TWO, WMAX, ZERO, ZZ, Z2 REAL AMAX1, SQRT, SNRM2, SASUM INTEGER IPIVOT(1), ITYPE(1), ISAMAX, IDOPE(8) LOGICAL HITCON, FEASBL, DONE, POS DATA ZERO /0.E0/, ONE /1.E0/, TWO /2.E0/, SRELPR /0.E0/ C C INITIALIZE-VARIABLES ASSIGN 10 TO IGO998 GO TO 180 C C PERFORM INITIAL TRIANGULARIZATION IN THE SUBMATRIX C CORRESPONDING TO THE UNCONSTRAINED VARIABLES USING C THE PROCEDURE INITIALLY-TRIANGULARIZE. 10 ASSIGN 20 TO IGO995 GO TO 280 C C PERFORM WNNLS ALGORITHM USING THE FOLLOWING STEPS. C C UNTIL(DONE) C C COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT C C WHEN (HITCON) ADD-CONSTRAINTS C C ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT C C FIN C C COMPUTE-FINAL-SOLUTION C 20 IF (DONE) GO TO 80 C ASSIGN 30 TO IGO991 GO TO 300 C C COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT C 30 IF (.NOT.(HITCON)) GO TO 50 ASSIGN 40 TO IGO986 GO TO 370 40 GO TO 70 C C WHEN (HITCON) ADD-CONSTRAINTS C 50 ASSIGN 60 TO IGO983 GO TO 640 60 CONTINUE C C ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT C 70 GO TO 20 C 80 ASSIGN 90 TO IGO980 GO TO 1000 C C COMPUTE-FINAL-SOLUTION C 90 RETURN 100 CONTINUE C C TO PROCESS-OPTION-VECTOR FAC = 1.E-4 C C THE NOMINAL TOLERANCE USED IN THE CODE, TAU = SQRT(SRELPR) C C THE NOMINAL BLOW-UP FACTOR USED IN THE CODE. BLOWUP = TAU C C THE NOMINAL COLUMN SCALING USED IN THE CODE IS C THE IDENTITY SCALING. D(1) = ONE CALL SCOPY(N, D, 0, D, 1) C C DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE. NOPT = 1000 C C DEFINE BOUND FOR POSITIVE VALUE OF LINK. NLINK = 100000 NTIMES = 0 LAST = 1 LINK = PRGOPT(1) IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 110 NERR = 3 IOPT = 1 CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR, * IOPT) MODE = 2 RETURN 110 IF (.NOT.(LINK.GT.1)) GO TO 160 NTIMES = NTIMES + 1 IF (.NOT.(NTIMES.GT.NOPT)) GO TO 120 NERR = 3 IOPT = 1 CALL XERROR( * 53HWNNLS( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 53, * NERR, IOPT) MODE = 2 RETURN 120 KEY = PRGOPT(LAST+1) IF (.NOT.(KEY.EQ.6 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 140 DO 130 J=1,N T = SNRM2(M,W(1,J),1) IF (T.NE.ZERO) T = ONE/T D(J) = T 130 CONTINUE 140 IF (KEY.EQ.7) CALL SCOPY(N, PRGOPT(LAST+2), 1, D, 1) IF (KEY.EQ.8) TAU = AMAX1(SRELPR,PRGOPT(LAST+2)) IF (KEY.EQ.9) BLOWUP = AMAX1(SRELPR,PRGOPT(LAST+2)) NEXT = PRGOPT(LINK) IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 150 NERR = 3 IOPT = 1 CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR, * IOPT) MODE = 2 RETURN 150 LAST = LINK LINK = NEXT GO TO 110 160 DO 170 J=1,N CALL SSCAL(M, D(J), W(1,J), 1) 170 CONTINUE GO TO 1260 180 CONTINUE C C TO INITIALIZE-VARIABLES C C SRELPR IS THE PRECISION FOR THE PARTICULAR MACHINE C BEING USED. THIS LOGIC AVOIDS RECOMPUTING IT EVERY ENTRY. IF (.NOT.(SRELPR.EQ.ZERO)) GO TO 210 SRELPR = ONE 190 IF (ONE+SRELPR.EQ.ONE) GO TO 200 SRELPR = SRELPR/TWO GO TO 190 200 SRELPR = SRELPR*TWO 210 M = MA + MME ME = MME MEP1 = ME + 1 ASSIGN 220 TO IGO977 GO TO 100 C C PROCESS-OPTION-VECTOR 220 DONE = .FALSE. ITER = 0 ITMAX = 3*(N-L) MODE = 0 LP1 = L + 1 NSOLN = L NSP1 = NSOLN + 1 NP1 = N + 1 NM1 = N - 1 L1 = MIN0(M,L) C C COMPUTE SCALE FACTOR TO APPLY TO EQUAL. CONSTRAINT EQUAS. DO 230 J=1,N WD(J) = SASUM(M,W(1,J),1) 230 CONTINUE IMAX = ISAMAX(N,WD,1) EANORM = WD(IMAX) BNORM = SASUM(M,W(1,NP1),1) ALAMDA = EANORM/(SRELPR*FAC) C C DEFINE SCALING DIAG MATRIX FOR MOD GIVENS USAGE AND C CLASSIFY EQUATION TYPES. ALSQ = ALAMDA**2 DO 260 I=1,M C C WHEN EQU I IS HEAVILY WEIGHTED ITYPE(I)=0, ELSE ITYPE(I)=1. IF (.NOT.(I.LE.ME)) GO TO 240 T = ALSQ ITEMP = 0 GO TO 250 240 T = ONE ITEMP = 1 250 SCALE(I) = T ITYPE(I) = ITEMP 260 CONTINUE C C SET THE SOLN VECTOR X(*) TO ZERO AND THE COL INTERCHANGE C MATRIX TO THE IDENTITY. X(1) = ZERO CALL SCOPY(N, X, 0, X, 1) DO 270 I=1,N IPIVOT(I) = I 270 CONTINUE GO TO 1230 280 CONTINUE C C TO INITIALLY-TRIANGULARIZE C C SET FIRST L COMPS. OF DUAL VECTOR TO ZERO BECAUSE C THESE CORRESPOND TO THE UNCONSTRAINED VARIABLES. IF (.NOT.(L.GT.0)) GO TO 290 WD(1) = ZERO CALL SCOPY(L, WD, 0, WD, 1) C C THE ARRAYS IDOPE(*) AND DOPE(*) ARE USED TO PASS C INFORMATION TO WNLIT(). THIS WAS DONE TO AVOID C A LONG CALLING SEQUENCE OR THE USE OF COMMON. 290 IDOPE(1) = ME IDOPE(2) = MEP1 IDOPE(3) = 0 IDOPE(4) = 1 IDOPE(5) = NSOLN IDOPE(6) = 0 IDOPE(7) = 1 IDOPE(8) = L1 C DOPE(1) = ALSQ DOPE(2) = EANORM DOPE(3) = FAC DOPE(4) = TAU CALL WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM, * IDOPE, DOPE, DONE) ME = IDOPE(1) MEP1 = IDOPE(2) KRANK = IDOPE(3) KRP1 = IDOPE(4) NSOLN = IDOPE(5) NIV = IDOPE(6) NIV1 = IDOPE(7) L1 = IDOPE(8) GO TO 1240 300 CONTINUE C C TO COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT C C SOLVE THE TRIANGULAR SYSTEM OF CURRENTLY NON-ACTIVE C VARIABLES AND STORE THE SOLUTION IN Z(*). C C SOLVE-SYSTEM ASSIGN 310 TO IGO958 GO TO 1110 C C INCREMENT ITERATION COUNTER AND CHECK AGAINST MAX. NUMBER C OF ITERATIONS. 310 ITER = ITER + 1 IF (.NOT.(ITER.GT.ITMAX)) GO TO 320 MODE = 1 DONE = .TRUE. C C CHECK TO SEE IF ANY CONSTRAINTS HAVE BECOME ACTIVE. C IF SO, CALCULATE AN INTERPOLATION FACTOR SO THAT ALL C ACTIVE CONSTRAINTS ARE REMOVED FROM THE BASIS. 320 ALPHA = TWO HITCON = .FALSE. IF (.NOT.(L.LT.NSOLN)) GO TO 360 DO 350 J=LP1,NSOLN ZZ = Z(J) IF (.NOT.(ZZ.LE.ZERO)) GO TO 340 T = X(J)/(X(J)-ZZ) IF (.NOT.(T.LT.ALPHA)) GO TO 330 ALPHA = T JCON = J 330 HITCON = .TRUE. 340 CONTINUE 350 CONTINUE 360 GO TO 1220 370 CONTINUE C C TO ADD-CONSTRAINTS C C USE COMPUTED ALPHA TO INTERPOLATE BETWEEN LAST C FEASIBLE SOLUTION X(*) AND CURRENT UNCONSTRAINED C (AND INFEASIBLE) SOLUTION Z(*). IF (.NOT.(LP1.LE.NSOLN)) GO TO 390 DO 380 J=LP1,NSOLN X(J) = X(J) + ALPHA*(Z(J)-X(J)) 380 CONTINUE 390 FEASBL = .FALSE. GO TO 410 400 IF (FEASBL) GO TO 610 C C REMOVE COL JCON AND SHIFT COLS JCON+1 THROUGH N TO THE C LEFT. SWAP COL JCON INTO THE N-TH POSITION. THIS ACHIEVES C UPPER HESSENBERG FORM FOR THE NONACTIVE CONSTRAINTS AND C LEAVES AN UPPER HESSENBERG MATRIX TO RETRIANGULARIZE. 410 DO 420 I=1,M T = W(I,JCON) CALL SCOPY(N-JCON, W(I,JCON+1), MDW, W(I,JCON), MDW) W(I,N) = T 420 CONTINUE C C UPDATE PERMUTED INDEX VECTOR TO REFLECT THIS SHIFT AND SWAP. ITEMP = IPIVOT(JCON) IF (.NOT.(JCON.LT.N)) GO TO 440 DO 430 I=JCON,NM1 IPIVOT(I) = IPIVOT(I+1) 430 CONTINUE 440 IPIVOT(N) = ITEMP C C SIMILARLY REPERMUTE X(*) VECTOR. CALL SCOPY(N-JCON, X(JCON+1), 1, X(JCON), 1) X(N) = ZERO NSP1 = NSOLN NSOLN = NSOLN - 1 NIV1 = NIV NIV = NIV - 1 C C RETRIANGULARIZE UPPER HESSENBERG MATRIX AFTER ADDING CONSTRAINTS. J = JCON I = KRANK + JCON - L 450 IF (.NOT.(J.LE.NSOLN)) GO TO 570 IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0)) GO TO 470 ASSIGN 460 TO IGO938 GO TO 620 C C (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0) ZERO-IP1-TO-I-IN-COL-J 460 GO TO 560 470 IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1)) GO TO 490 ASSIGN 480 TO IGO938 GO TO 620 C C (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1) ZERO-IP1-TO-I-IN-COL-J 480 GO TO 560 490 IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.0)) GO TO 510 CALL SSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW) CALL SSWAP(1, SCALE(I), 1, SCALE(I+1), 1) ITEMP = ITYPE(I+1) ITYPE(I+1) = ITYPE(I) ITYPE(I) = ITEMP C C SWAPPED ROW WAS FORMERLY A PIVOT ELT., SO IT WILL C BE LARGE ENOUGH TO PERFORM ELIM. ASSIGN 500 TO IGO938 GO TO 620 C C ZERO-IP1-TO-I-IN-COL-J 500 GO TO 560 510 IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.1)) GO TO 550 T = SCALE(I)*W(I,J)**2/ALSQ IF (.NOT.(T.GT.TAU**2*EANORM**2)) GO TO 530 ASSIGN 520 TO IGO938 GO TO 620 520 GO TO 540 530 CALL SSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW) CALL SSWAP(1, SCALE(I), 1, SCALE(I+1), 1) ITEMP = ITYPE(I+1) ITYPE(I+1) = ITYPE(I) ITYPE(I) = ITEMP W(I+1,J) = ZERO 540 CONTINUE 550 CONTINUE 560 I = I + 1 J = J + 1 GO TO 450 C C SEE IF THE REMAINING COEFFS IN THE SOLN SET ARE FEASIBLE. THEY C SHOULD BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. IF ANY ARE C INFEASIBLE IT IS DUE TO ROUNDOFF ERROR. ANY THAT ARE NON- C POSITIVE WILL BE SET TO ZERO AND REMOVED FROM THE SOLN SET. 570 IF (.NOT.(LP1.LE.NSOLN)) GO TO 590 DO 580 JCON=LP1,NSOLN IF (X(JCON).LE.ZERO) GO TO 600 580 CONTINUE 590 FEASBL = .TRUE. 600 CONTINUE GO TO 400 610 GO TO 1200 620 CONTINUE C C TO ZERO-IP1-TO-I-IN-COL-J IF (.NOT.(W(I+1,J).NE.ZERO)) GO TO 630 CALL SROTMG(SCALE(I), SCALE(I+1), W(I,J), W(I+1,J), SPARAM) W(I+1,J) = ZERO CALL SROTM(NP1-J, W(I,J+1), MDW, W(I+1,J+1), MDW, SPARAM) 630 GO TO 1290 640 CONTINUE C C TO PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT CALL SCOPY(NSOLN, Z, 1, X, 1) IF (.NOT.(NSOLN.LT.N)) GO TO 650 X(NSP1) = ZERO CALL SCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1) 650 I = NIV1 660 IF (.NOT.(I.LE.ME)) GO TO 690 C C RECLASSIFY LEAST SQUARES EQATIONS AS EQUALITIES AS C NECESSARY. IF (.NOT.(ITYPE(I).EQ.0)) GO TO 670 I = I + 1 GO TO 680 670 CALL SSWAP(NP1, W(I,1), MDW, W(ME,1), MDW) CALL SSWAP(1, SCALE(I), 1, SCALE(ME), 1) ITEMP = ITYPE(I) ITYPE(I) = ITYPE(ME) ITYPE(ME) = ITEMP MEP1 = ME ME = ME - 1 680 GO TO 660 C C FORM INNER PRODUCT VECTOR WD(*) OF DUAL COEFFS. 690 IF (.NOT.(NSP1.LE.N)) GO TO 730 DO 720 J=NSP1,N SM = ZERO IF (.NOT.(NSOLN.LT.M)) GO TO 710 DO 700 I=NSP1,M SM = SM + SCALE(I)*W(I,J)*W(I,NP1) 700 CONTINUE 710 WD(J) = SM 720 CONTINUE 730 GO TO 750 740 IF (POS .OR. DONE) GO TO 970 C C FIND J SUCH THAT WD(J)=WMAX IS MAXIMUM. THIS DETERMINES C THAT THE INCOMING COL J WILL REDUCE THE RESIDUAL VECTOR C AND BE POSITIVE. 750 WMAX = ZERO IWMAX = NSP1 IF (.NOT.(NSP1.LE.N)) GO TO 780 DO 770 J=NSP1,N IF (.NOT.(WD(J).GT.WMAX)) GO TO 760 WMAX = WD(J) IWMAX = J 760 CONTINUE 770 CONTINUE 780 IF (.NOT.(WMAX.LE.ZERO)) GO TO 790 DONE = .TRUE. GO TO 960 C C SET DUAL COEFF TO ZERO FOR INCOMING COL. 790 WD(IWMAX) = ZERO C C WMAX .GT. ZERO, SO OKAY TO MOVE COL IWMAX TO SOLN SET. C PERFORM TRANSFORMATION TO RETRIANGULARIZE, AND TEST C FOR NEAR LINEAR DEPENDENCE. C SWAP COL IWMAX INTO NSOLN-TH POSITION TO MAINTAIN UPPER C HESSENBERG FORM OF ADJACENT COLS, AND ADD NEW COL TO C TRIANGULAR DECOMPOSITION. NSOLN = NSP1 NSP1 = NSOLN + 1 NIV = NIV1 NIV1 = NIV + 1 IF (.NOT.(NSOLN.NE.IWMAX)) GO TO 800 CALL SSWAP(M, W(1,NSOLN), 1, W(1,IWMAX), 1) WD(IWMAX) = WD(NSOLN) WD(NSOLN) = ZERO ITEMP = IPIVOT(NSOLN) IPIVOT(NSOLN) = IPIVOT(IWMAX) IPIVOT(IWMAX) = ITEMP C C REDUCE COL NSOLN SO THAT THE MATRIX OF NONACTIVE C CONSTRAINTS VARIABLES IS TRIANGULAR. 800 J = M 810 IF (.NOT.(J.GT.NIV)) GO TO 870 JM1 = J - 1 JP = JM1 C C WHEN OPERATING NEAR THE ME LINE, TEST TO SEE IF THE PIVOT ELT. C IS NEAR ZERO. IF SO, USE THE LARGEST ELT. ABOVE IT AS THE PIVOT. C THIS IS TO MAINTAIN THE SHARP INTERFACE BETWEEN WEIGHTED AND C NON-WEIGHTED ROWS IN ALL CASES. IF (.NOT.(J.EQ.MEP1)) GO TO 850 IMAX = ME AMAX = SCALE(ME)*W(ME,NSOLN)**2 820 IF (.NOT.(JP.GE.NIV)) GO TO 840 T = SCALE(JP)*W(JP,NSOLN)**2 IF (.NOT.(T.GT.AMAX)) GO TO 830 IMAX = JP AMAX = T 830 JP = JP - 1 GO TO 820 840 JP = IMAX 850 IF (.NOT.(W(J,NSOLN).NE.ZERO)) GO TO 860 CALL SROTMG(SCALE(JP), SCALE(J), W(JP,NSOLN), W(J,NSOLN), SPARAM) W(J,NSOLN) = ZERO CALL SROTM(NP1-NSOLN, W(JP,NSP1), MDW, W(J,NSP1), MDW, SPARAM) 860 J = JM1 GO TO 810 C C SOLVE FOR Z(NSOLN)=PROPOSED NEW VALUE FOR X(NSOLN). C TEST IF THIS IS NONPOSITIVE OR TOO LARGE. C IF THIS WAS TRUE OR IF THE PIVOT TERM WAS ZERO REJECT C THE COL AS DEPENDENT. 870 IF (.NOT.(W(NIV,NSOLN).NE.ZERO)) GO TO 890 ISOL = NIV ASSIGN 880 TO IGO897 GO TO 980 C C TEST-PROPOSED-NEW-COMPONENT 880 GO TO 940 890 IF (.NOT.(NIV.LE.ME .AND. W(MEP1,NSOLN).NE.ZERO)) GO TO 920 C C TRY TO ADD ROW MEP1 AS AN ADDITIONAL EQUALITY CONSTRAINT. C CHECK SIZE OF PROPOSED NEW SOLN COMPONENT. C REJECT IT IF IT IS TOO LARGE. ISOL = MEP1 ASSIGN 900 TO IGO897 GO TO 980 C C TEST-PROPOSED-NEW-COMPONENT 900 IF (.NOT.(POS)) GO TO 910 C C SWAP ROWS MEP1 AND NIV, AND SCALE FACTORS FOR THESE ROWS. CALL SSWAP(NP1, W(MEP1,1), MDW, W(NIV,1), MDW) CALL SSWAP(1, SCALE(MEP1), 1, SCALE(NIV), 1) ITEMP = ITYPE(MEP1) ITYPE(MEP1) = ITYPE(NIV) ITYPE(NIV) = ITEMP ME = MEP1 MEP1 = ME + 1 910 GO TO 930 920 POS = .FALSE. 930 CONTINUE 940 IF (POS) GO TO 950 NSP1 = NSOLN NSOLN = NSOLN - 1 NIV1 = NIV NIV = NIV - 1 950 CONTINUE 960 GO TO 740 970 GO TO 1250 980 CONTINUE C C TO TEST-PROPOSED-NEW-COMPONENT Z2 = W(ISOL,NP1)/W(ISOL,NSOLN) Z(NSOLN) = Z2 POS = Z2.GT.ZERO IF (.NOT.(Z2*EANORM.GE.BNORM .AND. POS)) GO TO 990 POS = .NOT.(BLOWUP*Z2*EANORM.GE.BNORM) 990 GO TO 1280 1000 CONTINUE C TO COMPUTE-FINAL-SOLUTION C C SOLVE SYSTEM, STORE RESULTS IN X(*). C ASSIGN 1010 TO IGO958 GO TO 1110 C SOLVE-SYSTEM 1010 CALL SCOPY(NSOLN, Z, 1, X, 1) C C APPLY HOUSEHOLDER TRANSFORMATIONS TO X(*) IF KRANK.LT.L IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.L)) GO TO 1030 DO 1020 I=1,KRANK CALL H12(2, I, KRP1, L, W(I,1), MDW, H(I), X, 1, 1, 1) 1020 CONTINUE C C FILL IN TRAILING ZEROES FOR CONSTRAINED VARIABLES NOT IN SOLN. 1030 IF (.NOT.(NSOLN.LT.N)) GO TO 1040 X(NSP1) = ZERO CALL SCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1) C C REPERMUTE SOLN VECTOR TO NATURAL ORDER. 1040 DO 1070 I=1,N J = I 1050 IF (IPIVOT(J).EQ.I) GO TO 1060 J = J + 1 GO TO 1050 1060 IPIVOT(J) = IPIVOT(I) IPIVOT(I) = J CALL SSWAP(1, X(J), 1, X(I), 1) 1070 CONTINUE C C RESCALE THE SOLN USING THE COL SCALING. DO 1080 J=1,N X(J) = X(J)*D(J) 1080 CONTINUE C IF (.NOT.(NSOLN.LT.M)) GO TO 1100 REMK C DO 1090 I=NSP1,M REMK IF (.NOT.(NIV.LT.M)) GO TO 1100 DO 1090 I = NIV1,M T = W(I,NP1) IF (I.LE.ME) T = T/ALAMDA T = (SCALE(I)*T)*T RNORM = RNORM + T 1090 CONTINUE 1100 RNORM = SQRT(RNORM) GO TO 1210 C C TO SOLVE-SYSTEM C 1110 CONTINUE IF (.NOT.(DONE)) GO TO 1120 ISOL = 1 GO TO 1130 1120 ISOL = LP1 1130 IF (.NOT.(NSOLN.GE.ISOL)) GO TO 1190 C C COPY RT. HAND SIDE INTO TEMP VECTOR TO USE OVERWRITING METHOD. CALL SCOPY(NIV, W(1,NP1), 1, TEMP, 1) DO 1180 JJ=ISOL,NSOLN J = NSOLN - JJ + ISOL IF (.NOT.(J.GT.KRANK)) GO TO 1140 I = NIV - JJ + ISOL GO TO 1150 1140 I = J 1150 IF (.NOT.(J.GT.KRANK .AND. J.LE.L)) GO TO 1160 Z(J) = ZERO GO TO 1170 1160 Z(J) = TEMP(I)/W(I,J) CALL SAXPY(I-1, -Z(J), W(1,J), 1, TEMP, 1) 1170 CONTINUE 1180 CONTINUE 1190 GO TO 1270 1200 GO TO IGO986, (40) 1210 GO TO IGO980, (90) 1220 GO TO IGO991, (30) 1230 GO TO IGO998, (10) 1240 GO TO IGO995, (20) 1250 GO TO IGO983, (60) 1260 GO TO IGO977, (220) 1270 GO TO IGO958, (310, 1010) 1280 GO TO IGO897, (880, 900) 1290 GO TO IGO938, (460, 480, 500, 520) END SUBROUTINE WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM, WNL 10 * IDOPE, DOPE, DONE) C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (BEGIN CHANGES AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SCOPY/DCOPY/,/SROTM/DROTM/, C /SSCAL/DSCAL/,/SQRT/DSQRT/, REMK C /SSWAP/DSWAP/,/AMAX1/DMAX1/,/ISAMAX/IDAMAX/,/.E-/.D-/,/E0/D0/ C C THIS IS A COMPANION SUBPROGRAM TO WNNLS( ). C THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE C USAGE INSTRUCTIONS. C C NOTE THE M BY (N+1) MATRIX W( , ) CONTAINS THE RT. HAND SIDE C B AS THE (N+1)ST COL. C C C TRIANGULARIZE L1 BY L1 SUBSYSTEM, WHERE L1=MIN(M,L), WITH C COL INTERCHANGES. C REVISED MARCH 4, 1982 C C++ REAL W(MDW,1), H(1), SCALE(1), DOPE(4), SPARAM(5) C REAL ALSQ, AMAX, EANORM, FAC, FACTOR, HBAR, ONE, RN REAL ALSQ, EANORM, FACTOR, HBAR, ONE, RN REAL RNORM, SN, T, TAU, TENM3, ZERO REAL AMAX1 INTEGER ITYPE(1), IPIVOT(1), IDOPE(8) INTEGER ISAMAX LOGICAL INDEP, DONE, RECALC DATA TENM3 /1.E-3/, ZERO /0.E0/, ONE /1.E0/ C ME = IDOPE(1) MEP1 = IDOPE(2) KRANK = IDOPE(3) KRP1 = IDOPE(4) NSOLN = IDOPE(5) NIV = IDOPE(6) NIV1 = IDOPE(7) L1 = IDOPE(8) C ALSQ = DOPE(1) EANORM = DOPE(2) C FAC = DOPE(3) REMK TAU = DOPE(4) NP1 = N + 1 LB = MIN0(M-1,L) RECALC = .TRUE. RNORM = ZERO KRANK = 0 C WE SET FACTOR=1.E0 SO THAT THE HEAVY WEIGHT ALAMDA WILL BE C INCLUDED IN THE TEST FOR COL INDEPENDENCE. FACTOR = 1.E0 I = 1 IP1 = 2 LEND = L 10 IF (.NOT.(I.LE.LB)) GO TO 150 IF (.NOT.(I.LE.ME)) GO TO 130 C C SET IR TO POINT TO THE I-TH ROW. IR = I MEND = M ASSIGN 20 TO IGO996 GO TO 460 C C UPDATE-COL-SS-AND-FIND-PIVOT-COL 20 ASSIGN 30 TO IGO993 GO TO 560 C C PERFORM-COL-INTERCHANGE C C SET IC TO POINT TO I-TH COL. 30 IC = I ASSIGN 40 TO IGO990 GO TO 520 C C TEST-INDEP-OF-INCOMING-COL 40 IF (.NOT.(INDEP)) GO TO 110 C C ELIMINATE I-TH COL BELOW DIAG. USING MOD. GIVENS TRANSFORMATIONS C APPLIED TO (A B). J = M DO 100 JJ=IP1,M JM1 = J - 1 JP = JM1 C WHEN OPERATING NEAR THE ME LINE, USE THE LARGEST ELT. REMK C ABOVE IT AS THE PIVOT. REMK C IF (.NOT.(J.EQ.MEP1)) GO TO 80 REMK C IMAX = ME REMK C AMAX = SCALE(ME)*W(ME,I)**2 REMK C 50 IF (.NOT.(JP.GE.I)) GO TO 70 REMK C T = SCALE(JP)*W(JP,I)**2 REMK C IF (.NOT.(T.GT.AMAX)) GO TO 60 REMK C IMAX = JP REMK C AMAX = T REMK C 60 JP = JP - 1 REMK C GO TO 50 REMK C 70 JP = IMAX REMK IF (.NOT. (JJ.EQ.M)) GO TO 70 IF (.NOT. (I.LT.MEP1)) GO TO 80 J = MEP1 JP = I T = SCALE(JP)*W(JP,I)**2*TAU**2 IF (.NOT.(T.GT.SCALE(J)*W(J,I)**2)) GO TO 130 GO TO 80 70 IF (.NOT.(J.EQ.MEP1)) GO TO 80 J = JM1 JM1 = J - 1 JP = JM1 80 IF (.NOT.(W(J,I).NE.ZERO)) GO TO 90 CALL SROTMG(SCALE(JP), SCALE(J), W(JP,I), W(J,I), SPARAM) W(J,I) = ZERO CALL SROTM(NP1-I, W(JP,IP1), MDW, W(J,IP1), MDW, SPARAM) 90 J = JM1 100 CONTINUE GO TO 140 110 CONTINUE IF (.NOT.(LEND.GT.I)) GO TO 130 C C COL I IS DEPENDENT. SWAP WITH COL LEND. MAX = LEND C C PERFORM-COL-INTERCHANGE ASSIGN 120 TO IGO993 GO TO 560 120 CONTINUE LEND = LEND - 1 C C FIND COL IN REMAINING SET WITH LARGEST SS. MAX = ISAMAX(LEND-I+1,H(I),1) + I - 1 HBAR = H(MAX) GO TO 30 130 CONTINUE KRANK = I - 1 GO TO 160 140 I = IP1 IP1 = IP1 + 1 GO TO 10 150 KRANK = L1 160 CONTINUE KRP1 = KRANK + 1 IF (.NOT.(KRANK.LT.ME)) GO TO 290 FACTOR = ALSQ DO 170 I=KRP1,ME IF (L.GT.0) W(I,1) = ZERO CALL SCOPY(L, W(I,1), 0, W(I,1), MDW) 170 CONTINUE C C DETERMINE THE RANK OF THE REMAINING EQUALITY CONSTRAINT C EQUATIONS BY ELIMINATING WITHIN THE BLOCK OF CONSTRAINED C VARIABLES. REMOVE ANY REDUNDANT CONSTRAINTS. IR = KRP1 IF (.NOT. (L.LT.N)) GO TO 245 LP1 = L + 1 RECALC = .TRUE. LB = MIN0(L+ME-KRANK,N) I = LP1 IP1 = I + 1 180 IF (.NOT.(I.LE.LB)) GO TO 280 IR = KRANK + I - L LEND = N MEND = ME ASSIGN 190 TO IGO996 GO TO 460 C C UPDATE-COL-SS-AND-FIND-PIVOT-COL 190 ASSIGN 200 TO IGO993 GO TO 560 C C PERFORM-COL-INTERCHANGE C C ELIMINATE ELEMENTS IN THE I-TH COL. 200 J = ME 210 IF (.NOT.(J.GT.IR)) GO TO 230 JM1 = J - 1 IF (.NOT.(W(J,I).NE.ZERO)) GO TO 220 CALL SROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM) W(J,I) = ZERO CALL SROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM) 220 J = JM1 GO TO 210 C C SET IC=I=COL BEING ELIMINATED 230 IC = I ASSIGN 240 TO IGO990 GO TO 520 C C TEST-INDEP-OF-INCOMING-COL 240 IF (INDEP) GO TO 270 C C REMOVE ANY REDUNDANT OR DEPENDENT EQUALITY CONSTRAINTS. 245 CONTINUE JJ = IR 250 IF (.NOT.(IR.LE.ME)) GO TO 260 W(IR,1) = ZERO CALL SCOPY(N, W(IR,1), 0, W(IR,1), MDW) RNORM = RNORM + (SCALE(IR)*W(IR,NP1)/ALSQ)*W(IR,NP1) W(IR,NP1) = ZERO SCALE(IR) = ONE C RECLASSIFY THE ZEROED ROW AS A LEAST SQUARES EQUATION. ITYPE(IR) = 1 IR = IR + 1 GO TO 250 C C REDUCE ME TO REFLECT ANY DISCOVERED DEPENDENT EQUALITY C CONSTRAINTS. 260 CONTINUE ME = JJ - 1 MEP1 = ME + 1 GO TO 300 270 I = IP1 IP1 = IP1 + 1 GO TO 180 280 CONTINUE 290 CONTINUE 300 CONTINUE IF (.NOT.(KRANK.LT.L1)) GO TO 420 C C TRY TO DETERMINE THE VARIABLES KRANK+1 THROUGH L1 FROM THE C LEAST SQUARES EQUATIONS. CONTINUE THE TRIANGULARIZATION WITH C PIVOT ELEMENT W(MEP1,I). C RECALC = .TRUE. C C SET FACTOR=ALSQ TO REMOVE EFFECT OF HEAVY WEIGHT FROM C TEST FOR COL INDEPENDENCE. FACTOR = ALSQ KK = KRP1 I = KK IP1 = I + 1 310 IF (.NOT.(I.LE.L1)) GO TO 410 C C SET IR TO POINT TO THE MEP1-ST ROW. IR = MEP1 LEND = L MEND = M ASSIGN 320 TO IGO996 GO TO 460 C C UPDATE-COL-SS-AND-FIND-PIVOT-COL 320 ASSIGN 330 TO IGO993 GO TO 560 C C PERFORM-COL-INTERCHANGE C C ELIMINATE I-TH COL BELOW THE IR-TH ELEMENT. 330 IRP1 = IR + 1 IF (.NOT.(IRP1.LE.M)) GO TO 355 J = M DO 350 JJ=IRP1,M JM1 = J - 1 IF (.NOT.(W(J,I).NE.ZERO)) GO TO 340 CALL SROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM) W(J,I) = ZERO CALL SROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM) 340 J = JM1 350 CONTINUE 355 CONTINUE C C TEST IF NEW PIVOT ELEMENT IS NEAR ZERO. IF SO, THE COL IS C DEPENDENT. T = SCALE(IR)*W(IR,I)**2 INDEP = T.GT.TAU**2*EANORM**2 IF (.NOT.INDEP) GO TO 380 C C COL TEST PASSED. NOW MUST PASS ROW NORM TEST TO BE CLASSIFIED C AS INDEPENDENT. RN = ZERO DO 370 I1=IR,M DO 360 J1=IP1,N RN = AMAX1(RN,SCALE(I1)*W(I1,J1)**2) 360 CONTINUE 370 CONTINUE INDEP = T.GT.TAU**2*RN C C IF INDEPENDENT, SWAP THE IR-TH AND KRP1-ST ROWS TO MAINTAIN THE C TRIANGULAR FORM. UPDATE THE RANK INDICATOR KRANK AND THE C EQUALITY CONSTRAINT POINTER ME. 380 IF (.NOT.(INDEP)) GO TO 390 CALL SSWAP(NP1, W(KRP1,1), MDW, W(IR,1), MDW) CALL SSWAP(1, SCALE(KRP1), 1, SCALE(IR), 1) C RECLASSIFY THE LEAST SQ. EQUATION AS AN EQUALITY CONSTRAINT AND C RESCALE IT. ITYPE(IR) = 0 T = SQRT(SCALE(KRP1)) CALL SSCAL(NP1, T, W(KRP1,1), MDW) SCALE(KRP1) = ALSQ ME = MEP1 MEP1 = ME + 1 KRANK = KRP1 KRP1 = KRANK + 1 GO TO 400 390 GO TO 430 400 I = IP1 IP1 = IP1 + 1 GO TO 310 410 CONTINUE 420 CONTINUE 430 CONTINUE C C IF PSEUDORANK IS LESS THAN L, APPLY HOUSEHOLDER TRANS. C FROM RIGHT. IF (.NOT.(KRANK.LT.L)) GO TO 450 DO 440 I=1,KRANK J = KRP1 - I CALL H12(1, J, KRP1, L, W(J,1), MDW, H(J), W, MDW, 1, J-1) 440 CONTINUE 450 NIV = KRANK + NSOLN - L NIV1 = NIV + 1 IF (L.EQ.N) DONE = .TRUE. C C END OF INITIAL TRIANGULARIZATION. IDOPE(1) = ME IDOPE(2) = MEP1 IDOPE(3) = KRANK IDOPE(4) = KRP1 IDOPE(5) = NSOLN IDOPE(6) = NIV IDOPE(7) = NIV1 IDOPE(8) = L1 RETURN 460 CONTINUE C C TO UPDATE-COL-SS-AND-FIND-PIVOT-COL C C THE COL SS VECTOR WILL BE UPDATED AT EACH STEP. WHEN C NUMERICALLY NECESSARY, THESE VALUES WILL BE RECOMPUTED. C IF (.NOT.(IR.NE.1 .AND. (.NOT.RECALC))) GO TO 480 C UPDATE COL SS =SUM OF SQUARES. DO 470 J=I,LEND H(J) = H(J) - SCALE(IR-1)*W(IR-1,J)**2 470 CONTINUE C C TEST FOR NUMERICAL ACCURACY. MAX = ISAMAX(LEND-I+1,H(I),1) + I - 1 RECALC = HBAR + TENM3*H(MAX).EQ.HBAR C C IF REQUIRED, RECALCULATE COL SS, USING ROWS IR THROUGH MEND. 480 IF (.NOT.(RECALC)) GO TO 510 DO 500 J=I,LEND H(J) = ZERO DO 490 K=IR,MEND H(J) = H(J) + SCALE(K)*W(K,J)**2 490 CONTINUE 500 CONTINUE C C FIND COL WITH LARGEST SS. MAX = ISAMAX(LEND-I+1,H(I),1) + I - 1 HBAR = H(MAX) 510 GO TO 600 520 CONTINUE C C TO TEST-INDEP-OF-INCOMING-COL C C TEST THE COL IC TO DETERMINE IF IT IS LINEARLY INDEPENDENT C OF THE COLS ALREADY IN THE BASIS. IN THE INIT TRI C STEP, WE USUALLY WANT THE HEAVY WEIGHT ALAMDA TO C BE INCLUDED IN THE TEST FOR INDEPENDENCE. IN THIS CASE THE C VALUE OF FACTOR WILL HAVE BEEN SET TO 1.E0 BEFORE THIS C PROCEDURE IS INVOKED. IN THE POTENTIALLY RANK DEFICIENT C PROBLEM, THE VALUE OF FACTOR WILL HAVE BEEN C SET TO ALSQ=ALAMDA**2 TO REMOVE THE EFFECT OF THE HEAVY WEIGHT C FROM THE TEST FOR INDEPENDENCE. C C WRITE NEW COL AS PARTITIONED VECTOR C (A1) NUMBER OF COMPONENTS IN SOLN SO FAR = NIV C (A2) M-NIV COMPONENTS C AND COMPUTE SN = INVERSE WEIGHTED LENGTH OF A1 C RN = INVERSE WEIGHTED LENGTH OF A2 C CALL THE COL INDEPENDENT WHEN RN .GT. TAU*SN SN = ZERO RN = ZERO DO 550 J=1,MEND T = SCALE(J) IF (J.LE.ME) T = T/FACTOR T = T*W(J,IC)**2 IF (.NOT.(J.LT.IR)) GO TO 530 SN = SN + T GO TO 540 530 RN = RN + T 540 CONTINUE 550 CONTINUE INDEP = RN.GT.TAU**2*SN GO TO 590 560 CONTINUE C C TO PERFORM-COL-INTERCHANGE C IF (.NOT.(MAX.NE.I)) GO TO 570 C EXCHANGE ELEMENTS OF PERMUTED INDEX VECTOR AND PERFORM COL C INTERCHANGES. ITEMP = IPIVOT(I) IPIVOT(I) = IPIVOT(MAX) IPIVOT(MAX) = ITEMP CALL SSWAP(M, W(1,MAX), 1, W(1,I), 1) T = H(MAX) H(MAX) = H(I) H(I) = T 570 GO TO 580 580 GO TO IGO993, (30, 200, 330, 120) 590 GO TO IGO990, (40, 240) 600 GO TO IGO996, (20, 190, 320) END SUBROUTINE HFTI(A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, HFTI 10 * G, IP) C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (BEGIN CHANGES AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/ SQRT/ DSQRT/, C /, ABS/, DABS/,/ABS(/DABS(/,/E0/D0/ C C DIMENSION A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N) C C WRITTEN BY C. L. LAWSON AND R. J. HANSON. FROM THE BOOK SOLVING C LEAST SQUARES PROBLEMS, PRENTICE-HALL, INC. (1974). FOR FURTHER C ALGORITHMIC DETAILS SEE ALGORITHM HFTI IN CHAPTER 14. C C ABSTRACT C C THIS SUBROUTINE SOLVES A LINEAR LEAST SQUARES PROBLEM OR A SET OF C LINEAR LEAST SQUARES PROBLEMS HAVING THE SAME MATRIX BUT DIFFERENT C RIGHT-SIDE VECTORS. THE PROBLEM DATA CONSISTS OF AN M BY N MATRIX C A, AN M BY NB MATRIX B, AND AN ABSOLUTE TOLERANCE PARAMETER TAU C WHOSE USAGE IS DESCRIBED BELOW. THE NB COLUMN VECTORS OF B C REPRESENT RIGHT-SIDE VECTORS FOR NB DISTINCT LINEAR LEAST SQUARES C PROBLEMS. C C THIS SET OF PROBLEMS CAN ALSO BE WRITTEN AS THE MATRIX LEAST C SQUARES PROBLEM C C AX = B, C C WHERE X IS THE N BY NB SOLUTION MATRIX. C C NOTE THAT IF B IS THE M BY M IDENTITY MATRIX, THEN X WILL BE THE C PSEUDO-INVERSE OF A. C C THIS SUBROUTINE FIRST TRANSFORMS THE AUGMENTED MATRIX (A B) TO A C MATRIX (R C) USING PREMULTIPLYING HOUSEHOLDER TRANSFORMATIONS WITH C COLUMN INTERCHANGES. ALL SUBDIAGONAL ELEMENTS IN THE MATRIX R ARE C ZERO AND ITS DIAGONAL ELEMENTS SATISFY C C ABS(R(I,I)).GE.ABS(R(I+1,I+1)), C C I = 1,...,L-1, WHERE C C L = MIN(M,N). C C THE SUBROUTINE WILL COMPUTE AN INTEGER, KRANK, EQUAL TO THE NUMBER C OF DIAGONAL TERMS OF R THAT EXCEED TAU IN MAGNITUDE. THEN A C SOLUTION OF MINIMUM EUCLIDEAN LENGTH IS COMPUTED USING THE FIRST C KRANK ROWS OF (R C). C C TO BE SPECIFIC WE SUGGEST THAT THE USER CONSIDER AN EASILY C COMPUTABLE MATRIX NORM, SUCH AS, THE MAXIMUM OF ALL COLUMN SUMS OF C MAGNITUDES. C C NOW IF THE RELATIVE UNCERTAINTY OF B IS EPS, (NORM OF UNCERTAINTY/ C NORM OF B), IT IS SUGGESTED THAT TAU BE SET APPROXIMATELY EQUAL TO C EPS*(NORM OF A). C C THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST.. C A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N). THIS C PERMITS THE SOLUTION OF A RANGE OF PROBLEMS IN THE SAME ARRAY C SPACE. C C THE ENTIRE SET OF PARAMETERS FOR HFTI ARE C C INPUT.. C C A(*,*),MDA,M,N THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY N C MATRIX A OF THE LEAST SQUARES PROBLEM AX = B. C THE FIRST DIMENSIONING PARAMETER OF THE ARRAY C A(*,*) IS MDA, WHICH MUST SATISFY MDA.GE.M C EITHER M.GE.N OR M.LT.N IS PERMITTED. THERE C IS NO RESTRICTION ON THE RANK OF A. THE C CONDITION MDA.LT.M IS CONSIDERED AN ERROR. C C B(*),MDB,NB IF NB = 0 THE SUBROUTINE WILL PERFORM THE C ORTHOGONAL DECOMPOSITION BUT WILL MAKE NO C REFERENCES TO THE ARRAY B(*). IF NB.GT.0 C THE ARRAY B(*) MUST INITIALLY CONTAIN THE M BY C NB MATRIX B OF THE LEAST SQUARES PROBLEM AX = C B. IF NB.GE.2 THE ARRAY B(*) MUST BE DOUBLY C SUBSCRIPTED WITH FIRST DIMENSIONING PARAMETER C MDB.GE.MAX(M,N). IF NB = 1 THE ARRAY B(*) MAY C BE EITHER DOUBLY OR SINGLY SUBSCRIPTED. IN C THE LATTER CASE THE VALUE OF MDB IS ARBITRARY C BUT IT SHOULD BE SET TO SOME VALID INTEGER C VALUE SUCH AS MDB = M. C C THE CONDITION OF NB.GT.1.AND.MDB.LT. MAX(M,N) C IS CONSIDERED AN ERROR. C C TAU ABSOLUTE TOLERANCE PARAMETER PROVIDED BY USER C FOR PSEUDORANK DETERMINATION. C C H(*),G(*),IP(*) ARRAYS OF WORKING SPACE USED BY HFTI. C C OUTPUT.. C C A(*,*) THE CONTENTS OF THE ARRAY A(*,*) WILL BE C MODIFIED BY THE SUBROUTINE. THESE CONTENTS C ARE NOT GENERALLY REQUIRED BY THE USER. C C B(*) ON RETURN THE ARRAY B(*) WILL CONTAIN THE N BY C NB SOLUTION MATRIX X. C C KRANK SET BY THE SUBROUTINE TO INDICATE THE C PSEUDORANK OF A. C C RNORM(*) ON RETURN, RNORM(J) WILL CONTAIN THE EUCLIDEAN C NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM C DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY C B(*,*) FOR J = 1,...,NB. C C H(*),G(*) ON RETURN THESE ARRAYS RESPECTIVELY CONTAIN C ELEMENTS OF THE PRE- AND POST-MULTIPLYING C HOUSEHOLDER TRANSFORMATIONS USED TO COMPUTE C THE MINIMUM EUCLIDEAN LENGTH SOLUTION. C C IP(*) ARRAY IN WHICH THE SUBROUTINE RECORDS INDICES C DESCRIBING THE PERMUTATION OF COLUMN VECTORS. C THE CONTENTS OF ARRAYS H(*),G(*) AND IP(*) C ARE NOT GENERALLY REQUIRED BY THE USER. C C++ REAL A(MDA,N), B(MDB,1), H(N), G(N), RNORM(NB), TAU REAL FACTOR, HMAX, SM1, ZERO, SM, TMP REAL DIFF, SQRT, ABS INTEGER IP(N) ZERO = 0.E0 FACTOR = 0.001E0 C K = 0 LDIAG = MIN0(M,N) IF (LDIAG.LE.0) GO TO 310 IF (.NOT.MDA.LT.M) GO TO 10 NERR = 2 IOPT = 2 CALL XERROR(31HHFTI MDA.LT.M.. PROBABLE ERROR., 31, NERR, IOPT) RETURN 10 CONTINUE C IF (.NOT.(NB.GT.1 .AND. MAX0(M,N).GT.MDB)) GO TO 20 NERR = 2 IOPT = 2 CALL XERROR(49HHFTI MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR., * 49, NERR, IOPT) RETURN 20 CONTINUE C DO 100 J=1,LDIAG IF (J.EQ.1) GO TO 40 C C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX C .. LMAX = J DO 30 L=J,N H(L) = H(L) - A(J-1,L)**2 IF (H(L).GT.H(LMAX)) LMAX = L 30 CONTINUE IF (DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 40, 40, 70 C C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX C .. 40 LMAX = J DO 60 L=J,N H(L) = ZERO DO 50 I=J,M H(L) = H(L) + A(I,L)**2 50 CONTINUE IF (H(L).GT.H(LMAX)) LMAX = L 60 CONTINUE HMAX = H(LMAX) C .. C LMAX HAS BEEN DETERMINED C C DO COLUMN INTERCHANGES IF NEEDED. C .. 70 CONTINUE IP(J) = LMAX IF (IP(J).EQ.J) GO TO 90 DO 80 I=1,M TMP = A(I,J) A(I,J) = A(I,LMAX) A(I,LMAX) = TMP 80 CONTINUE H(LMAX) = H(J) 90 JCOL = MIN0(J+1,N) C C COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B. C .. CALL H12(1, J, J+1, M, A(1,J), 1, H(J), A(1,JCOL), 1, MDA, N-J) CALL H12(2, J, J+1, M, A(1,J), 1, H(J), B, 1, MDB, NB) 100 CONTINUE C C DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU. C .. DO 110 J=1,LDIAG IF (ABS(A(J,J)).LE.TAU) GO TO 120 110 CONTINUE K = LDIAG GO TO 130 120 K = J - 1 130 KP1 = K + 1 C C COMPUTE THE NORMS OF THE RESIDUAL VECTORS. C IF (NB.LE.0) GO TO 170 DO 160 JB=1,NB TMP = ZERO IF (KP1.GT.M) GO TO 150 DO 140 I=KP1,M TMP = TMP + B(I,JB)**2 140 CONTINUE 150 RNORM(JB) = SQRT(TMP) 160 CONTINUE 170 CONTINUE C SPECIAL FOR PSEUDORANK = 0 IF (K.GT.0) GO TO 200 IF (NB.LE.0) GO TO 310 DO 190 JB=1,NB DO 180 I=1,N B(I,JB) = ZERO 180 CONTINUE 190 CONTINUE GO TO 310 C C IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER C DECOMPOSITION OF FIRST K ROWS. C .. 200 IF (K.EQ.N) GO TO 220 DO 210 II=1,K I = KP1 - II CALL H12(1, I, KP1, N, A(I,1), MDA, G(I), A, MDA, 1, I-1) 210 CONTINUE 220 CONTINUE C C IF (NB.LE.0) GO TO 310 DO 300 JB=1,NB C C SOLVE THE K BY K TRIANGULAR SYSTEM. C .. DO 250 L=1,K SM = ZERO I = KP1 - L IF (I.EQ.K) GO TO 240 IP1 = I + 1 DO 230 J=IP1,K SM = SM + A(I,J)*B(J,JB) 230 CONTINUE 240 SM1 = SM B(I,JB) = (B(I,JB)-SM1)/A(I,I) 250 CONTINUE C C COMPLETE COMPUTATION OF SOLUTION VECTOR. C .. IF (K.EQ.N) GO TO 280 DO 260 J=KP1,N B(J,JB) = ZERO 260 CONTINUE DO 270 I=1,K CALL H12(2, I, KP1, N, A(I,1), MDA, G(I), B(1,JB), 1, MDB, 1) 270 CONTINUE C C RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE C COLUMN INTERCHANGES. C .. 280 DO 290 JJ=1,LDIAG J = LDIAG + 1 - JJ IF (IP(J).EQ.J) GO TO 290 L = IP(J) TMP = B(L,JB) B(L,JB) = B(J,JB) B(J,JB) = TMP 290 CONTINUE 300 CONTINUE C .. C THE SOLUTION VECTORS, X, ARE NOW C IN THE FIRST N ROWS OF THE ARRAY B(,). C 310 KRANK = K RETURN END SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) H 12 10 C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (START CHANGES AT LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SDOT/DDOT/,/ABS,/DABS,/, C /SSWAP/DSWAP/,/SQRT/DSQRT/,/ABS(/ DABS(/,/AMAX1/DMAX1/, C /SAXPY/DAXPY/,/E0/D0/ C C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 C C MODIFIED AT SANDIA LABS., MAY 1977, TO -- C C 1) REMOVE DOUBLE PRECISION ACCUMULATION, AND C 2) INCLUDE USAGE OF THE BASIC LINEAR ALGEBRA PACKAGE FOR C VECTORS LONGER THAN A PARTICULAR THRESHOLD. C C CONSTRUCTION AND/OR APPLICATION OF A SINGLE C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B C C MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 . C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. C U(),IUE,UP ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR. C IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS. C ON EXIT FROM H1 U() AND UP C CONTAIN QUANTITIES DEFINING THE VECTOR U OF THE C HOUSEHOLDER TRANSFORMATION. ON ENTRY TO H2 U() C AND UP SHOULD CONTAIN QUANTITIES PREVIOUSLY COMPUTED C BY H1. THESE WILL NOT BE MODIFIED BY H2. C C() ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE C REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER C TRANSFORMATION IS TO BE APPLIED. ON EXIT C() CONTAINS THE C SET OF TRANSFORMED VECTORS. C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). C ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 C NO OPERATIONS WILL BE DONE ON C(). C C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) C++ REAL U(IUE,M), C(1), UP REAL B, CL, CLINV, ONE, SM, UL1M1 REAL ABS, AMAX1, SQRT, SDOT ONE=1.E0 C IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN CL=ABS(U(1,LPIVOT)) IF (MODE.EQ.2) GO TO 60 C ****** CONSTRUCT THE TRANSFORMATION. ****** DO 10 J=L1,M 10 CL=AMAX1(ABS(U(1,J)),CL) IF (CL) 130,130,20 20 CLINV=ONE/CL SM=(U(1,LPIVOT)*CLINV)**2 DO 30 J=L1,M 30 SM=SM+(U(1,J)*CLINV)**2 CL=CL*SQRT(SM) IF (U(1,LPIVOT)) 50,50,40 40 CL=-CL 50 UP=U(1,LPIVOT)-CL U(1,LPIVOT)=CL GO TO 70 C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** C 60 IF (CL) 130,130,70 70 IF (NCV.LE.0) RETURN B=UP*U(1,LPIVOT) C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. C IF (B) 80,130,130 80 B=ONE/B MML1P2=M-L1+2 IF (MML1P2.GT.20) GO TO 140 I2=1-ICV+ICE*(LPIVOT-1) INCR=ICE*(L1-LPIVOT) DO 120 J=1,NCV I2=I2+ICV I3=I2+INCR I4=I3 SM=C(I2)*UP DO 90 I=L1,M SM=SM+C(I3)*U(1,I) 90 I3=I3+ICE IF (SM) 100,120,100 100 SM=SM*B C(I2)=C(I2)+SM*UP DO 110 I=L1,M C(I4)=C(I4)+SM*U(1,I) 110 I4=I4+ICE 120 CONTINUE 130 RETURN 140 CONTINUE L1M1=L1-1 KL1=1+(L1M1-1)*ICE KL2=KL1 KLP=1+(LPIVOT-1)*ICE UL1M1=U(1,L1M1) U(1,L1M1)=UP IF (LPIVOT.EQ.L1M1) GO TO 150 CALL SSWAP(NCV,C(KL1),ICV,C(KLP),ICV) 150 CONTINUE DO 160 J=1,NCV SM=SDOT(MML1P2,U(1,L1M1),IUE,C(KL1),ICE) SM=SM*B CALL SAXPY (MML1P2,SM,U(1,L1M1),IUE,C(KL1),ICE) KL1=KL1+ICV 160 CONTINUE U(1,L1M1)=UL1M1 IF (LPIVOT.EQ.L1M1) RETURN KL1=KL2 CALL SSWAP(NCV,C(KL1),ICV,C(KLP),ICV) RETURN END REAL FUNCTION DIFF(X,Y) DIFF 10 C C THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO C DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. C USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. C (APPLY CHANGES TO ENTIRE PROGRAM UNIT.) C /REAL (12 BLANKS)/DOUBLE PRECISION/ C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 REAL X, Y DIFF=X-Y RETURN END SUBROUTINE SROTMG (SD1,SD2,SX1,SY1,SPARAM) SBLA 10 C C CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS C THE SECOND COMPONENT OF THE 2-VECTOR (SQRT(SD1)*SX1,SQRT(SD2)* C SY2)**T. C WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. C C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0 C C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0) C H=( ) ( ) ( ) ( ) C (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0). C C END OF ABSTRACT DIMENSION SPARAM(5) C DATA ZERO,ONE,TWO /0.E0,1.E0,2.E0/,IFLAG/1/ DATA GAM,GAMSQ,RGAM,RGAMSQ/4096.E0,1.678E7,2.441E-4,5.960E-8/ C IF(.NOT. SD1 .LT. ZERO) GO TO 10 C GO ZERO-H-D-AND-SX1.. GO TO 60 10 CONTINUE C CASE-SD1-NONNEGATIVE SP2=SD2*SY1 IF(.NOT. SP2 .EQ. ZERO) GO TO 20 SFLAG=-TWO GO TO 260 C REGULAR-CASE.. 20 CONTINUE SP1=SD1*SX1 SQ2=SP2*SY1 SQ1=SP1*SX1 C IF(.NOT. ABS(SQ1) .GT. ABS(SQ2)) GO TO 40 SH21=-SY1/SX1 SH12=SP2/SP1 C SU=ONE-SH12*SH21 C IF(.NOT. SU .LE. ZERO) GO TO 30 C GO ZERO-H-D-AND-SX1.. GO TO 60 30 CONTINUE SFLAG=ZERO SD1=SD1/SU SD2=SD2/SU SX1=SX1*SU C GO SCALE-CHECK.. GO TO 100 40 CONTINUE IF(.NOT. SQ2 .LT. ZERO) GO TO 50 C GO ZERO-H-D-AND-SX1.. GO TO 60 50 CONTINUE SFLAG=ONE SH11=SP1/SP2 SH22=SX1/SY1 SU=ONE+SH11*SH22 STEMP=SD2/SU SD2=SD1/SU SD1=STEMP SX1=SY1*SU C GO SCALE-CHECK GO TO 100 C PROCEDURE..ZERO-H-D-AND-SX1.. 60 CONTINUE SFLAG=-ONE SH11=ZERO SH12=ZERO SH21=ZERO SH22=ZERO C SD1=ZERO SD2=ZERO SX1=ZERO C RETURN.. GO TO 220 C PROCEDURE..FIX-H.. 70 CONTINUE IF(.NOT. SFLAG .GE. ZERO) GO TO 90 C IF(.NOT. SFLAG .EQ. ZERO) GO TO 80 SH11=ONE SH22=ONE SFLAG=-ONE GO TO 90 80 CONTINUE SH21=-ONE SH12=ONE SFLAG=-ONE 90 CONTINUE GO TO IGO,(120,150,180,210) C PROCEDURE..SCALE-CHECK 100 CONTINUE IF(.NOT. IFLAG.EQ.1) GO TO 105 C C RECOMPUTE RESCALING PARAMETERS C MORE ACCURATELY.. C RGAM = ONE/GAM GAMSQ = GAM**2 RGAMSQ = RGAM**2 IFLAG = 2 105 CONTINUE 110 CONTINUE IF(.NOT. SD1 .LE. RGAMSQ) GO TO 130 IF(SD1 .EQ. ZERO) GO TO 160 ASSIGN 120 TO IGO C FIX-H.. GO TO 70 120 CONTINUE SD1=SD1*GAMSQ SX1=SX1*RGAM SH11=SH11*RGAM SH12=SH12*RGAM GO TO 110 130 CONTINUE 140 CONTINUE IF(.NOT. SD1 .GE. GAMSQ) GO TO 160 ASSIGN 150 TO IGO C FIX-H.. GO TO 70 150 CONTINUE SD1=SD1*RGAMSQ SX1=SX1*GAM SH11=SH11*GAM SH12=SH12*GAM GO TO 140 160 CONTINUE 170 CONTINUE IF(.NOT. ABS(SD2) .LE. RGAMSQ) GO TO 190 IF(SD2 .EQ. ZERO) GO TO 220 ASSIGN 180 TO IGO C FIX-H.. GO TO 70 180 CONTINUE SD2=SD2*GAMSQ SH21=SH21*RGAM SH22=SH22*RGAM GO TO 170 190 CONTINUE 200 CONTINUE IF(.NOT. ABS(SD2) .GE. GAMSQ) GO TO 220 ASSIGN 210 TO IGO C FIX-H.. GO TO 70 210 CONTINUE SD2=SD2*RGAMSQ SH21=SH21*GAM SH22=SH22*GAM GO TO 200 220 CONTINUE IF(SFLAG)250,230,240 230 CONTINUE SPARAM(3)=SH21 SPARAM(4)=SH12 GO TO 260 240 CONTINUE SPARAM(2)=SH11 SPARAM(5)=SH22 GO TO 260 250 CONTINUE SPARAM(2)=SH11 SPARAM(3)=SH21 SPARAM(4)=SH12 SPARAM(5)=SH22 260 CONTINUE SPARAM(1)=SFLAG RETURN END SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) SBLA 10 C C COPY SINGLE PRECISION SX TO SINGLE PRECISION SY. C C END OF ABSTRACT REAL SX(1),SY(1) C IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 20 M = N - (N/7)*7 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 SY(I) = SX(I) SY(I + 1) = SX(I + 1) SY(I + 2) = SX(I + 2) SY(I + 3) = SX(I + 3) SY(I + 4) = SX(I + 4) SY(I + 5) = SX(I + 5) SY(I + 6) = SX(I + 6) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SX(I) 70 CONTINUE RETURN END SUBROUTINE SSWAP (N,SX,INCX,SY,INCY) SBLA 10 C C INTERCHANGE SINGLE PRECISION SX AND SINGLE PRECISION SY. C C END OF ABSTRACT REAL SX(1),SY(1),STEMP1,STEMP2,STEMP3 C IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N STEMP1 = SX(IX) SX(IX) = SY(IY) SY(IY) = STEMP1 IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3. C 20 M = N - (N/3)*3 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M STEMP1 = SX(I) SX(I) = SY(I) SY(I) = STEMP1 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 STEMP1 = SX(I) STEMP2 = SX(I+1) STEMP3 = SX(I+2) SX(I) = SY(I) SX(I+1) = SY(I+1) SX(I+2) = SY(I+2) SY(I) = STEMP1 SY(I+1) = STEMP2 SY(I+2) = STEMP3 50 CONTINUE RETURN 60 CONTINUE C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C NS = N*INCX DO 70 I=1,NS,INCX STEMP1 = SX(I) SX(I) = SY(I) SY(I) = STEMP1 70 CONTINUE RETURN END REAL FUNCTION SNRM2 ( N, SX, INCX) SBLA 10 INTEGER NEXT REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE DATA ZERO, ONE /0.0E0, 1.0E0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C END OF ABSTRACT DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C IF(N .GT. 0) GO TO 10 SNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( SX(I) .EQ. ZERO) GO TO 200 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / SX(I)) / SX(I) 105 XMAX = ABS(SX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( ABS(SX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / SX(I))**2 XMAX = ABS(SX(I)) GO TO 200 C 115 SUM = SUM + (SX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(ABS(SX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + SX(J)**2 SNRM2 = SQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C SNRM2 = XMAX * SQRT(SUM) 300 CONTINUE RETURN END REAL FUNCTION SASUM(N,SX,INCX) SBLA 10 C C RETURNS SUM OF MAGNITUDES OF SINGLE PRECISION SX. C C END OF ABSTRACT REAL SX(1) C SASUM = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I=1,NS,INCX SASUM = SASUM + ABS(SX(I)) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6. C 20 M = N - (N/6)*6 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SASUM = SASUM + ABS(SX(I)) 30 CONTINUE IF( N .LT. 6 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 SASUM = SASUM + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2)) A + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5)) 50 CONTINUE RETURN END SUBROUTINE SSCAL(N,SA,SX,INCX) SBLA 10 C C REPLACE SINGLE PRECISION SX BY SINGLE PRECISION SA*SX. C C END OF ABSTRACT REAL SA,SX(1) C IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I = 1,NS,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = N - (N/5)*5 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SX(I) = SA*SX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SX(I) = SA*SX(I) SX(I + 1) = SA*SX(I + 1) SX(I + 2) = SA*SX(I + 2) SX(I + 3) = SA*SX(I + 3) SX(I + 4) = SA*SX(I + 4) 50 CONTINUE RETURN END INTEGER FUNCTION ISAMAX(N,SX,INCX) SBLA 10 C C FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF SINGLE PRECISION SX. C C END OF ABSTRACT REAL SX(1),SMAX,XMAG C ISAMAX = 0 IF(N.LE.0) RETURN ISAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C SMAX = ABS(SX(1)) NS = N*INCX II = 1 DO 10 I=1,NS,INCX XMAG = ABS(SX(I)) IF(XMAG.LE.SMAX) GO TO 5 ISAMAX = II SMAX = XMAG 5 II = II + 1 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C 20 SMAX = ABS(SX(1)) DO 30 I = 2,N XMAG = ABS(SX(I)) IF(XMAG.LE.SMAX) GO TO 30 ISAMAX = I SMAX = XMAG 30 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) SBLA 10 C C RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY. C REAL SX(1),SY(1) SDOT = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1)5,20,60 5 CONTINUE C C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SDOT = SDOT + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = N - (N/5)*5 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SDOT = SDOT + SX(I)*SY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) + * SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4) 50 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX SDOT = SDOT + SX(I)*SY(I) 70 CONTINUE RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) SBLA 10 C C OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY. C REAL SX(1),SY(1),SA IF(N.LE.0.OR.SA.EQ.0.E0) RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 20 M = N - (N/4)*4 IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I + 1) = SY(I + 1) + SA*SX(I + 1) SY(I + 2) = SY(I + 2) + SA*SX(I + 2) SY(I + 3) = SY(I + 3) + SA*SX(I + 3) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX SY(I) = SA*SX(I) + SY(I) 70 CONTINUE RETURN END SUBROUTINE SROTM (N,SX,INCX,SY,INCY,SPARAM) SBLA 10 C C APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX C C (SX(1) SX(N)) C ( ... ) C (SY(1) SY(N)) C C WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS.. C C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0 C C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0) C H=( ) ( ) ( ) ( ) C (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0). C DIMENSION SX(1),SY(1),SPARAM(5) DATA ZERO,TWO /0.E0,2.E0/ C SFLAG=SPARAM(1) IF(N .LE. 0 .OR.(SFLAG+TWO.EQ.ZERO)) GO TO 140 IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70 C NSTEPS=N*INCX IF(SFLAG) 50,10,30 10 CONTINUE SH12=SPARAM(4) SH21=SPARAM(3) DO 20 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W+Z*SH12 SY(I)=W*SH21+Z 20 CONTINUE GO TO 140 30 CONTINUE SH11=SPARAM(2) SH22=SPARAM(5) DO 40 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W*SH11+Z SY(I)=-W+SH22*Z 40 CONTINUE GO TO 140 50 CONTINUE SH11=SPARAM(2) SH12=SPARAM(4) SH21=SPARAM(3) SH22=SPARAM(5) DO 60 I=1,NSTEPS,INCX W=SX(I) Z=SY(I) SX(I)=W*SH11+Z*SH12 SY(I)=W*SH21+Z*SH22 60 CONTINUE GO TO 140 70 CONTINUE KX=1 KY=1 IF(INCX .LT. 0) KX=1+(1-N)*INCX IF(INCY .LT. 0) KY=1+(1-N)*INCY C IF(SFLAG)120,80,100 80 CONTINUE SH12=SPARAM(4) SH21=SPARAM(3) DO 90 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W+Z*SH12 SY(KY)=W*SH21+Z KX=KX+INCX KY=KY+INCY 90 CONTINUE GO TO 140 100 CONTINUE SH11=SPARAM(2) SH22=SPARAM(5) DO 110 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W*SH11+Z SY(KY)=-W+SH22*Z KX=KX+INCX KY=KY+INCY 110 CONTINUE GO TO 140 120 CONTINUE SH11=SPARAM(2) SH12=SPARAM(4) SH21=SPARAM(3) SH22=SPARAM(5) DO 130 I=1,N W=SX(KX) Z=SY(KY) SX(KX)=W*SH11+Z*SH12 SY(KY)=W*SH21+Z*SH22 KX=KX+INCX KY=KY+INCY 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE DROTMG (DD1,DD2,DX1,DY1,DPARAM) DBLA 10 C***BEGIN PROLOGUE DROTMG C***REVISION 3/1/80 C***CATEGORY NO. C***KEYWORD(S) C***AUTHOR--DATE WRITTEN C***PURPOSE C CONSTRUCT D.P. MODIFIED GIVENS TRANSFORMATION C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C C --INPUT-- C DD1 DOUBLE PRECISION SCALAR C DD2 DOUBLE PRECISION SCALAR C DX1 DOUBLE PRECISION SCALAR C DX2 DOUBLE PRECISION SCALAR C DPARAM D.P. 5-VECTOR. DPARAM(1)=DFLAG DEFINED BELOW. C ELEMENTS 2-5 DEFINE THE TRANSFORMATION MATRIX H. C C --OUTPUT-- C DD1 CHANGED TO REPRESENT THE EFFECT OF THE TRANSFORMATION C DD2 CHANGED TO REFLECT THE TRANSFORMATION C DX1 CHANGED TO REFLECT THE TRANSFORMATION C DX2 UNCHANGED C C CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS C THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)* C DY2)**T. C WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS.. C C DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0 C C (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0) C H=( ) ( ) ( ) ( ) C (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0). C LOCATIONS 2-5 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22 C RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE C VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.) C C***REFERENCE(S) C LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED NONE C***CARD COUNT IS 198 WITH 69 COMMENTS C***END PROLOGUE C DOUBLE PRECISION GAM,ONE,RGAMSQ,DD2,DH11,DH21,DPARAM,DP2, 1 DQ2,DU,DY1,ZERO,GAMSQ,DD1,DFLAG,DH12,DH22,DP1,DQ1, 2 DTEMP,DX1,TWO DIMENSION DPARAM(5) DATA ZERO,ONE,TWO /0.D0,1.D0,2.D0/ DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/ C***FIRST EXECUTABLE STATEMENT IF(.NOT. DD1 .LT. ZERO) GO TO 10 C GO ZERO-H-D-AND-DX1.. GO TO 60 10 CONTINUE C CASE-DD1-NONNEGATIVE DP2=DD2*DY1 IF(.NOT. DP2 .EQ. ZERO) GO TO 20 DFLAG=-TWO GO TO 260 C REGULAR-CASE.. 20 CONTINUE DP1=DD1*DX1 DQ2=DP2*DY1 DQ1=DP1*DX1 C IF(.NOT. DABS(DQ1) .GT. DABS(DQ2)) GO TO 40 DH21=-DY1/DX1 DH12=DP2/DP1 C DU=ONE-DH12*DH21 C IF(.NOT. DU .LE. ZERO) GO TO 30 C GO ZERO-H-D-AND-DX1.. GO TO 60 30 CONTINUE DFLAG=ZERO DD1=DD1/DU DD2=DD2/DU DX1=DX1*DU C GO SCALE-CHECK.. GO TO 100 40 CONTINUE IF(.NOT. DQ2 .LT. ZERO) GO TO 50 C GO ZERO-H-D-AND-DX1.. GO TO 60 50 CONTINUE DFLAG=ONE DH11=DP1/DP2 DH22=DX1/DY1 DU=ONE+DH11*DH22 DTEMP=DD2/DU DD2=DD1/DU DD1=DTEMP DX1=DY1*DU C GO SCALE-CHECK GO TO 100 C PROCEDURE..ZERO-H-D-AND-DX1.. 60 CONTINUE DFLAG=-ONE DH11=ZERO DH12=ZERO DH21=ZERO DH22=ZERO C DD1=ZERO DD2=ZERO DX1=ZERO C RETURN.. GO TO 220 C PROCEDURE..FIX-H.. 70 CONTINUE IF(.NOT. DFLAG .GE. ZERO) GO TO 90 C IF(.NOT. DFLAG .EQ. ZERO) GO TO 80 DH11=ONE DH22=ONE DFLAG=-ONE GO TO 90 80 CONTINUE DH21=-ONE DH12=ONE DFLAG=-ONE 90 CONTINUE GO TO IGO,(120,150,180,210) C PROCEDURE..SCALE-CHECK 100 CONTINUE 110 CONTINUE IF(.NOT. DD1 .LE. RGAMSQ) GO TO 130 IF(DD1 .EQ. ZERO) GO TO 160 ASSIGN 120 TO IGO C FIX-H.. GO TO 70 120 CONTINUE DD1=DD1*GAM**2 DX1=DX1/GAM DH11=DH11/GAM DH12=DH12/GAM GO TO 110 130 CONTINUE 140 CONTINUE IF(.NOT. DD1 .GE. GAMSQ) GO TO 160 ASSIGN 150 TO IGO C FIX-H.. GO TO 70 150 CONTINUE DD1=DD1/GAM**2 DX1=DX1*GAM DH11=DH11*GAM DH12=DH12*GAM GO TO 140 160 CONTINUE 170 CONTINUE IF(.NOT. DABS(DD2) .LE. RGAMSQ) GO TO 190 IF(DD2 .EQ. ZERO) GO TO 220 ASSIGN 180 TO IGO C FIX-H.. GO TO 70 180 CONTINUE DD2=DD2*GAM**2 DH21=DH21/GAM DH22=DH22/GAM GO TO 170 190 CONTINUE 200 CONTINUE IF(.NOT. DABS(DD2) .GE. GAMSQ) GO TO 220 ASSIGN 210 TO IGO C FIX-H.. GO TO 70 210 CONTINUE DD2=DD2/GAM**2 DH21=DH21*GAM DH22=DH22*GAM GO TO 200 220 CONTINUE IF(DFLAG)250,230,240 230 CONTINUE DPARAM(3)=DH21 DPARAM(4)=DH12 GO TO 260 240 CONTINUE DPARAM(2)=DH11 DPARAM(5)=DH22 GO TO 260 250 CONTINUE DPARAM(2)=DH11 DPARAM(3)=DH21 DPARAM(4)=DH12 DPARAM(5)=DH22 260 CONTINUE DPARAM(1)=DFLAG RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) DBLA 10 C***BEGIN PROLOGUE DCOPY C***REVISION 3/1/80 C***CATEGORY NO. C***KEYWORD(S) C***AUTHOR--DATE WRITTEN 10/79,LAWSON C. (JPL),HANSON R. (SLA), C KINCAID D. (U TEXAS), KROGH F. (JPL) C***PURPOSE C D.P. VECTOR COPY Y = X C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C C --OUTPUT-- C DY COPY OF VECTOR DX (UNCHANGED IF N.LE.0) C C COPY DOUBLE PRECISION DX TO DOUBLE PRECISION DY. C FOR I = 0 TO N-1, COPY DX(LX+I*INCX) TO DY(LY+I*INCY), C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS C DEFINED IN A SIMILAR WAY USING INCY. C C***REFERENCE(S) C LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED NONE C***CARD COUNT IS 88 WITH 49 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1),DY(1) C***FIRST EXECUTABLE STATEMENT IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX DY(I) = DX(I) 70 CONTINUE RETURN END SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) DBLA 10 C***BEGIN PROLOGUE DSWAP C***REVISION 3/1/80 C***CATEGORY NO. C***KEYWORD(S) C***AUTHOR--DATE WRITTEN 10/79,LAWSON C. (JPL),HANSON R. (SLA), C KINCAID D. (U TEXAS), KROGH F. (JPL) C***PURPOSE C INTERCHANGE D.P. VECTORS C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C DY DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCY STORAGE SPACING BETWEEN ELEMENTS OF DY C C --OUTPUT-- C DX INPUT VECTOR DY (UNCHANGED IF N.LE.0) C DY INPUT VECTOR DX (UNCHANGED IF N.LE.0) C C INTERCHANGE DOUBLE PRECISION DX AND DOUBLE PRECISION DY. C FOR I = 0 TO N-1, INTERCHANGE DX(LX+I*INCX) AND DY(LY+I*INCY), C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS C DEFINED IN A SIMILAR WAY USING INCY. C C***REFERENCE(S) C LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED NONE C***CARD COUNT IS 97 WITH 50 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1),DY(1),DTEMP1,DTEMP2,DTEMP3 C***FIRST EXECUTABLE STATEMENT IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP1 = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP1 IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3. C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP1 = DX(I) DTEMP2 = DX(I+1) DTEMP3 = DX(I+2) DX(I) = DY(I) DX(I+1) = DY(I+1) DX(I+2) = DY(I+2) DY(I) = DTEMP1 DY(I+1) = DTEMP2 DY(I+2) = DTEMP3 50 CONTINUE RETURN 60 CONTINUE C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C NS = N*INCX DO 70 I=1,NS,INCX DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) DBLA 10 C***BEGIN PROLOGUE DNRM2 C***REVISION 3/1/80 C***CATEGORY NO. C***KEYWORD(S) C***AUTHOR--DATE WRITTEN 10/79,LAWSON C. (JPL),HANSON R. (SLA), C KINCAID D. (U TEXAS), KROGH F. (JPL) C***PURPOSE C EUCLIDEAN LENGTH (L2 NORM) OF D.P. VECTOR C***DESCRIPTION C B L A S SUBPROGRAM C DESCRIPTION OF PARAMETERS C C --INPUT-- C N NUMBER OF ELEMENTS IN INPUT VECTOR(S) C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C C --OUTPUT-- C DNRM2 DOUBLE PRECISION RESULT (ZERO IF N.LE.0) C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS