home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / SCHDD.FOR < prev    next >
Text File  |  1985-01-13  |  6KB  |  185 lines
  1.       SUBROUTINE SCHDD(R,LDR,P,X,Z,LDZ,NZ,Y,RHO,C,S,INFO)
  2.       INTEGER LDR,P,LDZ,NZ,INFO
  3.       REAL R(LDR,1),X(1),Z(LDZ,1),Y(1),S(1)
  4.       REAL RHO(1),C(1)
  5. C
  6. C     SCHDD DOWNDATES AN AUGMENTED CHOLESKY DECOMPOSITION OR THE
  7. C     TRIANGULAR FACTOR OF AN AUGMENTED QR DECOMPOSITION.
  8. C     SPECIFICALLY, GIVEN AN UPPER TRIANGULAR MATRIX R OF ORDER P,  A
  9. C     ROW VECTOR X, A COLUMN VECTOR Z, AND A SCALAR Y, SCHDD
  10. C     DETERMINEDS A ORTHOGONAL MATRIX U AND A SCALAR ZETA SUCH THAT
  11. C
  12. C                        (R   Z )     (RR  ZZ)
  13. C                    U * (      )  =  (      ) ,
  14. C                        (0 ZETA)     ( X   Y)
  15. C
  16. C     WHERE RR IS UPPER TRIANGULAR.  IF R AND Z HAVE BEEN OBTAINED
  17. C     FROM THE FACTORIZATION OF A LEAST SQUARES PROBLEM, THEN
  18. C     RR AND ZZ ARE THE FACTORS CORRESPONDING TO THE PROBLEM
  19. C     WITH THE OBSERVATION (X,Y) REMOVED.  IN THIS CASE, IF RHO
  20. C     IS THE NORM OF THE RESIDUAL VECTOR, THEN THE NORM OF
  21. C     THE RESIDUAL VECTOR OF THE DOWNDATED PROBLEM IS
  22. C     SQRT(RHO**2 - ZETA**2). SCHDD WILL SIMULTANEOUSLY DOWNDATE
  23. C     SEVERAL TRIPLETS (Z,Y,RHO) ALONG WITH R.
  24. C     FOR A LESS TERSE DESCRIPTION OF WHAT SCHDD DOES AND HOW
  25. C     IT MAY BE APPLIED, SEE THE LINPACK GUIDE.
  26. C
  27. C     THE MATRIX U IS DETERMINED AS THE PRODUCT U(1)*...*U(P)
  28. C     WHERE U(I) IS A ROTATION IN THE (P+1,I)-PLANE OF THE
  29. C     FORM
  30. C
  31. C                       ( C(I)     -S(I)     )
  32. C                       (                    ) .
  33. C                       ( S(I)       C(I)    )
  34. C
  35. C     THE ROTATIONS ARE CHOSEN SO THAT C(I) IS REAL.
  36. C
  37. C     THE USER IS WARNED THAT A GIVEN DOWNDATING PROBLEM MAY
  38. C     BE IMPOSSIBLE TO ACCOMPLISH OR MAY PRODUCE
  39. C     INACCURATE RESULTS.  FOR EXAMPLE, THIS CAN HAPPEN
  40. C     IF X IS NEAR A VECTOR WHOSE REMOVAL WILL REDUCE THE
  41. C     RANK OF R.  BEWARE.
  42. C
  43. C     ON ENTRY
  44. C
  45. C         R      REAL(LDR,P), WHERE LDR .GE. P.
  46. C                R CONTAINS THE UPPER TRIANGULAR MATRIX
  47. C                THAT IS TO BE DOWNDATED.  THE PART OF  R
  48. C                BELOW THE DIAGONAL IS NOT REFERENCED.
  49. C
  50. C         LDR    INTEGER.
  51. C                LDR IS THE LEADING DIMENSION FO THE ARRAY R.
  52. C
  53. C         P      INTEGER.
  54. C                P IS THE ORDER OF THE MATRIX R.
  55. C
  56. C         X      REAL(P).
  57. C                X CONTAINS THE ROW VECTOR THAT IS TO
  58. C                BE REMOVED FROM R.  X IS NOT ALTERED BY SCHDD.
  59. C
  60. C         Z      REAL(LDZ,NZ), WHERE LDZ .GE. P.
  61. C                Z IS AN ARRAY OF NZ P-VECTORS WHICH
  62. C                ARE TO BE DOWNDATED ALONG WITH R.
  63. C
  64. C         LDZ    INTEGER.
  65. C                LDZ IS THE LEADING DIMENSION OF THE ARRAY Z.
  66. C
  67. C         NZ     INTEGER.
  68. C                NZ IS THE NUMBER OF VECTORS TO BE DOWNDATED
  69. C                NZ MAY BE ZERO, IN WHICH CASE Z, Y, AND RHO
  70. C                ARE NOT REFERENCED.
  71. C
  72. C         Y      REAL(NZ).
  73. C                Y CONTAINS THE SCALARS FOR THE DOWNDATING
  74. C                OF THE VECTORS Z.  Y IS NOT ALTERED BY SCHDD.
  75. C
  76. C         RHO    REAL(NZ).
  77. C                RHO CONTAINS THE NORMS OF THE RESIDUAL
  78. C                VECTORS THAT ARE TO BE DOWNDATED.
  79. C
  80. C     ON RETURN
  81. C
  82. C         R
  83. C         Z      CONTAIN THE DOWNDATED QUANTITIES.
  84. C         RHO
  85. C
  86. C         C      REAL(P).
  87. C                C CONTAINS THE COSINES OF THE TRANSFORMING
  88. C                ROTATIONS.
  89. C
  90. C         S      REAL(P).
  91. C                S CONTAINS THE SINES OF THE TRANSFORMING
  92. C                ROTATIONS.
  93. C
  94. C         INFO   INTEGER.
  95. C                INFO IS SET AS FOLLOWS.
  96. C
  97. C                   INFO = 0  IF THE ENTIRE DOWNDATING
  98. C                             WAS SUCCESSFUL.
  99. C
  100. C                   INFO =-1  IF R COULD NOT BE DOWNDATED.
  101. C                             IN THIS CASE, ALL QUANTITIES
  102. C                             ARE LEFT UNALTERED.
  103. C
  104. C                   INFO = 1  IF SOME RHO COULD NOT BE
  105. C                             DOWNDATED.  THE OFFENDING RHOS ARE
  106. C                             SET TO -1.
  107. C
  108. C     LINPACK. THIS VERSION DATED 08/14/78 .
  109. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  110. C
  111. C     SCHDD USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  112. C
  113. C     FORTRAN ABS
  114. C     BLAS SDOT, SNRM2
  115. C
  116.       INTEGER I,II,J
  117.       REAL A,ALPHA,AZETA,NORM,SNRM2
  118.       REAL SDOT,T,ZETA,B,XX
  119. C
  120. C     SOLVE THE SYSTEM TRANS(R)*A = X, PLACING THE RESULT
  121. C     IN THE ARRAY S.
  122. C
  123.       INFO = 0
  124.       S(1) = X(1)/R(1,1)
  125.       IF (P .LT. 2) GO TO 20
  126.       DO 10 J = 2, P
  127.          S(J) = X(J) - SDOT(J-1,R(1,J),1,S,1)
  128.          S(J) = S(J)/R(J,J)
  129.    10 CONTINUE
  130.    20 CONTINUE
  131.       NORM = SNRM2(P,S,1)
  132.       IF (NORM .LT. 1.0E0) GO TO 30
  133.          INFO = -1
  134.       GO TO 120
  135.    30 CONTINUE
  136.          ALPHA = SQRT(1.0E0-NORM**2)
  137. C
  138. C        DETERMINE THE TRANSFORMATIONS.
  139. C
  140.          DO 40 II = 1, P
  141.             I = P - II + 1
  142.             SCALE = ALPHA + ABS(S(I))
  143.             A = ALPHA/SCALE
  144.             B = S(I)/SCALE
  145.             NORM = SQRT(A**2+B**2+0.0E0**2)
  146.             C(I) = A/NORM
  147.             S(I) = B/NORM
  148.             ALPHA = SCALE*NORM
  149.    40    CONTINUE
  150. C
  151. C        APPLY THE TRANSFORMATIONS TO R.
  152. C
  153.          DO 60 J = 1, P
  154.             XX = 0.0E0
  155.             DO 50 II = 1, J
  156.                I = J - II + 1
  157.                T = C(I)*XX + S(I)*R(I,J)
  158.                R(I,J) = C(I)*R(I,J) - S(I)*XX
  159.                XX = T
  160.    50       CONTINUE
  161.    60    CONTINUE
  162. C
  163. C        IF REQUIRED, DOWNDATE Z AND RHO.
  164. C
  165.          IF (NZ .LT. 1) GO TO 110
  166.          DO 100 J = 1, NZ
  167.             ZETA = Y(J)
  168.             DO 70 I = 1, P
  169.                Z(I,J) = (Z(I,J) - S(I)*ZETA)/C(I)
  170.                ZETA = C(I)*ZETA - S(I)*Z(I,J)
  171.    70       CONTINUE
  172.             AZETA = ABS(ZETA)
  173.             IF (AZETA .LE. RHO(J)) GO TO 80
  174.                INFO = 1
  175.                RHO(J) = -1.0E0
  176.             GO TO 90
  177.    80       CONTINUE
  178.                RHO(J) = RHO(J)*SQRT(1.0E0-(AZETA/RHO(J))**2)
  179.    90       CONTINUE
  180.   100    CONTINUE
  181.   110    CONTINUE
  182.   120 CONTINUE
  183.       RETURN
  184.       END
  185.