home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / H12.FOR < prev    next >
Encoding:
Text File  |  1984-02-23  |  6.4 KB  |  81 lines
  1. C     SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)          H1200100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 H1200200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974H1200300
  4. C                                                                       H1200400
  5. C     CONSTRUCTION AND/OR APPLICATION OF A SINGLE                       H1200500
  6. C     HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B               H1200600
  7. C                                                                       H1200700
  8. C     MODE    = 1 OR 2   TO SELECT ALGORITHM  H1  OR  H2 .              H1200800
  9. C     LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.                         H1200900
  10. C     L1,M   IF L1 .LE. M   THE TRANSFORMATION WILL BE CONSTRUCTED TO   H1201000
  11. C            ZERO ELEMENTS INDEXED FROM L1 THROUGH M.   IF L1 GT. M     H1201100
  12. C            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.            H1201200
  13. C     U(),IUE,UP    ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR.       H1201300
  14. C                   IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS.      H1201400
  15. C                                       ON EXIT FROM H1 U() AND UP      H1201500
  16. C                   CONTAIN QUANTITIES DEFINING THE VECTOR U OF THE     H1201600
  17. C                   HOUSEHOLDER TRANSFORMATION.   ON ENTRY TO H2 U()    H1201700
  18. C                   AND UP SHOULD CONTAIN QUANTITIES PREVIOUSLY COMPUTEDH1201800
  19. C                   BY H1.  THESE WILL NOT BE MODIFIED BY H2.           H1201900
  20. C     C()    ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE   H1202000
  21. C            REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER      H1202100
  22. C            TRANSFORMATION IS TO BE APPLIED.  ON EXIT C() CONTAINS THE H1202200
  23. C            SET OF TRANSFORMED VECTORS.                                H1202300
  24. C     ICE    STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().      H1202400
  25. C     ICV    STORAGE INCREMENT BETWEEN VECTORS IN C().                  H1202500
  26. C     NCV    NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0  H1202600
  27. C            NO OPERATIONS WILL BE DONE ON C().                         H1202700
  28. C                                                                       H1202800
  29.       SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)          H1202900
  30.       DIMENSION U(IUE,M), C(1)                                          H1203000
  31.       DOUBLE PRECISION SM,B                                             H1203100
  32.       ONE=1.                                                            H1203200
  33. C                                                                       H1203300
  34.       IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN                H1203400
  35.       CL=ABS(U(1,LPIVOT))                                               H1203500
  36.       IF (MODE.EQ.2) GO TO 60                                           H1203600
  37. C                            ****** CONSTRUCT THE TRANSFORMATION. ******H1203700
  38.           DO 10 J=L1,M                                                  H1203800
  39.    10     CL=AMAX1(ABS(U(1,J)),CL)                                      H1203900
  40.       IF (CL) 130,130,20                                                H1204000
  41.    20 CLINV=ONE/CL                                                      H1204100
  42.       SM=(DBLE(U(1,LPIVOT))*CLINV)**2                                   H1204200
  43.           DO 30 J=L1,M                                                  H1204300
  44.    30     SM=SM+(DBLE(U(1,J))*CLINV)**2                                 H1204400
  45. C                              CONVERT DBLE. PREC. SM TO SNGL. PREC. SM1H1204500
  46.       SM1=SM                                                            H1204600
  47.       CL=CL*SQRT(SM1)                                                   H1204700
  48.       IF (U(1,LPIVOT)) 50,50,40                                         H1204800
  49.    40 CL=-CL                                                            H1204900
  50.    50 UP=U(1,LPIVOT)-CL                                                 H1205000
  51.       U(1,LPIVOT)=CL                                                    H1205100
  52.       GO TO 70                                                          H1205200
  53. C            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******H1205300
  54. C                                                                       H1205400
  55.    60 IF (CL) 130,130,70                                                H1205500
  56.    70 IF (NCV.LE.0) RETURN                                              H1205600
  57.       B=DBLE(UP)*U(1,LPIVOT)                                            H1205700
  58. C                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.H1205800
  59. C                                                                       H1205900
  60.       IF (B) 80,130,130                                                 H1206000
  61.    80 B=ONE/B                                                           H1206100
  62.       I2=1-ICV+ICE*(LPIVOT-1)                                           H1206200
  63.       INCR=ICE*(L1-LPIVOT)                                              H1206300
  64.           DO 120 J=1,NCV                                                H1206400
  65.           I2=I2+ICV                                                     H1206500
  66.           I3=I2+INCR                                                    H1206600
  67.           I4=I3                                                         H1206700
  68.           SM=C(I2)*DBLE(UP)                                             H1206800
  69.               DO 90 I=L1,M                                              H1206900
  70.               SM=SM+C(I3)*DBLE(U(1,I))                                  H1207000
  71.    90         I3=I3+ICE                                                 H1207100
  72.           IF (SM) 100,120,100                                           H1207200
  73.   100     SM=SM*B                                                       H1207300
  74.           C(I2)=C(I2)+SM*DBLE(UP)                                       H1207400
  75.               DO 110 I=L1,M                                             H1207500
  76.               C(I4)=C(I4)+SM*DBLE(U(1,I))                               H1207600
  77.   110         I4=I4+ICE                                                 H1207700
  78.   120     CONTINUE                                                      H1207800
  79.   130 RETURN                                                            H1207900
  80.       END                                                               H1208000
  81.