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