C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) H1200100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 H1200200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974H1200300 C H1200400 C CONSTRUCTION AND/OR APPLICATION OF A SINGLE H1200500 C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B H1200600 C H1200700 C MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 . H1200800 C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. H1200900 C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO H1201000 C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M H1201100 C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. H1201200 C U(),IUE,UP ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR. H1201300 C IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS. H1201400 C ON EXIT FROM H1 U() AND UP H1201500 C CONTAIN QUANTITIES DEFINING THE VECTOR U OF THE H1201600 C HOUSEHOLDER TRANSFORMATION. ON ENTRY TO H2 U() H1201700 C AND UP SHOULD CONTAIN QUANTITIES PREVIOUSLY COMPUTEDH1201800 C BY H1. THESE WILL NOT BE MODIFIED BY H2. H1201900 C C() ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE H1202000 C REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER H1202100 C TRANSFORMATION IS TO BE APPLIED. ON EXIT C() CONTAINS THE H1202200 C SET OF TRANSFORMED VECTORS. H1202300 C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). H1202400 C ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). H1202500 C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 H1202600 C NO OPERATIONS WILL BE DONE ON C(). H1202700 C H1202800 SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) H1202900 DIMENSION U(IUE,M), C(1) H1203000 DOUBLE PRECISION SM,B H1203100 ONE=1. H1203200 C H1203300 IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN H1203400 CL=ABS(U(1,LPIVOT)) H1203500 IF (MODE.EQ.2) GO TO 60 H1203600 C ****** CONSTRUCT THE TRANSFORMATION. ******H1203700 DO 10 J=L1,M H1203800 10 CL=AMAX1(ABS(U(1,J)),CL) H1203900 IF (CL) 130,130,20 H1204000 20 CLINV=ONE/CL H1204100 SM=(DBLE(U(1,LPIVOT))*CLINV)**2 H1204200 DO 30 J=L1,M H1204300 30 SM=SM+(DBLE(U(1,J))*CLINV)**2 H1204400 C CONVERT DBLE. PREC. SM TO SNGL. PREC. SM1H1204500 SM1=SM H1204600 CL=CL*SQRT(SM1) H1204700 IF (U(1,LPIVOT)) 50,50,40 H1204800 40 CL=-CL H1204900 50 UP=U(1,LPIVOT)-CL H1205000 U(1,LPIVOT)=CL H1205100 GO TO 70 H1205200 C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******H1205300 C H1205400 60 IF (CL) 130,130,70 H1205500 70 IF (NCV.LE.0) RETURN H1205600 B=DBLE(UP)*U(1,LPIVOT) H1205700 C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.H1205800 C H1205900 IF (B) 80,130,130 H1206000 80 B=ONE/B H1206100 I2=1-ICV+ICE*(LPIVOT-1) H1206200 INCR=ICE*(L1-LPIVOT) H1206300 DO 120 J=1,NCV H1206400 I2=I2+ICV H1206500 I3=I2+INCR H1206600 I4=I3 H1206700 SM=C(I2)*DBLE(UP) H1206800 DO 90 I=L1,M H1206900 SM=SM+C(I3)*DBLE(U(1,I)) H1207000 90 I3=I3+ICE H1207100 IF (SM) 100,120,100 H1207200 100 SM=SM*B H1207300 C(I2)=C(I2)+SM*DBLE(UP) H1207400 DO 110 I=L1,M H1207500 C(I4)=C(I4)+SM*DBLE(U(1,I)) H1207600 110 I4=I4+ICE H1207700 120 CONTINUE H1207800 130 RETURN H1207900 END H1208000