home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpklib.arc / SQRSL.FOR < prev    next >
Encoding:
Text File  |  1984-01-05  |  9.0 KB  |  274 lines
  1.       SUBROUTINE SQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO)
  2.       INTEGER LDX,N,K,JOB,INFO
  3.       REAL X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1)
  4. C
  5. C     SQRSL APPLIES THE OUTPUT OF SQRDC TO COMPUTE COORDINATE
  6. C     TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS.
  7. C     FOR K .LE. MIN(N,P), LET XK BE THE MATRIX
  8. C
  9. C            XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K)))
  10. C
  11. C     FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL
  12. C     N X P MATRIX X THAT WAS INPUT TO SQRDC (IF NO PIVOTING WAS
  13. C     DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR
  14. C     ORIGINAL ORDER).  SQRDC PRODUCES A FACTORED ORTHOGONAL MATRIX Q
  15. C     AND AN UPPER TRIANGULAR MATRIX R SUCH THAT
  16. C
  17. C              XK = Q * (R)
  18. C                       (0)
  19. C
  20. C     THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS
  21. C     X AND QRAUX.
  22. C
  23. C     ON ENTRY
  24. C
  25. C        X      REAL(LDX,P).
  26. C               X CONTAINS THE OUTPUT OF SQRDC.
  27. C
  28. C        LDX    INTEGER.
  29. C               LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  30. C
  31. C        N      INTEGER.
  32. C               N IS THE NUMBER OF ROWS OF THE MATRIX XK.  IT MUST
  33. C               HAVE THE SAME VALUE AS N IN SQRDC.
  34. C
  35. C        K      INTEGER.
  36. C               K IS THE NUMBER OF COLUMNS OF THE MATRIX XK.  K
  37. C               MUST NNOT BE GREATER THAN MIN(N,P), WHERE P IS THE
  38. C               SAME AS IN THE CALLING SEQUENCE TO SQRDC.
  39. C
  40. C        QRAUX  REAL(P).
  41. C               QRAUX CONTAINS THE AUXILIARY OUTPUT FROM SQRDC.
  42. C
  43. C        Y      REAL(N)
  44. C               Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED
  45. C               BY SQRSL.
  46. C
  47. C        JOB    INTEGER.
  48. C               JOB SPECIFIES WHAT IS TO BE COMPUTED.  JOB HAS
  49. C               THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING
  50. C               MEANING.
  51. C
  52. C                    IF A.NE.0, COMPUTE QY.
  53. C                    IF B,C,D, OR E .NE. 0, COMPUTE QTY.
  54. C                    IF C.NE.0, COMPUTE B.
  55. C                    IF D.NE.0, COMPUTE RSD.
  56. C                    IF E.NE.0, COMPUTE XB.
  57. C
  58. C               NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB
  59. C               AUTOMATICALLY TRIGGERS THE COMPUTATION OF QTY, FOR
  60. C               WHICH AN ARRAY MUST BE PROVIDED IN THE CALLING
  61. C               SEQUENCE.
  62. C
  63. C     ON RETURN
  64. C
  65. C        QY     REAL(N).
  66. C               QY CONNTAINS Q*Y, IF ITS COMPUTATION HAS BEEN
  67. C               REQUESTED.
  68. C
  69. C        QTY    REAL(N).
  70. C               QTY CONTAINS TRANS(Q)*Y, IF ITS COMPUTATION HAS
  71. C               BEEN REQUESTED.  HERE TRANS(Q) IS THE
  72. C               TRANSPOSE OF THE MATRIX Q.
  73. C
  74. C        B      REAL(K)
  75. C               B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM
  76. C
  77. C                    MINIMIZE NORM2(Y - XK*B),
  78. C
  79. C               IF ITS COMPUTATION HAS BEEN REQUESTED.  (NOTE THAT
  80. C               IF PIVOTING WAS REQUESTED IN SQRDC, THE J-TH
  81. C               COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J)
  82. C               OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO SQRDC.)
  83. C
  84. C        RSD    REAL(N).
  85. C               RSD CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B,
  86. C               IF ITS COMPUTATION HAS BEEN REQUESTED.  RSD IS
  87. C               ALSO THE ORTHOGONAL PROJECTION OF Y ONTO THE
  88. C               ORTHOGONAL COMPLEMENT OF THE COLUMN SPACE OF XK.
  89. C
  90. C        XB     REAL(N).
  91. C               XB CONTAINS THE LEAST SQUARES APPROXIMATION XK*B,
  92. C               IF ITS COMPUTATION HAS BEEN REQUESTED.  XB IS ALSO
  93. C               THE ORTHOGONAL PROJECTION OF Y ONTO THE COLUMN SPACE
  94. C               OF X.
  95. C
  96. C        INFO   INTEGER.
  97. C               INFO IS ZERO UNLESS THE COMPUTATION OF B HAS
  98. C               BEEN REQUESTED AND R IS EXACTLY SINGULAR.  IN
  99. C               THIS CASE, INFO IS THE INDEX OF THE FIRST ZERO
  100. C               DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED.
  101. C
  102. C     THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED
  103. C     IF THEIR COMPUTATION IS NOT REQUESTED AND IN THIS CASE
  104. C     CAN BE REPLACED BY DUMMY VARIABLES IN THE CALLING PROGRAM.
  105. C     TO SAVE STORAGE, THE USER MAY IN SOME CASES USE THE SAME
  106. C     ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE.  A
  107. C     FREQUENTLY OCCURING EXAMPLE IS WHEN ONE WISHES TO COMPUTE
  108. C     ANY OF B, RSD, OR XB AND DOES NOT NEED Y OR QTY.  IN THIS
  109. C     CASE ONE MAY IDENTIFY Y, QTY, AND ONE OF B, RSD, OR XB, WHILE
  110. C     PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE THAT IS TO BE
  111. C     COMPUTED.  THUS THE CALLING SEQUENCE
  112. C
  113. C          CALL SQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO)
  114. C
  115. C     WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD
  116. C     OVERWRITING Y.  MORE GENERALLY, EACH ITEM IN THE FOLLOWING
  117. C     LIST CONTAINS GROUPS OF PERMISSIBLE IDENTIFICATIONS FOR
  118. C     A SINGLE CALLINNG SEQUENCE.
  119. C
  120. C          1. (Y,QTY,B) (RSD) (XB) (QY)
  121. C
  122. C          2. (Y,QTY,RSD) (B) (XB) (QY)
  123. C
  124. C          3. (Y,QTY,XB) (B) (RSD) (QY)
  125. C
  126. C          4. (Y,QY) (QTY,B) (RSD) (XB)
  127. C
  128. C          5. (Y,QY) (QTY,RSD) (B) (XB)
  129. C
  130. C          6. (Y,QY) (QTY,XB) (B) (RSD)
  131. C
  132. C     IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO
  133. C     THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP.
  134. C
  135. C     LINPACK. THIS VERSION DATED 08/14/78 .
  136. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  137. C
  138. C     SQRSL USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  139. C
  140. C     BLAS SAXPY,SCOPY,SDOT
  141. C     FORTRAN ABS,MIN0,MOD
  142. C
  143. C     INTERNAL VARIABLES
  144. C
  145.       INTEGER I,J,JJ,JU,KP1
  146.       REAL SDOT,T,TEMP
  147.       LOGICAL CB,CQY,CQTY,CR,CXB
  148. C
  149. C
  150. C     SET INFO FLAG.
  151. C
  152.       INFO = 0
  153. C
  154. C     DETERMINE WHAT IS TO BE COMPUTED.
  155. C
  156.       CQY = JOB/10000 .NE. 0
  157.       CQTY = MOD(JOB,10000) .NE. 0
  158.       CB = MOD(JOB,1000)/100 .NE. 0
  159.       CR = MOD(JOB,100)/10 .NE. 0
  160.       CXB = MOD(JOB,10) .NE. 0
  161.       JU = MIN0(K,N-1)
  162. C
  163. C     SPECIAL ACTION WHEN N=1.
  164. C
  165.       IF (JU .NE. 0) GO TO 40
  166.          IF (CQY) QY(1) = Y(1)
  167.          IF (CQTY) QTY(1) = Y(1)
  168.          IF (CXB) XB(1) = Y(1)
  169.          IF (.NOT.CB) GO TO 30
  170.             IF (X(1,1) .NE. 0.0E0) GO TO 10
  171.                INFO = 1
  172.             GO TO 20
  173.    10       CONTINUE
  174.                B(1) = Y(1)/X(1,1)
  175.    20       CONTINUE
  176.    30    CONTINUE
  177.          IF (CR) RSD(1) = 0.0E0
  178.       GO TO 250
  179.    40 CONTINUE
  180. C
  181. C        SET UP TO COMPUTE QY OR QTY.
  182. C
  183.          IF (CQY) CALL SCOPY(N,Y,1,QY,1)
  184.          IF (CQTY) CALL SCOPY(N,Y,1,QTY,1)
  185.          IF (.NOT.CQY) GO TO 70
  186. C
  187. C           COMPUTE QY.
  188. C
  189.             DO 60 JJ = 1, JU
  190.                J = JU - JJ + 1
  191.                IF (QRAUX(J) .EQ. 0.0E0) GO TO 50
  192.                   TEMP = X(J,J)
  193.                   X(J,J) = QRAUX(J)
  194.                   T = -SDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J)
  195.                   CALL SAXPY(N-J+1,T,X(J,J),1,QY(J),1)
  196.                   X(J,J) = TEMP
  197.    50          CONTINUE
  198.    60       CONTINUE
  199.    70    CONTINUE
  200.          IF (.NOT.CQTY) GO TO 100
  201. C
  202. C           COMPUTE TRANS(Q)*Y.
  203. C
  204.             DO 90 J = 1, JU
  205.                IF (QRAUX(J) .EQ. 0.0E0) GO TO 80
  206.                   TEMP = X(J,J)
  207.                   X(J,J) = QRAUX(J)
  208.                   T = -SDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J)
  209.                   CALL SAXPY(N-J+1,T,X(J,J),1,QTY(J),1)
  210.                   X(J,J) = TEMP
  211.    80          CONTINUE
  212.    90       CONTINUE
  213.   100    CONTINUE
  214. C
  215. C        SET UP TO COMPUTE B, RSD, OR XB.
  216. C
  217.          IF (CB) CALL SCOPY(K,QTY,1,B,1)
  218.          KP1 = K + 1
  219.          IF (CXB) CALL SCOPY(K,QTY,1,XB,1)
  220.          IF (CR .AND. K .LT. N) CALL SCOPY(N-K,QTY(KP1),1,RSD(KP1),1)
  221.          IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120
  222.             DO 110 I = KP1, N
  223.                XB(I) = 0.0E0
  224.   110       CONTINUE
  225.   120    CONTINUE
  226.          IF (.NOT.CR) GO TO 140
  227.             DO 130 I = 1, K
  228.                RSD(I) = 0.0E0
  229.   130       CONTINUE
  230.   140    CONTINUE
  231.          IF (.NOT.CB) GO TO 190
  232. C
  233. C           COMPUTE B.
  234. C
  235.             DO 170 JJ = 1, K
  236.                J = K - JJ + 1
  237.                IF (X(J,J) .NE. 0.0E0) GO TO 150
  238.                   INFO = J
  239. C           ......EXIT
  240.                   GO TO 180
  241.   150          CONTINUE
  242.                B(J) = B(J)/X(J,J)
  243.                IF (J .EQ. 1) GO TO 160
  244.                   T = -B(J)
  245.                   CALL SAXPY(J-1,T,X(1,J),1,B,1)
  246.   160          CONTINUE
  247.   170       CONTINUE
  248.   180       CONTINUE
  249.   190    CONTINUE
  250.          IF (.NOT.CR .AND. .NOT.CXB) GO TO 240
  251. C
  252. C           COMPUTE RSD OR XB AS REQUIRED.
  253. C
  254.             DO 230 JJ = 1, JU
  255.                J = JU - JJ + 1
  256.                IF (QRAUX(J) .EQ. 0.0E0) GO TO 220
  257.                   TEMP = X(J,J)
  258.                   X(J,J) = QRAUX(J)
  259.                   IF (.NOT.CR) GO TO 200
  260.                      T = -SDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J)
  261.                      CALL SAXPY(N-J+1,T,X(J,J),1,RSD(J),1)
  262.   200             CONTINUE
  263.                   IF (.NOT.CXB) GO TO 210
  264.                      T = -SDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J)
  265.                      CALL SAXPY(N-J+1,T,X(J,J),1,XB(J),1)
  266.   210             CONTINUE
  267.                   X(J,J) = TEMP
  268.   220          CONTINUE
  269.   230       CONTINUE
  270.   240    CONTINUE
  271.   250 CONTINUE
  272.       RETURN
  273.       END
  274.