*DECK WNLIT SUBROUTINE WNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM, + IDOPE, DOPE, DONE) C***BEGIN PROLOGUE WNLIT C***SUBSIDIARY C***PURPOSE Subsidiary to WNNLS C***LIBRARY SLATEC C***TYPE SINGLE PRECISION (WNLIT-S, DWNLIT-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION C C This is a companion subprogram to WNNLS( ). C The documentation for WNNLS( ) has complete 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 Triangularize L1 by L1 subsystem, where L1=MIN(M,L), with C col interchanges. C C***SEE ALSO WNNLS C***ROUTINES CALLED H12, ISAMAX, SCOPY, SROTM, SROTMG, SSCAL, SSWAP, C WNLT1, WNLT2, WNLT3 C***REVISION HISTORY (YYMMDD) C 790701 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890618 Completely restructured and revised. (WRB & RWC) C 890620 Revised to make WNLT1, WNLT2, and WNLT3 subroutines. (RWC) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) C***END PROLOGUE WNLIT INTEGER IDOPE(*), IPIVOT(*), ITYPE(*), L, M, MDW, N REAL DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*) LOGICAL DONE C EXTERNAL H12, ISAMAX, SCOPY, SROTM, SROTMG, SSCAL, SSWAP, WNLT1, * WNLT2, WNLT3 INTEGER ISAMAX LOGICAL WNLT2 C REAL ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5), * T, TAU INTEGER I, I1, IMAX, IR, J, J1, JJ, JP, KRANK, L1, LB, LEND, ME, * MEND, NIV, NSOLN LOGICAL INDEP, RECALC C C***FIRST EXECUTABLE STATEMENT WNLIT ME = IDOPE(1) NSOLN = IDOPE(2) L1 = IDOPE(3) C ALSQ = DOPE(1) EANORM = DOPE(2) TAU = DOPE(3) C LB = MIN(M-1,L) RECALC = .TRUE. RNORM = 0.E0 KRANK = 0 C C We set FACTOR=1.0 so that the heavy weight ALAMDA will be C included in the test for column independence. C FACTOR = 1.E0 LEND = L DO 180 I=1,LB C C Set IR to point to the I-th row. C IR = I MEND = M CALL WNLT1 (I, LEND, M, IR, MDW, RECALC, IMAX, HBAR, H, SCALE, + W) C C Update column SS and find pivot column. C CALL WNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C C Perform column interchange. C Test independence of incoming column. C 130 IF (WNLT2(ME, MEND, IR, FACTOR, TAU, SCALE, W(1,I))) THEN C C Eliminate I-th column below diagonal using modified Givens C transformations applied to (A B). C C When operating near the ME line, use the largest element C above it as the pivot. C DO 160 J=M,I+1,-1 JP = J-1 IF (J.EQ.ME+1) THEN IMAX = ME AMAX = SCALE(ME)*W(ME,I)**2 DO 150 JP=J-1,I,-1 T = SCALE(JP)*W(JP,I)**2 IF (T.GT.AMAX) THEN IMAX = JP AMAX = T ENDIF 150 CONTINUE JP = IMAX ENDIF C IF (W(J,I).NE.0.E0) THEN CALL SROTMG (SCALE(JP), SCALE(J), W(JP,I), W(J,I), + SPARAM) W(J,I) = 0.E0 CALL SROTM (N+1-I, W(JP,I+1), MDW, W(J,I+1), MDW, + SPARAM) ENDIF 160 CONTINUE ELSE IF (LEND.GT.I) THEN C C Column I is dependent. Swap with column LEND. C Perform column interchange, C and find column in remaining set with largest SS. C CALL WNLT3 (I, LEND, M, MDW, IPIVOT, H, W) LEND = LEND - 1 IMAX = ISAMAX(LEND-I+1, H(I), 1) + I - 1 HBAR = H(IMAX) GO TO 130 ELSE KRANK = I - 1 GO TO 190 ENDIF 180 CONTINUE KRANK = L1 C 190 IF (KRANK.LT.ME) THEN FACTOR = ALSQ DO 200 I=KRANK+1,ME CALL SCOPY (L, 0.E0, 0, W(I,1), MDW) 200 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. C RECALC = .TRUE. LB = MIN(L+ME-KRANK, N) DO 270 I=L+1,LB IR = KRANK + I - L LEND = N MEND = ME CALL WNLT1 (I, LEND, ME, IR, MDW, RECALC, IMAX, HBAR, H, + SCALE, W) C C Update col ss and find pivot col C CALL WNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C C Perform column interchange C Eliminate elements in the I-th col. C DO 240 J=ME,IR+1,-1 IF (W(J,I).NE.0.E0) THEN CALL SROTMG (SCALE(J-1), SCALE(J), W(J-1,I), W(J,I), + SPARAM) W(J,I) = 0.E0 CALL SROTM (N+1-I, W(J-1,I+1), MDW,W(J,I+1), MDW, + SPARAM) ENDIF 240 CONTINUE C C I=column being eliminated. C Test independence of incoming column. C Remove any redundant or dependent equality constraints. C IF (.NOT.WNLT2(ME, MEND, IR, FACTOR,TAU,SCALE,W(1,I))) THEN JJ = IR DO 260 IR=JJ,ME CALL SCOPY (N, 0.E0, 0, W(IR,1), MDW) RNORM = RNORM + (SCALE(IR)*W(IR,N+1)/ALSQ)*W(IR,N+1) W(IR,N+1) = 0.E0 SCALE(IR) = 1.E0 C C Reclassify the zeroed row as a least squares equation. C ITYPE(IR) = 1 260 CONTINUE C C Reduce ME to reflect any discovered dependent equality C constraints. C ME = JJ - 1 GO TO 280 ENDIF 270 CONTINUE ENDIF 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(ME+1,I). C 280 IF (KRANK.LT.L1) THEN RECALC = .TRUE. C C Set FACTOR=ALSQ to remove effect of heavy weight from C test for column independence. C FACTOR = ALSQ DO 350 I=KRANK+1,L1 C C Set IR to point to the ME+1-st row. C IR = ME+1 LEND = L MEND = M CALL WNLT1 (I, L, M, IR, MDW, RECALC, IMAX, HBAR, H, SCALE, + W) C C Update column SS and find pivot column. C CALL WNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C C Perform column interchange. C Eliminate I-th column below the IR-th element. C DO 320 J=M,IR+1,-1 IF (W(J,I).NE.0.E0) THEN CALL SROTMG (SCALE(J-1), SCALE(J), W(J-1,I), W(J,I), + SPARAM) W(J,I) = 0.E0 CALL SROTM (N+1-I, W(J-1,I+1), MDW, W(J,I+1), MDW, + SPARAM) ENDIF 320 CONTINUE C C Test if new pivot element is near zero. C If so, the column is dependent. C Then check row norm test to be classified as independent. C T = SCALE(IR)*W(IR,I)**2 INDEP = T .GT. (TAU*EANORM)**2 IF (INDEP) THEN RN = 0.E0 DO 340 I1=IR,M DO 330 J1=I+1,N RN = MAX(RN, SCALE(I1)*W(I1,J1)**2) 330 CONTINUE 340 CONTINUE INDEP = T .GT. RN*TAU**2 ENDIF C C If independent, swap the IR-th and KRANK+1-th rows to C maintain the triangular form. Update the rank indicator C KRANK and the equality constraint pointer ME. C IF (.NOT.INDEP) GO TO 360 CALL SSWAP(N+1, W(KRANK+1,1), MDW, W(IR,1), MDW) CALL SSWAP(1, SCALE(KRANK+1), 1, SCALE(IR), 1) C C Reclassify the least square equation as an equality C constraint and rescale it. C ITYPE(IR) = 0 T = SQRT(SCALE(KRANK+1)) CALL SSCAL(N+1, T, W(KRANK+1,1), MDW) SCALE(KRANK+1) = ALSQ ME = ME+1 KRANK = KRANK+1 350 CONTINUE ENDIF C C If pseudorank is less than L, apply Householder transformation. C from right. C 360 IF (KRANK.LT.L) THEN DO 370 J=KRANK,1,-1 CALL H12 (1, J, KRANK+1, L, W(J,1), MDW, H(J), W, MDW, 1, + J-1) 370 CONTINUE ENDIF C NIV = KRANK + NSOLN - L IF (L.EQ.N) DONE = .TRUE. C C End of initial triangularization. C IDOPE(1) = ME IDOPE(2) = KRANK IDOPE(3) = NIV RETURN END