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 THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(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 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 151 WITH 105 COMMENTS C***END PROLOGUE INTEGER NEXT DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C***FIRST EXECUTABLE STATEMENT IF(N .GT. 0) GO TO 10 DNRM2 = 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( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(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 / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(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( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(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(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C TERM. OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) DBLA 10 C***BEGIN PROLOGUE DASUM 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 SUM OF MAGNITUDES OF D.P. VECTOR COMPONENTS 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 DASUM DOUBLE PRECISION RESULT (ZERO IF N.LE.0) C C RETURNS SUM OF MAGNITUDES OF DOUBLE PRECISION DX. C DASUM = SUM FROM 0 TO N-1 OF DABS(DX(1+I*INCX)) 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 65 WITH 42 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1) C***FIRST EXECUTABLE STATEMENT DASUM = 0.D0 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 DASUM = DASUM + DABS(DX(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 = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DASUM = DASUM + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 DASUM = DASUM + DABS(DX(I)) + DABS(DX(I+1)) + DABS(DX(I+2)) $ + DABS(DX(I+3)) + DABS(DX(I+4)) + DABS(DX(I+5)) 50 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) DBLA 10 C***BEGIN PROLOGUE DSCAL 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 SCALE X = A*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 DA DOUBLE PRECISION SCALE FACTOR C DX DOUBLE PRECISION VECTOR WITH N ELEMENTS C INCX STORAGE SPACING BETWEEN ELEMENTS OF DX C C --OUTPUT-- C DX DOUBLE PRECISION RESULT (UNCHANGED IF N.LE.0) C C REPLACE DOUBLE PRECISION DX BY DOUBLE PRECISION DA*DX. C FOR I = 0 TO N-1, REPLACE DX(1+I*INCX) WITH DA * DX(1+I*INCX) 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 68 WITH 43 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DA,DX(1) C***FIRST EXECUTABLE STATEMENT 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 DX(I) = DA*DX(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 = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END INTEGER FUNCTION IDAMAX(N,DX,INCX) DBLA 10 C***BEGIN PROLOGUE IDAMAX 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 FIND LARGEST COMPONENT 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 IDAMAX SMALLEST INDEX (ZERO IF N.LE.0) C C FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF DOUBLE PRECISION DX. C IDAMAX = FIRST I, I = 1 TO N, TO MINIMIZE ABS(DX(1-INCX+I*INCX) 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 66 WITH 39 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1),DMAX,XMAG C***FIRST EXECUTABLE STATEMENT IDAMAX = 0 IF(N.LE.0) RETURN IDAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C DMAX = DABS(DX(1)) NS = N*INCX II = 1 DO 10 I = 1,NS,INCX XMAG = DABS(DX(I)) IF(XMAG.LE.DMAX) GO TO 5 IDAMAX = II DMAX = XMAG 5 II = II + 1 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C 20 DMAX = DABS(DX(1)) DO 30 I = 2,N XMAG = DABS(DX(I)) IF(XMAG.LE.DMAX) GO TO 30 IDAMAX = I DMAX = XMAG 30 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) DBLA 10 C***BEGIN PROLOGUE DDOT 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. INNER PRODUCT OF 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 DDOT DOUBLE PRECISION DOT PRODUCT (ZERO IF N.LE.0) C C RETURNS THE DOT PRODUCT OF DOUBLE PRECISION DX AND DY. C DDOT = SUM FOR I = 0 TO N-1 OF DX(LX+I*INCX) * 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 84 WITH 49 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1),DY(1) C***FIRST EXECUTABLE STATEMENT DDOT = 0.D0 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 DDOT = DDOT + DX(IX)*DY(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 = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DDOT = DDOT + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + $ DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(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 DDOT = DDOT + DX(I)*DY(I) 70 CONTINUE RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) DBLA 10 C***BEGIN PROLOGUE DAXPY 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 COMPUTATION Y = A*X + Y 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 DA DOUBLE PRECISION SCALAR MULTIPLIER 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 DOUBLE PRECISION RESULT (UNCHANGED IF N.LE.0) C C OVERWRITE DOUBLE PRECISION DY WITH DOUBLE PRECISION DA*DX + DY. C FOR I = 0 TO N-1, REPLACE DY(LY+I*INCY) WITH DA*DX(LX+I*INCX) + C DY(LY+I*INCY), WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N C AND LY IS 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 86 WITH 50 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DX(1),DY(1),DA C***FIRST EXECUTABLE STATEMENT IF(N.LE.0.OR.DA.EQ.0.D0) 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 DY(IY) = DY(IY) + DA*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 4. C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(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 DY(I) = DA*DX(I) + DY(I) 70 CONTINUE RETURN END SUBROUTINE DROTM (N,DX,INCX,DY,INCY,DPARAM) DBLA 10 C***BEGIN PROLOGUE DROTM 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 APPLY D.P. MODIFIED GIVENS TRANSFORMATION 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 DPARAM 5-ELEMENT D.P. VECTOR. DPARAM(1) IS DFLAG DESCRIBED BELOW C ELEMENTS 2-5 FORM THE TRANSFORMATION MATRIX H. C C --OUTPUT-- C DX ROTATED VECTOR (UNCHANGED IF N.LE.0) C DY ROTATED VECTOR (UNCHANGED IF N.LE.0) C C APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX C C (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN C (DY**T) C C DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE C LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY. 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 SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM. C 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 142 WITH 54 COMMENTS C***END PROLOGUE C DOUBLE PRECISION DFLAG,DH12,DH22,DX,TWO,Z,DH11,DH21, 1 DPARAM,DY,W,ZERO DIMENSION DX(1),DY(1),DPARAM(5) DATA ZERO,TWO/0.D0,2.D0/ C***FIRST EXECUTABLE STATEMENT DFLAG=DPARAM(1) IF(N .LE. 0 .OR.(DFLAG+TWO.EQ.ZERO)) GO TO 140 IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70 C NSTEPS=N*INCX IF(DFLAG) 50,10,30 10 CONTINUE DH12=DPARAM(4) DH21=DPARAM(3) DO 20 I=1,NSTEPS,INCX W=DX(I) Z=DY(I) DX(I)=W+Z*DH12 DY(I)=W*DH21+Z 20 CONTINUE GO TO 140 30 CONTINUE DH11=DPARAM(2) DH22=DPARAM(5) DO 40 I=1,NSTEPS,INCX W=DX(I) Z=DY(I) DX(I)=W*DH11+Z DY(I)=-W+DH22*Z 40 CONTINUE GO TO 140 50 CONTINUE DH11=DPARAM(2) DH12=DPARAM(4) DH21=DPARAM(3) DH22=DPARAM(5) DO 60 I=1,NSTEPS,INCX W=DX(I) Z=DY(I) DX(I)=W*DH11+Z*DH12 DY(I)=W*DH21+Z*DH22 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(DFLAG)120,80,100 80 CONTINUE DH12=DPARAM(4) DH21=DPARAM(3) DO 90 I=1,N W=DX(KX) Z=DY(KY) DX(KX)=W+Z*DH12 DY(KY)=W*DH21+Z KX=KX+INCX KY=KY+INCY 90 CONTINUE GO TO 140 100 CONTINUE DH11=DPARAM(2) DH22=DPARAM(5) DO 110 I=1,N W=DX(KX) Z=DY(KY) DX(KX)=W*DH11+Z DY(KY)=-W+DH22*Z KX=KX+INCX KY=KY+INCY 110 CONTINUE GO TO 140 120 CONTINUE DH11=DPARAM(2) DH12=DPARAM(4) DH21=DPARAM(3) DH22=DPARAM(5) DO 130 I=1,N W=DX(KX) Z=DY(KY) DX(KX)=W*DH11+Z*DH12 DY(KY)=W*DH21+Z*DH22 KX=KX+INCX KY=KY+INCY 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE FDUMP XER 10 C ABSTRACT C ***NOTE*** MACHINE DEPENDENT ROUTINE C FDUMP IS INTENDED TO BE REPLACED BY A LOCALLY WRITTEN C VERSION WHICH PRODUCES A SYMBOLIC DUMP. FAILING THIS, C IT SHOULD BE REPLACED BY A VERSION WHICH PRINTS THE C SUBPROGRAM NESTING LIST. NOTE THAT THIS DUMP MUST BE C PRINTED ON EACH OF UP TO FIVE FILES, AS INDICATED BY THE C XGETUA ROUTINE. SEE XSETUA AND XGETUA FOR DETAILS. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 23 MAY 1979 C RETURN END SUBROUTINE XERABT(MESSG,NMESSG) XER 10 C C ABSTRACT C ***NOTE*** MACHINE DEPENDENT ROUTINE C XERABT ABORTS THE EXECUTION OF THE PROGRAM. C THE ERROR MESSAGE CAUSING THE ABORT IS GIVEN IN THE CALLING C SEQUENCE IN CASE ONE NEEDS IT FOR PRINTING ON A DAYFILE, C FOR EXAMPLE. C C DESCRIPTION OF PARAMETERS C MESSG AND NMESSG ARE AS IN XERROR, EXCEPT THAT NMESSG MAY C BE ZERO, IN WHICH CASE NO MESSAGE IS BEING SUPPLIED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C DIMENSION MESSG(NMESSG) IF (.TRUE.) STOP RETURN END FUNCTION J4SAVE(IWHICH,IVALUE,ISET) XER 10 C C ABSTRACT C J4SAVE SAVES AND RECALLS SEVERAL GLOBAL VARIABLES NEEDED C BY THE LIBRARY ERROR HANDLING ROUTINES. C C DESCRIPTION OF PARAMETERS C --INPUT-- C IWHICH - INDEX OF ITEM DESIRED. C = 1 REFERS TO CURRENT ERROR NUMBER. C = 2 REFERS TO CURRENT ERROR CONTROL FLAG. C = 3 REFERS TO CURRENT UNIT NUMBER TO WHICH ERROR C MESSAGES ARE TO BE SENT. (0 MEANS USE STANDARD.) C = 4 REFERS TO THE MAXIMUM NUMBER OF TIMES ANY C MESSAGE IS TO BE PRINTED (AS SET BY XERMAX). C = 5 REFERS TO THE TOTAL NUMBER OF UNITS TO WHICH C EACH ERROR MESSAGE IS TO BE WRITTEN. C = 6 REFERS TO THE 2ND UNIT FOR ERROR MESSAGES C = 7 REFERS TO THE 3RD UNIT FOR ERROR MESSAGES C = 8 REFERS TO THE 4TH UNIT FOR ERROR MESSAGES C = 9 REFERS TO THE 5TH UNIT FOR ERROR MESSAGES C IVALUE - THE VALUE TO BE SET FOR THE IWHICH-TH PARAMETER, C IF ISET IS .TRUE. . C ISET - IF ISET=.TRUE., THE IWHICH-TH PARAMETER WILL BE C GIVEN THE VALUE, IVALUE. IF ISET=.FALSE., THE C IWHICH-TH PARAMETER WILL BE UNCHANGED, AND IVALUE C IS A DUMMY PARAMETER. C --OUTPUT-- C THE (OLD) VALUE OF THE IWHICH-TH PARAMETER WILL BE RETURNED C IN THE FUNCTION VALUE, J4SAVE. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C ADAPTED FROM BELL LABORATORIES PORT LIBRARY ERROR HANDLER C END OF ABSTRACT C LATEST REVISION --- 23 MAY 1979 C LOGICAL ISET INTEGER IPARAM(9) DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/ DATA IPARAM(5)/1/ DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ J4SAVE = IPARAM(IWHICH) IF (ISET) IPARAM(IWHICH) = IVALUE RETURN END FUNCTION NUMXER(NERR) XER 10 C C ABSTRACT C NUMXER RETURNS THE MOST RECENT ERROR NUMBER, C IN BOTH NUMXER AND THE PARAMETER NERR. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C NERR = J4SAVE(1,0,.FALSE.) NUMXER = NERR RETURN END SUBROUTINE S88FMT(N,IVALUE,IFMT) XER 10 C C ABSTRACT C S88FMT REPLACES IFMT(1), ... ,IFMT(N) WITH THE C CHARACTERS CORRESPONDING TO THE N LEAST SIGNIFICANT C DIGITS OF IVALUE. C C TAKEN FROM THE BELL LABORATORIES PORT LIBRARY ERROR HANDLER C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C DIMENSION IFMT(N),IDIGIT(10) DATA IDIGIT(1),IDIGIT(2),IDIGIT(3),IDIGIT(4),IDIGIT(5), 1 IDIGIT(6),IDIGIT(7),IDIGIT(8),IDIGIT(9),IDIGIT(10) 2 /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ NT = N IT = IVALUE 10 IF (NT .EQ. 0) RETURN INDEX = MOD(IT,10) IFMT(NT) = IDIGIT(INDEX+1) IT = IT/10 NT = NT - 1 GO TO 10 END SUBROUTINE XERCLR XER 10 C C ABSTRACT C THIS ROUTINE SIMPLY RESETS THE CURRENT ERROR NUMBER TO ZERO. C THIS MAY BE NECESSARY TO DO IN ORDER TO DETERMINE THAT C A CERTAIN ERROR HAS OCCURRED AGAIN SINCE THE LAST TIME C NUMXER WAS REFERENCED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C JUNK = J4SAVE(1,0,.TRUE.) RETURN END SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL) XER 10 C C ABSTRACT C ALLOWS USER CONTROL OVER HANDLING OF INDIVIDUAL ERRORS. C JUST AFTER EACH MESSAGE IS RECORDED, BUT BEFORE IT IS C PROCESSED ANY FURTHER (I.E., BEFORE IT IS PRINTED OR C A DECISION TO ABORT IS MADE) A CALL IS MADE TO XERCTL. C IF THE USER HAS PROVIDED HIS OWN VERSION OF XERCTL, HE C CAN THEN OVERRIDE THE VALUE OF KONTROL USED IN PROCESSING C THIS MESSAGE BY REDEFINING ITS VALUE. C KONTRL MAY BE SET TO ANY VALUE FROM -2 TO 2. C THE MEANINGS FOR KONTRL ARE THE SAME AS IN XSETF, EXCEPT C THAT THE VALUE OF KONTRL CHANGES ONLY FOR THIS MESSAGE. C IF KONTRL IS SET TO A VALUE OUTSIDE THE RANGE FROM -2 TO 2, C IT WILL BE MOVED BACK INTO THAT RANGE. C C DESCRIPTION OF PARAMETERS C C --INPUT-- C MESSG1 - THE FIRST WORD (ONLY) OF THE ERROR MESSAGE. C NMESSG - SAME AS IN THE CALL TO XERROR OR XERRWV. C NERR - SAME AS IN THE CALL TO XERROR OR XERRWV. C LEVEL - SAME AS IN THE CALL TO XERROR OR XERRWV. C KONTRL - THE CURRENT VALUE OF THE CONTROL FLAG AS SET C BY A CALL TO XSETF. C C --OUTPUT-- C KONTRL - THE NEW VALUE OF KONTRL. IF KONTRL IS NOT C DEFINED, IT WILL REMAIN AT ITS ORIGINAL VALUE. C THIS CHANGED VALUE OF CONTROL AFFECTS ONLY C THE CURRENT OCCURRENCE OF THE CURRENT MESSAGE. C END OF ABSTRACT C RETURN END SUBROUTINE XERDMP XER 10 C C ABSTRACT C XERDMP PRINTS AN ERROR TABLE SHOWING ALL ERRORS WHICH C HAVE OCCURRED DURING THE CURRENT EXECUTION, OR SINCE XERDMP C WAS LAST CALLED. AFTER PRINTING, THE ERROR TABLE IS CLEARED, C AND IF PROGRAM EXECUTION IS CONTINUED ACCUMULATION OF THE C ERROR TABLE BEGINS AT ZERO. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C CALL XERSAV(1H ,0,0,0,KOUNT) RETURN END SUBROUTINE XERMAX(MAX) XER 10 C C ABSTRACT C XERMAX SETS THE MAXIMUM NUMBER OF TIMES ANY MESSAGE C IS TO BE PRINTED. THAT IS, NON-FATAL MESSAGES ARE C NOT TO BE PRINTED AFTER THEY HAVE OCCURED MAX TIMES. C SUCH NON-FATAL MESSAGES MAY BE PRINTED LESS THAN C MAX TIMES EVEN IF THEY OCCUR MAX TIMES, IF ERROR C SUPPRESSION MODE (KONTRL=0) IS EVER IN EFFECT. C C THE DEFAULT VALUE FOR MAX IS 10. C C DESCRIPTION OF PARAMETER C --INPUT-- C MAX - THE MAXIMUM NUMBER OF TIMES ANY ONE MESSAGE C IS TO BE PRINTED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C JUNK = J4SAVE(4,MAX,.TRUE.) RETURN END SUBROUTINE XERPRT(MESSG,NMESSG) XER 10 C C ABSTRACT C PRINT THE HOLLERITH MESSAGE IN MESSG, OF LENGTH MESSG, C ON EACH FILE INDICATED BY XGETUA. C THIS VERSION PRINTS EXACTLY THE RIGHT NUMBER OF CHARACTERS, C NOT A NUMBER OF WORDS, AND THUS SHOULD WORK ON MACHINES C WHICH DO NOT BLANK FILL THE LAST WORD OF THE HOLLERITH. C C RON JONES, JUNE 1980 C END OF ABSTRACT C INTEGER F(10),G(14),LUN(5) DIMENSION MESSG(NMESSG) DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10) 1 / 1H( ,1H1 ,1HX ,1H, ,1H ,1H ,1HA ,1H ,1H ,1H) / DATA G(1),G(2),G(3),G(4),G(5),G(6),G(7),G(8),G(9),G(10) 1 / 1H( ,1H1 ,1HX ,1H ,1H ,1H ,1H ,1H ,1H ,1H / DATA G(11),G(12),G(13),G(14) 1 / 1H ,1H ,1H ,1H) / DATA LA/1HA/,LCOM/1H,/,LBLANK/1H / C PREPARE FORMAT FOR WHOLE LINES NCHAR = I1MACH(6) NFIELD = 72/NCHAR CALL S88FMT(2,NFIELD,F(5)) CALL S88FMT(2,NCHAR,F(8)) C PREPARE FORMAT FOR LAST, PARTIAL LINE, IF NEEDED NCHARL = NFIELD*NCHAR NLINES = NMESSG/NCHARL NWORD = NLINES*NFIELD NCHREM = NMESSG - NLINES*NCHARL IF (NCHREM.LE.0) GO TO 40 DO 10 I=4,13 10 G(I) = LBLANK NFIELD = NCHREM/NCHAR IF (NFIELD.LE.0) GO TO 20 C PREPARE WHOLE WORD FIELDS G(4) = LCOM CALL S88FMT(2,NFIELD,G(5)) G(7) = LA CALL S88FMT(2,NCHAR,G(8)) 20 CONTINUE NCHLST = MOD(NCHREM,NCHAR) IF (NCHLST.LE.0) GO TO 30 C PREPARE PARTIAL WORD FIELD G(10) = LCOM G(11) = LA CALL S88FMT(2,NCHLST,G(12)) 30 CONTINUE 40 CONTINUE C PRINT THE MESSAGE NWORD1 = NWORD+1 NWORD2 = (NMESSG+NCHAR-1)/NCHAR CALL XGETUA(LUN,NUNIT) DO 50 KUNIT = 1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NWORD.GT.0) WRITE (IUNIT,F) (MESSG(I),I=1,NWORD) IF (NCHREM.GT.0) WRITE (IUNIT,G) (MESSG(I),I=NWORD1,NWORD2) 50 CONTINUE RETURN END SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) XER 10 C C ABSTRACT C XERROR PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER C DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE C OF THE LIBRARY ERROR CONTROL FLAG, KONTRL. C (SEE SUBROUTINE XSETF FOR DETAILS.) C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED, CONTAINING C NO MORE THAN 72 CHARACTERS. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C NERR MUST NOT BE ZERO. C LEVEL - ERROR CATEGORY. C =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR. C =1 MEANS THIS IS A RECOVERABLE ERROR. (I.E., IT IS C NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.) C =0 MEANS THIS IS A WARNING MESSAGE ONLY. C =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE C PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY C TIMES THIS CALL IS EXECUTED. C C EXAMPLES C CALL XERROR(23HSMOOTH -- NUM WAS ZERO.,23,1,2) C CALL XERROR(43HINTEG -- LESS THAN FULL ACCURACY ACHIEVED., C 43,2,1) C CALL XERROR(65HROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL C 1 FULLY COLLAPSED.,65,3,0) C CALL XERROR(39HEXP -- UNDERFLOWS BEING SET TO ZERO.,39,1,-1) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C REVISED BY K HASKELL TO CHECK INPUT ARGS, 2/18/80 C DIMENSION MESSG(NMESSG) C CHECK FOR VALID INPUT LKNTRL = J4SAVE (2,0,.FALSE.) IF (NMESSG.GT.0) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (33HXERROR -- NMESSG MUST BE POSITIVE,33) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV (1H ,0,0,0,KDUMMY) CALL XERABT (23HXERROR -- INVALID INPUT,23) RETURN 10 CONTINUE IF (NERR.NE.0) GO TO 15 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (28HXERROR -- NERR=0 IS AN ERROR,28) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV (1H ,0,0,0,KDUMMY) CALL XERABT (23HXERROR -- INVALID INPUT,23) RETURN 15 CONTINUE IF ((LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 20 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (32HXERROR -- INVALID VALUE OF LEVEL,32) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV (1H ,0,0,0,KDUMMY) CALL XERABT (23HXERROR -- INVALID INPUT,23) RETURN 20 CONTINUE CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) RETURN END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) XER 10 C C ABSTRACT C XERRWV PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER C DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE C OF THE LIBRARY ERROR CONTROL FLAG, KONTRL. C (SEE SUBROUTINE XSETF FOR DETAILS.) C IN ADDITION, UP TO TWO INTEGER VALUES AND TWO REAL C VALUES MAY BE PRINTED ALONG WITH THE MESSAGE. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED. C NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG. C NERR - THE ERROR NUMBER ASSOCIATED WITH THIS MESSAGE. C NERR MUST NOT BE ZERO. C LEVEL - ERROR CATEGORY. C =2 MEANS THIS IS AN UNCONDITIONALLY FATAL ERROR. C =1 MEANS THIS IS A RECOVERABLE ERROR. (I.E., IT IS C NON-FATAL IF XSETF HAS BEEN APPROPRIATELY CALLED.) C =0 MEANS THIS IS A WARNING MESSAGE ONLY. C =-1 MEANS THIS IS A WARNING MESSAGE WHICH IS TO BE C PRINTED AT MOST ONCE, REGARDLESS OF HOW MANY C TIMES THIS CALL IS EXECUTED. C NI - NUMBER OF INTEGER VALUES TO BE PRINTED. (O TO 2) C I1 - FIRST INTEGER VALUE. C I2 - SECOND INTEGER VALUE. C NR - NUMBER OF REAL VALUES TO BE PRINTED. (0 TO 2) C R1 - FIRST REAL VALUE. C R2 - SECOND REAL VALUE. C C EXAMPLES C CALL XERRWV(29HSMOOTH -- NUM (=I1) WAS ZERO.,29,1,2, C 1 1,NUM,0,0,0.,0.) C CALL XERRWV(54HQUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM C 1 (R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 19 MAR 1980 C REVISED BY K HASKELL TO CHECK INPUT ARGS, 2/18/80 C DIMENSION MESSG(NMESSG),LUN(5) C GET FLAGS LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) C CHECK FOR VALID INPUT IF (NMESSG.GT.0) GO TO 2 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (33HXERRWV -- NMESSG MUST BE POSITIVE,33) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY) CALL XERABT (23HXERRWV -- INVALID INPUT,23) RETURN 2 CONTINUE IF (NERR.NE.0) GO TO 4 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (28HXERRWV -- NERR=0 IS AN ERROR,28) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY) CALL XERABT (23HXERRWV -- INVALID INPUT,23) RETURN 4 CONTINUE IF ((LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT(17HFATAL ERROR IN...,17) CALL XERPRT (32HXERRWV -- INVALID VALUE OF LEVEL,32) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT(29HJOB ABORT DUE TO FATAL ERROR., 1 29) IF (LKNTRL.GT.0) CALL XERSAV(1H ,0,0,0,KDUMMY) CALL XERABT(23HXERROR -- INVALID INPUT,23) RETURN 10 CONTINUE C RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) C LET USER OVERRIDE LFIRST = MESSG(1) LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) C RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) C DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(1H ,1) C INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT 1(57HWARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.,57) IF (LLEVEL.EQ.0) CALL XERPRT(13HWARNING IN...,13) IF (LLEVEL.EQ.1) CALL XERPRT 1 (23HRECOVERABLE ERROR IN...,23) IF (LLEVEL.EQ.2) CALL XERPRT(17HFATAL ERROR IN...,17) 20 CONTINUE C MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) IF (NI.GE.1) WRITE (IUNIT,22) I1 IF (NI.GE.2) WRITE (IUNIT,23) I2 IF (NR.GE.1) WRITE (IUNIT,24) R1 IF (NR.GE.2) WRITE (IUNIT,25) R2 22 FORMAT (11X,21HIN ABOVE MESSAGE, I1=,I10) 23 FORMAT (11X,21HIN ABOVE MESSAGE, I2=,I10) 24 FORMAT (11X,21HIN ABOVE MESSAGE, R1=,E20.10) 25 FORMAT (11X,21HIN ABOVE MESSAGE, R2=,E20.10) IF (LKNTRL.LE.0) GO TO 40 C ERROR NUMBER WRITE (IUNIT,30) LERR 30 FORMAT (15H ERROR NUMBER =,I10) 40 CONTINUE 50 CONTINUE C TRACE-BACK CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) 1IFATAL = 1 C QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF (LKNTRL.LE.0) GO TO 120 C PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT 1 (35HJOB ABORT DUE TO UNRECOVERED ERROR.,35) IF (LLEVEL.EQ.2) CALL XERPRT 1 (29HJOB ABORT DUE TO FATAL ERROR.,29) C PRINT ERROR SUMMARY CALL XERSAV(1H ,0,0,0,KDUMMY) 120 CONTINUE C ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) XER 10 C C ABSTRACT C RECORD THAT THIS ERROR OCCURRED. C C DESCRIPTION OF PARAMETERS C --INPUT-- C MESSG, NMESSG, NERR, LEVEL ARE AS IN XERROR, C EXCEPT THAT WHEN NMESSG=0 THE TABLES WILL BE C DUMPED AND CLEARED, AND WHEN NMESSG IS LESS THAN ZERO THE C TABLES WILL BE DUMPED AND NOT CLEARED. C --OUTPUT-- C ICOUNT WILL BE THE NUMBER OF TIMES THIS MESSAGE HAS C BEEN SEEN, OR ZERO IF THE TABLE HAS OVERFLOWED AND C DOES NOT CONTAIN THIS MESSAGE SPECIFICALLY. C WHEN NMESSG=0, ICOUNT WILL NOT BE ALTERED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C INTEGER F(17),LUN(5) DIMENSION MESSG(1) DIMENSION MESTAB(10),NERTAB(10),LEVTAB(10),KOUNT(10) C NEXT THREE DATA STATEMENTS ARE NEEDED MERELY TO SATISFY C CERTAIN CONVENTIONS FOR COMPILERS WHICH DYNAMICALLY C ALLOCATE STORAGE. DATA MESTAB(1),MESTAB(2),MESTAB(3),MESTAB(4),MESTAB(5), 1 MESTAB(6),MESTAB(7),MESTAB(8),MESTAB(9),MESTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA NERTAB(1),NERTAB(2),NERTAB(3),NERTAB(4),NERTAB(5), 1 NERTAB(6),NERTAB(7),NERTAB(8),NERTAB(9),NERTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA LEVTAB(1),LEVTAB(2),LEVTAB(3),LEVTAB(4),LEVTAB(5), 1 LEVTAB(6),LEVTAB(7),LEVTAB(8),LEVTAB(9),LEVTAB(10) 2 /0,0,0,0,0,0,0,0,0,0/ C NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK C ERROR TABLE INITIALLY DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNTX/0/ C NEXT DATA STATEMENT SETS UP OUTPUT FORMAT DATA F(1),F(2),F(3),F(4),F(5),F(6),F(7),F(8),F(9),F(10), 1 F(11),F(12),F(13),F(14),F(15),F(16),F(17) 2 /1H( ,1H1 ,1HX ,1H, ,1HA ,1H ,1H ,1H, ,1HI ,1H , 3 1H ,1H, ,1H2 ,1HI ,1H1 ,1H0 ,1H) / IF (NMESSG.GT.0) GO TO 80 C DUMP THE TABLE IF (KOUNT(1).EQ.0) RETURN C PREPARE FORMAT NCHAR = I1MACH(6) CALL S88FMT(2,NCHAR,F(6)) NCOL = 20 - NCHAR CALL S88FMT(2,NCOL,F(10)) C PRINT TO EACH UNIT CALL XGETUA(LUN,NUNIT) DO 60 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C PRINT TABLE HEADER WRITE (IUNIT,10) 10 FORMAT (32H0 ERROR MESSAGE SUMMARY/ 1 41H FIRST WORD NERR LEVEL COUNT) C PRINT BODY OF TABLE DO 20 I=1,10 IF (KOUNT(I).EQ.0) GO TO 30 WRITE (IUNIT,F) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) 20 CONTINUE 30 CONTINUE C PRINT NUMBER OF OTHER ERRORS IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX 40 FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10) WRITE (IUNIT,50) 50 FORMAT (1X) 60 CONTINUE IF (NMESSG.LT.0) RETURN C CLEAR THE ERROR TABLES DO 70 I=1,10 70 KOUNT(I) = 0 KOUNTX = 0 RETURN 80 CONTINUE C PROCESS A MESSAGE... C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. DO 90 I=1,10 II = I IF (KOUNT(I).EQ.0) GO TO 110 IF (MESSG(1).NE.MESTAB(I)) GO TO 90 IF (NERR.NE.NERTAB(I)) GO TO 90 IF (LEVEL.NE.LEVTAB(I)) GO TO 90 GO TO 100 90 CONTINUE C THREE POSSIBLE CASES... C TABLE IS FULL KOUNTX = KOUNTX+1 ICOUNT = 1 RETURN C MESSAGE FOUND IN TABLE 100 KOUNT(II) = KOUNT(II) + 1 ICOUNT = KOUNT(II) RETURN C EMPTY SLOT FOUND FOR NEW MESSAGE 110 MESTAB(II) = MESSG(1) NERTAB(II) = NERR LEVTAB(II) = LEVEL KOUNT(II) = 1 ICOUNT = 1 RETURN END SUBROUTINE XGETF(KONTRL) XER 10 C C ABSTRACT C XGETF RETURNS THE CURRENT VALUE OF THE ERROR CONTROL FLAG C IN KONTRL. SEE SUBROUTINE XSETF FOR FLAG VALUE MEANINGS. C (KONTRL IS AN OUTPUT PARAMETER ONLY.) C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C KONTRL = J4SAVE(2,0,.FALSE.) RETURN END SUBROUTINE XGETUA(IUNIT,N) XER 10 C C ABSTRACT C XGETUA MAY BE CALLED TO DETERMINE THE UNIT NUMBER OR NUMBERS C TO WHICH ERROR MESSAGES ARE BEING SENT. C THESE UNIT NUMBERS MAY HAVE BEEN SET BY A CALL TO XSETUN, C OR A CALL TO XSETUA, OR MAY BE A DEFAULT VALUE. C C DESCRIPTION OF PARAMETERS C --OUTPUT-- C IUNIT - AN ARRAY OF ONE TO FIVE UNIT NUMBERS, DEPENDING C ON THE VALUE OF N. A VALUE OF ZERO REFERS TO THE C DEFAULT UNIT, AS DEFINED BY THE I1MACH MACHINE C CONSTANT ROUTINE. ONLY IUNIT(1),...,IUNIT(N) ARE C DEFINED BY XGETUA. THE VALUES OF IUNIT(N+1),..., C IUNIT(5) ARE NOT DEFINED (FOR N.LT.5) OR ALTERED C IN ANY WAY BY XGETUA. C N - THE NUMBER OF UNITS TO WHICH COPIES OF THE C ERROR MESSAGES ARE BEING SENT. N WILL BE IN THE C RANGE FROM 1 TO 5. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C DIMENSION IUNIT(5) N = J4SAVE(5,0,.FALSE.) DO 30 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 IUNIT(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END SUBROUTINE XGETUN(IUNIT) XER 10 C C ABSTRACT C XGETUN GETS THE (FIRST) OUTPUT FILE TO WHICH ERROR MESSAGES C ARE BEING SENT. TO FIND OUT IF MORE THAN ONE FILE IS BEING C USED, ONE MUST USE THE XGETUA ROUTINE. C C DESCRIPTION OF PARAMETER C --OUTPUT-- C IUNIT - THE LOGICAL UNIT NUMBER OF THE (FIRST) UNIT TO C WHICH ERROR MESSAGES ARE BEING SENT. C A VALUE OF ZERO MEANS THAT THE DEFAULT FILE, AS C DEFINED BY THE I1MACH ROUTINE, IS BEING USED. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 23 MAY 1979 C IUNIT = J4SAVE(3,0,.FALSE.) RETURN END SUBROUTINE XSETF(KONTRL) XER 10 C C ABSTRACT C XSETF SETS THE ERROR CONTROL FLAG VALUE TO KONTRL. C (KONTRL IS AN INPUT PARAMETER ONLY.) C THE FOLLOWING TABLE SHOWS HOW EACH MESSAGE IS TREATED, C DEPENDING ON THE VALUES OF KONTRL AND LEVEL. (SEE XERROR C FOR DESCRIPTION OF LEVEL.) C C IF KONTRL IS ZERO OR NEGATIVE, NO INFORMATION OTHER THAN THE C MESSAGE ITSELF (INCLUDING NUMERIC VALUES, IF ANY) WILL BE C PRINTED. IF KONTRL IS POSITIVE, INTRODUCTORY MESSAGES, C TRACE-BACKS, ETC., WILL BE PRINTED IN ADDITION TO THE MESSAGE. C C IABS(KONTRL) C LEVEL 0 1 2 C VALUE C 2 FATAL FATAL FATAL C C 1 NOT PRINTED PRINTED FATAL C C 0 NOT PRINTED PRINTED PRINTED C C -1 NOT PRINTED PRINTED PRINTED C ONLY ONLY C ONCE ONCE C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 23 MAY 1979 C IF ((KONTRL.GE.(-2)).AND.(KONTRL.LE.2)) GO TO 10 CALL XERRWV(39HXSETF -- INVALID VALUE OF KONTRL (I1).,33,1,2, 1 1,KONTRL,0,0,0.,0.) RETURN 10 JUNK = J4SAVE(2,KONTRL,.TRUE.) RETURN END SUBROUTINE XSETUA(IUNIT,N) XER 10 C C ABSTRACT C XSETUA MAY BE CALLED TO DECLARE A LIST OF UP TO FIVE C LOGICAL UNITS, EACH OF WHICH IS TO RECEIVE A COPY OF C EACH ERROR MESSAGE PROCESSED BY THIS PACKAGE. C THE PURPOSE OF XSETUA IS TO ALLOW SIMULTANEOUS PRINTING C OF EACH ERROR MESSAGE ON, SAY, A MAIN OUTPUT FILE, C AN INTERACTIVE TERMINAL, AND OTHER FILES SUCH AS GRAPHICS C COMMUNICATION FILES. C C DESCRIPTION OF PARAMETERS C --INPUT-- C IUNIT - AN ARRAY OF UP TO FIVE UNIT NUMBERS. C NORMALLY THESE NUMBERS SHOULD ALL BE DIFFERENT. C (BUT DUPLICATES ARE NOT PROHIBITED.) C N - THE NUMBER OF UNIT NUMBERS PROVIDED IN IUNIT. C MUST HAVE 1 .LE. N .LE. 5. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 23 MAY 1979 C DIMENSION IUNIT(5) IF ((N.GE.1).AND.(N.LE.5)) GO TO 10 CALL XERRWV(34HXSETUA -- INVALID VALUE OF N (I1).,34,1,2, 1 1,N,0,0,0.,0.) RETURN 10 CONTINUE DO 20 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 JUNK = J4SAVE(INDEX,IUNIT(I),.TRUE.) 20 CONTINUE JUNK = J4SAVE(5,N,.TRUE.) RETURN END SUBROUTINE XSETUN(IUNIT) XER 10 C C ABSTRACT C XSETUN SETS THE OUTPUT FILE TO WHICH ERROR MESSAGES ARE TO C BE SENT. ONLY ONE FILE WILL BE USED. SEE XSETUA FOR C HOW TO DECLARE MORE THAN ONE FILE. C C DESCRIPTION OF PARAMETER C --INPUT-- C IUNIT - AN INPUT PARAMETER GIVING THE LOGICAL UNIT NUMBER C TO WHICH ERROR MESSAGES ARE TO BE SENT. C C WRITTEN BY RON JONES, WITH SLATEC COMMON MATH LIBRARY SUBCOMMITTEE C END OF ABSTRACT C LATEST REVISION --- 7 JUNE 1978 C JUNK = J4SAVE(3,IUNIT,.TRUE.) JUNK = J4SAVE(5,1,.TRUE.) RETURN END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***REVISION DATE 811015 (YYMMDD) C***CATEGORY NO. Q C***KEYWORDS MACHINE CONSTANTS,INTEGER C***DATE WRITTEN 1975 C***AUTHOR FOX P.A.,HALL A.D.,SCHRYER N.L. (BELL LABS) C***PURPOSE C RETURNS INTEGER MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C THESE MACHINE CONSTANT ROUTINES MUST BE ACTIVATED FOR C A PARTICULAR ENVIRONMENT. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C I1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C K = I1MACH(I) C C WHERE I=1,...,16. THE (OUTPUT) VALUE OF K ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C I/O UNIT NUMBERS. C I1MACH( 1) = THE STANDARD INPUT UNIT. C I1MACH( 2) = THE STANDARD OUTPUT UNIT. C I1MACH( 3) = THE STANDARD PUNCH UNIT. C I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT. C C WORDS. C I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT. C I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT. C C INTEGERS. C ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM C C SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1. C I1MACH( 7) = A, THE BASE. C I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS. C I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT, C BASE-B FORM C SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, C 0 .LT. X(1), AND EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, THE BASE. C C SINGLE-PRECISION C I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(13) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(16) = EMAX, THE LARGEST EXPONENT E. C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. C C***REFERENCES C FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A PORTABLE LIBRARY*, C ACM TRANSACTION ON MATHEMATICAL SOFTWARE, VOL. 4, NO. 2, C JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT C C EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) /6LOUTPUT/ C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA IMACH( 1) / 100 / C DATA IMACH( 2) / 101 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 101 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -8192 / C DATA IMACH(13) / 8191 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -8192 / C DATA IMACH(16) / 8191 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH(1) / 5/ C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 39 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 55 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN S SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN S SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 5 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024/ C DATA IMACH(16) / 1023 / C C C MACHINE CONSTANTS FOR THE VAX 11/780 C C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 5 / C DATA IMACH(4) / 6 / C DATA IMACH(5) / 32 / C DATA IMACH(6) / 4 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 31 / C DATA IMACH(9) /2147483647 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 24 / C DATA IMACH(12)/ -127 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 56 / C DATA IMACH(15)/ -127 / C DATA IMACH(16)/ 127 / C C***FIRST EXECUTABLE STATEMENT I1MACH C IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH=IMACH(I) RETURN C 10 CONTINUE WRITE(OUTPUT,9000) 9000 FORMAT(39H1ERROR 1 IN I1MACH - I OUT OF BOUNDS ) C C CALL FDUMP C C STOP END SUBROUTINE CLSTP(KLOG, COND, ISTAT) CLS 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 (BEGIN THE CHANGES AT THE LINE WITH C++ IN COLS. 1-3.) C /REAL (12 BLANKS)/DOUBLE PRECISION/,/SCOPY/DCOPY/,/SDOT/DDOT/, C /SNRM2/DNRM2/,/SQRT/DSQRT/,/E0/D0/,/SSCAL/DSCAL/,/SAXPY/DAXPY/, C /SRELPR/DRELPR/,/SSWAP/DSWAP/ C C REVISED 820305-2000 C REVISED YYMMDD-HHMM C C THIS SUBROUTINE EXERCISES MOST OF THE MATHEMATICAL FEATURES OF THE C CONSTRAINED LEAST SQUARES SUBPROGRAMS WNNLS( ) AND LSEI( ). C THE PROBLEM THAT IS SOLVED HERE IS OF THE FORM C C A*X=B (LEAST SQUARES, A MA BY N), C C C SUBJECT TO CONSTRAINTS C C E*X=F (CONSTRAINT EQUATIONS, ME BY N), C AND G*X .GE. H (INEQUALITY CONSTRAINTS, MG BY N). C C THE CLASS OF PROBLEMS THAT IS SOLVED HERE IS GENERATED WITH C HADAMARD MATRICES OF ORDER=POWER OF 2. EACH OF THE MATRICES C A,E, AND G HAVE A SPECIFIED CONDITION NUMBER. FOR EXAMPLE C A=HADAMARD MATRIX * DIAGONAL MATRIX * HADAMARD MATRIX. C DIAGONAL TERMS OF THE DIAGONAL MATRIX ARE CHOSEN SO THAT A C HAS A GIVEN CONDITION NUMBER. THE MATRICES E AND G ARE C CONSTRUCTED IN SIMILIAR WAYS. FURTHER, THE PROBLEM IS CONSTRUCTED C SO THAT THE TRUE SOLUTION IS X=(1,...,1) (TRANSPOSED). C THIS REQUIRES COMPUTING THE RIGHT HAND SIDE VECTORS B,F AND C H. THE VECTOR B=A*X+COMPONENT ORTHOGONAL TO COL. SPACE OF C A, F=E*X, AND H=G*H-SLACK COMPONENTS. THESE SLACK COMPONENTS C ARE CHOSEN SO THAT THE FIRST MI OF THE INEQUALITIES ARE C STRICT INEQUALITIES. C C THE PROBLEMS DIMENSIONS ARE SPECIFIED BY C C MA = 2**KA C ME = 2**KE C MG = 2**KG C MI = 2**KI C N = 2**KN C C WHERE KA, KE, KG, KI, AND KN ARE INPUT TO THE SUBROUTINE AS C DISCUSSED BELOW. C C THE SUBROUTINE ARGUMENTS ARE AS FOLLOWS C C I N P U T C C KLOG(*) - AN INTEGER ARRAY WHOSE DIMENSION IS AT LEAST 5. THE C ENTRIES CORRESPOND TO THE POWERS OF 2 NECESSARY C TO SPECIFY THE PROBLEM DIMENSIONS. REFERRING TO C THE ABOVE DISCUSSION, THE ENTRIES OF KLOG(*) C SHOULD BE SET AS FOLLOWS C C KLOG(1) = KA C KLOG(2) = KE C KLOG(3) = KG C KLOG(4) = KI C KLOG(5) = KN C C IF KA, KE, KG, OR KI IS LESS THAN ZERO, THE C CORRESPONDING DIMENSION WILL BE SET TO ZERO. C C KN.LT.0 WILL CAUSE THE SUBROUTINE TO SIMPLY RETURN. C C COND(*) - AN ARRAY WHOSE DIMENSION IS AT LEAST 3. THE C ENTRIES COND(1), COND(2), AND COND(3) RESPECTIVELY C SPECIFY THE CONDITION NUMBER FOR THE LEAST SQUARES C MATRIX, THE EQUALITY CONSTRAINT MATRIX, AND THE C INEQUALITY CONSTRAINT MATRIX. C C O U T P U T C C ISTAT - AN INTEGER FLAG WHICH INDICATES WHETHER THE SOLUTION C WAS CORRECTLY COMPUTED. C =1 NEITHER WNNLS( ) NOR LSEI( ) PASSED THE TEST. C =2 WNNLS( ) PASSED BUT LSEI( ) FAILED THE TEST. C =3 LSEI( ) PASSED BUT WNNLS( ) FAILED THE TEST. C =4 BOTH WNNLS( ) AND LSEI( ) PASSED THE TEST. C C THE DIMENSION STATEMENTS BELOW ARE SET UP TO SOLVE PROBLEMS FOR C WHICH NONE OF THE ABOVE LOGARITHMS IS GREATER THAN 5. TO CHANGE C THESE DIMENSIONS TO SOLVE LARGER PROBLEMS, USE THE FOLLOWING C FORMULAS C C DIMENSION W(MA+ME+MG,N+MG+1),X(N+MG),HH(MMAX,MMAX),GG(MMAX,MMAX) C DIMENSION WORK(2*(ME+N)+K+(MG+2)*(N+7)),IWORK(ME+MA+2*MG+N) C DIMENSION SA(MIN(MA,N)),SE(MIN(ME,N)),SG(MIN(MG,N)) C C WHERE MMAX = MAX(MA,ME,MG,N) C K = MAX(MA+MG,N) C C NOTE THAT IF THE DIMENSIONS ARE CHANGED, THE VALUES ASSIGNED TO C MDW, MDH, AND MDG BELOW MUST BE ALTERED APPROPRIATELY. THESE C ARE THE RESPECTIVE ROW DIMENSIONS OF THE ARRAYS W(*,*), HH(*,*), C AND GG(*,*). C++ REAL ANSR, BETA, BNORM, CONDA, CONDE, CONDG, * DXBYX, GAM, GNORM, ONE, PHI, RHO, RNORME, * RNORML, SOLERR, SRELPR, T, TWO, ZERO REAL SNRM2, SDOT, SQRT REAL W(96,65), X(64), HH(32,32), GG(32,32) REAL WORK(1518) INTEGER IWORK(0640) REAL SA(32), SE(32), SG(32) C C THE FOLLOWING DIMENSION STATEMENTS NEED NOT BE ALTERED TO C SOLVE LARGER PROBLEMS. REAL COND(3), PRGOPT(4) INTEGER KLOG(5) LOGICAL DONE C MDW = 96 MDH = 32 MDG = 32 ZERO = 0.E0 ONE = 1.E0 TWO = 2.E0 ISTAT = 1 C COMPUTE THE RELATIVE MACHINE PRECISION. SRELPR = ONE 10 IF (ONE+SRELPR.EQ.ONE) GO TO 20 SRELPR = SRELPR/TWO GO TO 10 20 SRELPR = SRELPR*TWO C C SET THE OUTPUT UNIT TO WHICH ERROR MESSAGES WILL BE PRINTED. LOUT = I1MACH(4) C C SET UP THE PROBLEM DIMENSIONS KA = KLOG(1) KE = KLOG(2) KG = KLOG(3) KI = KLOG(4) KN = KLOG(5) CONDA = ONE CONDE = ONE CONDG = ONE DONE = KN.LT.0 IF (.NOT.(DONE)) GO TO 30 RETURN 30 MA = 0 ME = 0 MG = 0 N = 0 C C SET NOISE-TO-SIGNAL RATIO PARAMETER FOR LEAST SQUARES EQUAS. C THIS ESSENTIALLY SPECIFIES THE RATIO C C NORM(RESIDUAL VECTOR) C --------------------- C NORM(A*X) ANSR = 0.01 C C SET UP THE CONDITION NUMBERS FOR THE MATRICES A, E, AND G. IF (KA.GE.0) CONDA = COND(1) IF (KE.GE.0) CONDE = COND(2) IF (KG.GE.0) CONDG = COND(3) C C CHECK THE VALIDITY OF THE INPUT IF (.NOT.(.NOT.(KA.LE.5 .AND. KE.LE.5 .AND. KG.LE.5 .AND. KI.LE.5 * .AND. KN.LE.5))) GO TO 40 WRITE (LOUT,99999) RETURN 40 IF (.NOT.(CONDA.LT.ONE .OR. CONDE.LT.ONE .OR. CONDG.LT.ONE)) GO * TO 50 WRITE (LOUT,99998) 99999 FORMAT (/42H KA, KE, KG, KI, AND KN MUST ALL BE .LE. 5, 7H AS REQ, * 53HUIRED BY THE CURRENT SUBPROGRAM DIMENSION STATEMENTS.) 99998 FORMAT (/46H CONDA, CONDE, AND CONDG MUST ALL BE .GE. ONE.) RETURN 50 CONTINUE ICASE = 1 60 CONTINUE ISEED = 100001 T = RAN(ISEED) ISEED = 0 C C COMPUTE THE PRE-MULTIPLYING HADAMARD MATRIX FOR E. K = KE ASSIGN 70 TO NGO GO TO 900 70 ME = NN C C SAVE THE HADAMARD MATRIX. J = 1 N20011 = ME GO TO 90 80 J = J + 1 90 IF ((N20011-J).LT.0) GO TO 100 CALL SCOPY(ME, HH(1,J), 1, GG(1,J), 1) GO TO 80 C C NOW FORM THE POST-MULTIPLYING HADAMARD MATRIX. 100 K = KN ASSIGN 110 TO NGO GO TO 900 110 N = NN C C COMPUTE THE SINGULAR VALUES OF THE MATRIX E. C DISTRIBUTE THEM UNIFORMLY BETWEEN 1. AND CONDE. SE(1) = CONDE MNE = MAX0(1,MIN0(ME,N)) SE(MNE) = ONE I = 2 N20016 = MNE - 1 GO TO 130 120 I = I + 1 130 IF ((N20016-I).LT.0) GO TO 140 SE(I) = ONE + RAN(ISEED)*(CONDE-ONE) GO TO 120 140 J = 1 N20020 = MNE GO TO 160 150 J = J + 1 160 IF ((N20020-J).LT.0) GO TO 170 CALL SSCAL(ME, SE(J), GG(1,J), 1) GO TO 150 170 J = 1 N20024 = N GO TO 190 180 J = J + 1 190 IF ((N20024-J).LT.0) GO TO 230 IF (ME.GT.0) W(1,J) = ZERO CALL SCOPY(ME, W(1,J), 0, W(1,J), 1) I = 1 N20028 = MNE GO TO 210 200 I = I + 1 210 IF ((N20028-I).LT.0) GO TO 220 CALL SAXPY(ME, HH(I,J), GG(1,I), 1, W(1,J), 1) GO TO 200 220 GO TO 180 C C COMPUTE E*X AND STORE IN W(*,N+1). 230 I = 1 N20032 = ME GO TO 250 240 I = I + 1 250 IF ((N20032-I).LT.0) GO TO 260 X(1) = ONE W(I,N+1) = SDOT(N,X,0,W(I,1),MDW) GO TO 240 C C COMPUTE THE PRE-MULTIPLYING HADAMARD MATRIX FOR A. 260 K = KA ASSIGN 270 TO NGO GO TO 900 270 MA = NN C C SAVE THE HADAMARD MATRIX. J = 1 N20037 = MA GO TO 290 280 J = J + 1 290 IF ((N20037-J).LT.0) GO TO 300 CALL SCOPY(MA, HH(1,J), 1, GG(1,J), 1) GO TO 280 C C NOW FORM THE POST-MULTIPLYING HADAMARD MATRIX. 300 K = KN ASSIGN 310 TO NGO GO TO 900 310 N = NN C C COMPUTE THE SINGULAR VALUES OF THE MATRIX A. C DISTRUBUTE THEM UNIFORMLY BETWEEN 1. AND CONDA. SA(1) = CONDA MNA = MAX0(1,MIN0(MA,N)) SA(MNA) = ONE I = 2 N20042 = MNA - 1 GO TO 330 320 I = I + 1 330 IF ((N20042-I).LT.0) GO TO 340 SA(I) = ONE + RAN(ISEED)*(CONDA-ONE) GO TO 320 340 J = 1 N20046 = MNA GO TO 360 350 J = J + 1 360 IF ((N20046-J).LT.0) GO TO 370 CALL SSCAL(MA, SA(J), GG(1,J), 1) GO TO 350 370 J = 1 N20050 = N GO TO 390 380 J = J + 1 390 IF ((N20050-J).LT.0) GO TO 430 C C COMPUTE THE PRODUCT IN PLACE INTO W(*,*). IF (MA.GT.0) W(ME+1,J) = ZERO CALL SCOPY(MA, W(ME+1,J), 0, W(ME+1,J), 1) I = 1 N20054 = MNA GO TO 410 400 I = I + 1 410 IF ((N20054-I).LT.0) GO TO 420 CALL SAXPY(MA, HH(I,J), GG(1,I), 1, W(ME+1,J), 1) GO TO 400 420 GO TO 380 C C COMPUTE A*X AND STORE IN W(*,N+1). 430 I = 1 N20058 = MA GO TO 450 440 I = I + 1 450 IF ((N20058-I).LT.0) GO TO 460 MEPI = ME + I X(1) = ONE W(MEPI,N+1) = SDOT(N,X,0,W(MEPI,1),MDW) GO TO 440 460 BNORM = SNRM2(MA,W(ME+1,N+1),1) C C ADD COMPONENTS TO RIGHT SIDE THAT ARE ORTHOGONAL TO COL. C SPACE OF A. K = KA ASSIGN 470 TO NGO GO TO 900 470 MA = NN I = N + 1 N20063 = MA GO TO 490 480 I = I + 1 490 IF ((N20063-I).LT.0) GO TO 500 T = RAN(ISEED)*BNORM*ANSR CALL SAXPY(MA, T, HH(1,I), 1, W(ME+1,N+1), 1) GO TO 480 C C COMPUTE THE PRE-MULTIPLYING HADAMARD MATRIX FOR G. 500 K = KG ASSIGN 510 TO NGO GO TO 900 510 MG = NN C C SAVE THE HADAMARD MATRIX. J = 1 N20068 = MG GO TO 530 520 J = J + 1 530 IF ((N20068-J).LT.0) GO TO 540 CALL SCOPY(MG, HH(1,J), 1, GG(1,J), 1) GO TO 520 C C NOW FORM THE POST-MULTIPLYING HADAMARD MATRIX. 540 K = KN ASSIGN 550 TO NGO GO TO 900 550 N = NN C C COMPUTE THE SINGULAR VALUES OF G. C DISTRIBUTE THEM UNIFORMLY BETWEEN 1. AND CONDG. SG(1) = CONDG MNG = MAX0(1,MIN0(MG,N)) SG(MNG) = ONE I = 2 N20073 = MNG - 1 GO TO 570 560 I = I + 1 570 IF ((N20073-I).LT.0) GO TO 580 SG(I) = ONE + RAN(ISEED)*(CONDG-ONE) GO TO 560 580 J = 1 N20077 = MNG GO TO 600 590 J = J + 1 600 IF ((N20077-J).LT.0) GO TO 610 CALL SSCAL(MG, SG(J), GG(1,J), 1) GO TO 590 610 J = 1 N20081 = N GO TO 630 620 J = J + 1 630 IF ((N20081-J).LT.0) GO TO 670 MEPMA = ME + MA IF (MG.GT.0) W(MEPMA+1,J) = ZERO CALL SCOPY(MG, W(MEPMA+1,J), 0, W(MEPMA+1,J), 1) I = 1 N20085 = MNG GO TO 650 640 I = I + 1 650 IF ((N20085-I).LT.0) GO TO 660 CALL SAXPY(MG, HH(I,J), GG(1,I), 1, W(MEPMA+1,J), 1) GO TO 640 660 GO TO 620 C C COMPUTE G*X AND STORE IN W(*,N+1). 670 I = 1 N20089 = MG GO TO 690 680 I = I + 1 690 IF ((N20089-I).LT.0) GO TO 700 IROW = I + MEPMA X(1) = ONE W(IROW,N+1) = SDOT(N,X,0,W(IROW,1),MDW) GO TO 680 C C MAKE FIRST MI=2**KI OF THE INEQUALITIES STRICT. 700 IF (.NOT.(KI.GE.0)) GO TO 710 MI = 1 GO TO 720 710 MI = 0 720 K = 1 N20096 = KI GO TO 740 730 K = K + 1 740 IF ((N20096-K).LT.0) GO TO 750 MI = MI + MI GO TO 730 750 GNORM = SNRM2(MG,W(MEPMA+1,N+1),1) I = 1 N20100 = MIN0(MI,MG) GO TO 770 760 I = I + 1 770 IF ((N20100-I).LT.0) GO TO 780 IROW = I + MEPMA W(IROW,N+1) = W(IROW,N+1) - RAN(ISEED)*GNORM GO TO 760 C C OBTAIN THE CONSTRAINED LEAST SQUARES SOLUTION. C C NOTE THE LENGTHS OF THE WORK ARRAYS IN IWORK(*). 780 IWORK(1) = 1518 IWORK(2) = 0640 IF (.NOT.(ICASE.EQ.1)) GO TO 810 C C EXCHANGE POSITIONS OF THE ROWS (A B) AND (G H). NP1 = N + 1 DO 790 J=1,NP1 IROW = ME + (MA+MG+2)/2 CALL SSWAP((MA+MG+1)/2, W(ME+1,J), 1, W(IROW,J), -1) 790 CONTINUE C C MOVE RT-SIDE TO W(*,N+MG+1). JCOL = N + MG + 1 CALL SCOPY(ME+MA+MG, W(1,N+1), 1, W(1,JCOL), 1) C C PUT IN SLACK VARIABLE COLS. AS REQUIRED. IF (.NOT.(MG.GT.0)) GO TO 800 W(1,N+1) = ZERO CALL SCOPY(MDW*MG, W(1,N+1), 0, W(1,N+1), 1) W(ME+1,N+1) = -ONE CALL SCOPY(MG, W(ME+1,N+1), 0, W(ME+1,N+1), MDW+1) 800 CONTINUE C C SET THE OPTION (NUMBER 6) FOR WNNLS( ) TO SCALE THE C COLUMNS OF THE MATRIX TO HAVE UNIT LENGTH. PRGOPT(1) = 4 PRGOPT(2) = 6 PRGOPT(3) = 1 PRGOPT(4) = 1 CALL WNNLS(W, MDW, ME+MG, MA, N+MG, N, PRGOPT, X, RNORML, MODE, * IWORK, WORK) GO TO 820 810 CONTINUE C C SET THE OPTION (NUMBER 2) FOR LSEI( ) TO SCALE THE C COLUMNS OF THE MATRIX TO HAVE UNIT LENGTH. PRGOPT(1) = 4 PRGOPT(2) = 2 PRGOPT(3) = 1 PRGOPT(4) = 1 CALL LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, MODE, * WORK, IWORK) 820 CONTINUE C C COMPUTE THE RESIDUAL SUM OF SQUARES OF ERRORS. SOLERR = ZERO I = 1 N20104 = N GO TO 840 830 I = I + 1 840 IF ((N20104-I).LT.0) GO TO 850 X(I) = ONE - X(I) SOLERR = SOLERR + X(I)**2 GO TO 830 850 SOLERR = SQRT(SOLERR) C C TEST SIZE OF ERROR (REF. LAWSON-HANSON, PAGE 51 AND CH. 16) PHI = 100. T = N DXBYX = SOLERR/SQRT(T) RHO = ONE IF (BNORM.NE.ZERO) RHO = RNORML/BNORM GAM = (6*MAX0(MA,N)-3*MIN0(MA,N))*MIN0(MA,N) BETA = CONDA*(ONE+CONDA*RHO)*GAM*PHI IF (.NOT.(DXBYX+BETA.EQ.BETA)) GO TO 860 ISTAT = ISTAT + ICASE GO TO 870 860 CONTINUE IF (MA.EQ.0 .AND. MODE.EQ.0) ISTAT = ISTAT + ICASE 870 CONTINUE IF (.NOT.(ICASE.EQ.1)) GO TO 880 WRITE (LOUT,99997) RNORML 99997 FORMAT (18H0LEAST SQS. RESID., E012.5, 2X, 12HFOR WNNLS( )/ * 30H ERRORS, 1-X(I), FOR WNNLS( ).) WRITE (LOUT,99996) (I,X(I),I=1,N) 99996 FORMAT (4(I4, E012.5)) GO TO 890 880 CONTINUE WRITE (LOUT,99995) RNORML, IWORK(1), IWORK(2) 99995 FORMAT (18H0LEAST SQS. RESID., E012.5, 2X, 11HFOR LSEI( )/ * 41H COMP. RANK OF E, COMP. RANK OF REDUCED A, 2I4/11H ERRORS, 1-, * 18HX(I), FOR LSEI( ).) WRITE (LOUT,99996) (I,X(I),I=1,N) 890 CONTINUE IF (ICASE.EQ.2) RETURN ICASE = 2 GO TO 60 C C PROCEDURE (GET HADAMARD MATRIX) 900 NN = 0 IF (.NOT.(K.GE.0)) GO TO 940 NN = 1 I = 1 N20114 = K GO TO 920 910 I = I + 1 920 IF ((N20114-I).LT.0) GO TO 930 NN = NN + NN GO TO 910 930 GO TO 950 940 NN = 0 GO TO 1080 C C PLACE THE SYMMETRIC HADAMARD MATRIX IN THE ARRAY HH(*,*). 950 HH(1,1) = ONE NN = 1 L = 1 N20118 = K GO TO 970 960 L = L + 1 970 IF ((N20118-L).LT.0) GO TO 1040 J = 1 N20122 = NN GO TO 990 980 J = J + 1 990 IF ((N20122-J).LT.0) GO TO 1000 CALL SCOPY(NN, HH(1,J), 1, HH(NN+1,J), 1) GO TO 980 1000 J = 1 N20126 = NN GO TO 1020 1010 J = J + 1 1020 IF ((N20126-J).LT.0) GO TO 1030 JPNN = J + NN CALL SCOPY(2*NN, HH(1,J), 1, HH(1,JPNN), 1) CALL SSCAL(NN, -ONE, HH(NN+1,JPNN), 1) GO TO 1010 1030 NN = NN + NN GO TO 960 C C MAKE THE MATRIX ORTHOGONAL BY SCALING THE ENTRIES. 1040 T = NN T = ONE/SQRT(T) J = 1 N20130 = NN GO TO 1060 1050 J = J + 1 1060 IF ((N20130-J).LT.0) GO TO 1070 CALL SSCAL(NN, T, HH(1,J), 1) GO TO 1050 1070 CONTINUE 1080 GO TO NGO, (70, 110, 270, 310, 470, 510, 550) END C PROGRAM CLSTST (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) MAN 10 C MAN 20 C TO CONVERT THE CONSTRAINED LEAST SQUARES PACKAGE OF SUBPROGRAMS MAN 30 C TO DOUBLE PRECISION, FOLLOW THE EDITING INSTRUCTIONS NEAR THE MAN 40 C BEGINNING OF THE FOLLOWING SUBPROGRAMS MAN 50 C LSEI(), LSI(), LPDP(), WNNLS(), WNLSM(), WNLIT(), MAN 60 C HFTI(), H12(), DIFF(), AND THE TESTING SUBPROGRAM CLSTP(). MAN 70 C REVISED 820305-2000 MAN 80 C REVISED YYMMDD-HHMM MAN 90 C MAN 100 C MAN 110 C CHANGE STRING /REAL /DOUBLE PRECISION/, AND MAN 120 C CHANGE STRING /E12/D12/ MAN 130 REAL COND(3) MAN 140 INTEGER KLOG(5) MAN 150 LOGICAL DONE MAN 160 C MAN 170 LIN = 5 MAN 180 LOUT = I1MACH(4) MAN 190 C MAN 200 C MAN 210 C READ IN THE LOGS OF THE VARIOUS DIMENSIONS. MAN 220 10 READ (LIN,99999) (KLOG(I),I=1,5) MAN 230 99999 FORMAT (5I5) MAN 240 DONE = KLOG(5).LT.0 MAN 250 IF (.NOT.(DONE)) GO TO 20 MAN 260 GO TO 30 MAN 270 C MAN 280 C READ THE CONDITION NUMBERS FOR THE THREE MATRICES. MAN 290 20 READ (LIN,99998) (COND(I),I=1,3) MAN 300 99998 FORMAT (3F10.0) MAN 310 C MAN 320 WRITE (LOUT,99997) KLOG, COND MAN 330 99997 FORMAT (/////17H KLOG(*) ARRAY = , 5I5/17H COND(*) ARRAY = , MAN 340 * 3E12.4) MAN 350 C MAN 360 C CALL TESTING SUBPROGRAM TO FORM CONSTRAINED LEAST SQUARES MAN 370 C SYSTEM AND SOLVE IT. MAN 380 CALL CLSTP(KLOG, COND, ISTAT) MAN 390 C MAN 400 WRITE (LOUT,99996) ISTAT MAN 410 99996 FORMAT (/23H FROM CLSTP( ), ISTAT =, I2, 2X, 16H(ISTAT=4 IS GOOD/ MAN 420 * 27X, 24H ISTAT=1,2,3 MAY BE BAD)) MAN 430 C MAN 440 GO TO 10 MAN 450 30 STOP MAN 460 END MAN 470 REAL FUNCTION RAN(K) RAN 10 C C RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 C BY PIKE AND HILL (MODIFIED BY HANSSON) C COLLECTED ALG. FROM CACM. C C THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH C FIXED POINT WORDLENGTH OF AT LEAST 29 BITS. IT IS C BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST C 29 BITS. C INTEGER IY,K DATA IY/100001/ C IF(K.GT.0) IY = K IY = IY * 125 IY = IY - (IY/2796203) * 2796203 RAN = FLOAT(IY) / 2796203.0E0 RETURN C ---------- LAST CARD OF RAN ---------- END 3 1 2 0 2 DATA 10 1. 1000. 1000. DATA 20 3 1 2 1 2 DATA 30 1. 1000. 1000. DATA 40 3 1 2 2 2 DATA 50 1000. 1000. 1000. DATA 60 4 2 2 2 3 DATA 70 1000. 1000. 1000. DATA 80 5 3 2 2 4 DATA 90 10000. 1000. 1000. DATA 100 0 0 0 0 -1 DATA 110