home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / SCHUD.FOR < prev    next >
Text File  |  1985-01-13  |  5KB  |  142 lines
  1.       SUBROUTINE SCHUD(R,LDR,P,X,Z,LDZ,NZ,Y,RHO,C,S)
  2.       INTEGER LDR,P,LDZ,NZ
  3.       REAL RHO(1),C(1)
  4.       REAL R(LDR,1),X(1),Z(LDZ,1),Y(1),S(1)
  5. C
  6. C     SCHUD UPDATES AN AUGMENTED CHOLESKY DECOMPOSITION OF THE
  7. C     TRIANGULAR PART OF AN AUGMENTED QR DECOMPOSITION.  SPECIFICALLY,
  8. C     GIVEN AN UPPER TRIANGULAR MATRIX R OF ORDER P, A ROW VECTOR
  9. C     X, A COLUMN VECTOR Z, AND A SCALAR Y, SCHUD DETERMINES A
  10. C     UNTIARY MATRIX U AND A SCALAR ZETA SUCH THAT
  11. C
  12. C
  13. C                              (R  Z)     (RR   ZZ )
  14. C                         U  * (    )  =  (        ) ,
  15. C                              (X  Y)     ( 0  ZETA)
  16. C
  17. C     WHERE RR IS UPPER TRIANGULAR.  IF R AND Z HAVE BEEN
  18. C     OBTAINED FROM THE FACTORIZATION OF A LEAST SQUARES
  19. C     PROBLEM, THEN RR AND ZZ ARE THE FACTORS CORRESPONDING TO
  20. C     THE PROBLEM WITH THE OBSERVATION (X,Y) APPENDED.  IN THIS
  21. C     CASE, IF RHO IS THE NORM OF THE RESIDUAL VECTOR, THEN THE
  22. C     NORM OF THE RESIDUAL VECTOR OF THE UPDATED PROBLEM IS
  23. C     SQRT(RHO**2 + ZETA**2).  SCHUD WILL SIMULTANEOUSLY UPDATE
  24. C     SEVERAL TRIPLETS (Z,Y,RHO).
  25. C     FOR A LESS TERSE DESCRIPTION OF WHAT SCHUD DOES AND HOW
  26. C     IT MAY BE APPLIED, SEE THE LINPACK GUIDE.
  27. C
  28. C     THE MATRIX U IS DETERMINED AS THE PRODUCT U(P)*...*U(1),
  29. C     WHERE U(I) IS A ROTATION IN THE (I,P+1) PLANE OF THE
  30. C     FORM
  31. C
  32. C                       (     C(I)      S(I) )
  33. C                       (                    ) .
  34. C                       (    -S(I)      C(I) )
  35. C
  36. C     THE ROTATIONS ARE CHOSEN SO THAT C(I) IS REAL.
  37. C
  38. C     ON ENTRY
  39. C
  40. C         R      REAL(LDR,P), WHERE LDR .GE. P.
  41. C                R CONTAINS THE UPPER TRIANGULAR MATRIX
  42. C                THAT IS TO BE UPDATED.  THE PART OF R
  43. C                BELOW THE DIAGONAL IS NOT REFERENCED.
  44. C
  45. C         LDR    INTEGER.
  46. C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
  47. C
  48. C         P      INTEGER.
  49. C                P IS THE ORDER OF THE MATRIX R.
  50. C
  51. C         X      REAL(P).
  52. C                X CONTAINS THE ROW TO BE ADDED TO R.  X IS
  53. C                NOT ALTERED BY SCHUD.
  54. C
  55. C         Z      REAL(LDZ,NZ), WHERE LDZ .GE. P.
  56. C                Z IS AN ARRAY CONTAINING NZ P-VECTORS TO
  57. C                BE UPDATED WITH R.
  58. C
  59. C         LDZ    INTEGER.
  60. C                LDZ IS THE LEADING DIMENSION OF THE ARRAY Z.
  61. C
  62. C         NZ     INTEGER.
  63. C                NZ IS THE NUMBER OF VECTORS TO BE UPDATED
  64. C                NZ MAY BE ZERO, IN WHICH CASE Z, Y, AND RHO
  65. C                ARE NOT REFERENCED.
  66. C
  67. C         Y      REAL(NZ).
  68. C                Y CONTAINS THE SCALARS FOR UPDATING THE VECTORS
  69. C                Z.  Y IS NOT ALTERED BY SCHUD.
  70. C
  71. C         RHO    REAL(NZ).
  72. C                RHO CONTAINS THE NORMS OF THE RESIDUAL
  73. C                VECTORS THAT ARE TO BE UPDATED.  IF RHO(J)
  74. C                IS NEGATIVE, IT IS LEFT UNALTERED.
  75. C
  76. C     ON RETURN
  77. C
  78. C         RC
  79. C         RHO    CONTAIN THE UPDATED QUANTITIES.
  80. C         Z
  81. C
  82. C         C      REAL(P).
  83. C                C CONTAINS THE COSINES OF THE TRANSFORMING
  84. C                ROTATIONS.
  85. C
  86. C         S      REAL(P).
  87. C                S CONTAINS THE SINES OF THE TRANSFORMING
  88. C                ROTATIONS.
  89. C
  90. C     LINPACK. THIS VERSION DATED 08/14/78 .
  91. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  92. C
  93. C     SCHUD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
  94. C
  95. C     EXTENDED BLAS SROTG
  96. C     FORTRAN SQRT
  97. C
  98.       INTEGER I,J,JM1
  99.       REAL AZETA,SCALE
  100.       REAL T,XJ,ZETA
  101. C
  102. C     UPDATE R.
  103. C
  104.       DO 30 J = 1, P
  105.          XJ = X(J)
  106. C
  107. C        APPLY THE PREVIOUS ROTATIONS.
  108. C
  109.          JM1 = J - 1
  110.          IF (JM1 .LT. 1) GO TO 20
  111.          DO 10 I = 1, JM1
  112.             T = C(I)*R(I,J) + S(I)*XJ
  113.             XJ = C(I)*XJ - S(I)*R(I,J)
  114.             R(I,J) = T
  115.    10    CONTINUE
  116.    20    CONTINUE
  117. C
  118. C        COMPUTE THE NEXT ROTATION.
  119. C
  120.          CALL SROTG(R(J,J),XJ,C(J),S(J))
  121.    30 CONTINUE
  122. C
  123. C     IF REQUIRED, UPDATE Z AND RHO.
  124. C
  125.       IF (NZ .LT. 1) GO TO 70
  126.       DO 60 J = 1, NZ
  127.          ZETA = Y(J)
  128.          DO 40 I = 1, P
  129.             T = C(I)*Z(I,J) + S(I)*ZETA
  130.             ZETA = C(I)*ZETA - S(I)*Z(I,J)
  131.             Z(I,J) = T
  132.    40    CONTINUE
  133.          AZETA = ABS(ZETA)
  134.          IF (AZETA .EQ. 0.0E0 .OR. RHO(J) .LT. 0.0E0) GO TO 50
  135.             SCALE = AZETA + RHO(J)
  136.             RHO(J) = SCALE*SQRT((AZETA/SCALE)**2+(RHO(J)/SCALE)**2)
  137.    50    CONTINUE
  138.    60 CONTINUE
  139.    70 CONTINUE
  140.       RETURN
  141.       END
  142.