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 / WQRDC.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  6.7 KB  |  187 lines
  1.       SUBROUTINE WQRDC (XR, XI, LDX, N, P, QRAUXR, QRAUXI,
  2.      .                  JPVT, WORKR, WORKI, JOB)
  3.       IMPLICIT NONE
  4. C
  5. C USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR FACTORIZATION OF
  6. C  AN N BY P MATRIX X.  COLUMN PIVOTING BASED ON THE 2-NORMS OF THE REDUCED
  7. C  COLUMNS MAY BE PERFORMED AT THE USERS OPTION
  8. C
  9. C ON ENTRY:
  10. C    X     DOUBLE-COMPLEX(LDX,P), WHERE LDX .GE. N.   CONTAINS THE MATRIX
  11. C           WHOSE DECOMPOSITION IS TO BE COMPUTED
  12. C    LDX   INTEGER,  LEADING DIMENSION OF THE ARRAY X
  13. C    N     INTEGER,  NUMBER OF ROWS OF THE MATRIX X
  14. C    P     INTEGER,  NUMBER OF COLUMNS OF THE MATRIX X
  15. C    JPVT  INTEGER(P),  CONTAINS INTEGERS THAT CONTROL THE SELECTION OF THE
  16. C           PIVOT COLUMNS.  THE K-TH COLUMN X(K) OF X IS PLACED IN ONE OF
  17. C           THREE CLASSES ACCORDING TO THE VALUE OF JPVT(K):
  18. C             IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL COLUMN
  19. C             IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN
  20. C             IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN
  21. C           BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS ARE MOVED
  22. C           TO THE BEGINNING OF THE ARRAY X AND FINAL COLUMNS TO THE END.
  23. C           BOTH INITIAL AND FINAL COLUMNS ARE FROZEN IN PLACE DURING THE
  24. C           COMPUTATION AND ONLY FREE COLUMNS ARE MOVED.  AT THE K-TH STAGE
  25. C           OF THE REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN IT IS
  26. C           INTERCHANGED WITH THE FREE COLUMN OF LARGEST REDUCED NORM.
  27. C           JPVT IS NOT REFERENCED IF JOB .EQ. 0
  28. C    WORK  DOUBLE-COMPLEX(P),  WORK ARRAY.  NOT REFERENCED IF JOB .EQ. 0
  29. C    JOB   INTEGER,  INITIATES COLUMN PIVOTING:
  30. C           IF JOB .EQ. 0, NO PIVOTING IS DONE
  31. C           IF JOB .NE. 0, PIVOTING IS DONE
  32. C
  33. C ON RETURN:
  34. C    X      CONTAINS IN ITS UPPER TRIANGLE THE UPPER TRIANGULAR MATRIX R OF
  35. C            THE QR FACTORIZATION.  BELOW ITS DIAGONAL X CONTAINS INFORMATION
  36. C            FROM WHICH THE UNITARY PART OF THE DECOMPOSITION CAN BE RECOVERED.
  37. C            NOTE THAT IF PIVOTING HAS BEEN REQUESTED, THE DECOMPOSITION IS
  38. C            NOT THAT OF THE ORIGINAL MATRIX X BUT THAT OF X WITH ITS COLUMNS
  39. C            PERMUTED AS DESCRIBED BY JPVT
  40. C    QRAUX  DOUBLE-COMPLEX(P),  CONTAINS FURTHER INFORMATION REQUIRED TO
  41. C            RECOVER THE UNITARY PART OF THE DECOMPOSITION
  42. C    JPVT   JPVT(K) =  INDEX OF THE COLUMN OF THE ORIGINAL MATRIX THAT HAS
  43. C            BEEN INTERCHANGED INTO THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
  44. C
  45. C VERSION 07/03/79 C.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  46. C
  47.       INTEGER LDX, N, P, JOB, JPVT(*)
  48.       DOUBLE PRECISION XR(LDX,*), XI(LDX,*), QRAUXR(*),
  49.      .                 QRAUXI(*), WORKR(*), WORKI(*)
  50. C
  51.       INTEGER J, JP, L, LP1, LUP, MAXJ, PL, PU, JJ
  52.       DOUBLE PRECISION MAXNRM, TT, NRMXLR, NRMXLI, TR, TI
  53.       LOGICAL NEGJ, SWAPJ
  54.       DOUBLE PRECISION ZDUMR, ZDUMI
  55. C
  56.       DOUBLE PRECISION CABS1, PYTHAG, WDOTCR, WDOTCI, FLOP, WNRM2
  57. C
  58. C
  59.       CABS1 (ZDUMR, ZDUMI) = DABS (ZDUMR)+DABS (ZDUMI)
  60. C
  61. C
  62.       PL = 1
  63.       PU = 0
  64.       IF (JOB.EQ.0) GO TO 60
  65. C
  66. C ***      PIVOTING HAS BEEN REQUESTED.  REARRANGE COLUMNS ACCORDING TO JPVT
  67.       DO 20 J = 1, P
  68.         SWAPJ = JPVT(J).GT.0
  69.         NEGJ = JPVT(J).LT.0
  70.         JPVT(J) = J
  71.         IF (NEGJ) JPVT(J) = -J
  72.         IF (.NOT.SWAPJ) GO TO 10
  73.         IF (J.NE.PL) CALL WSWAP (N, XR(1,PL), XI(1,PL),
  74.      .                           1, XR(1,J), XI(1,J), 1)
  75.         JPVT(J) = JPVT(PL)
  76.         JPVT(PL) = J
  77.         PL = PL+1
  78. 10      CONTINUE
  79. 20    CONTINUE
  80.       PU = P
  81.       DO 50 JJ = 1, P
  82.         J = P-JJ+1
  83.         IF (JPVT(J).GE.0) GO TO 40
  84.         JPVT(J) = -JPVT(J)
  85.         IF (J.EQ.PU) GO TO 30
  86.         CALL WSWAP (N, XR(1,PU), XI(1,PU), 1, XR(1,J), XI(1,J), 1)
  87.         JP = JPVT(PU)
  88.         JPVT(PU) = JPVT(J)
  89.         JPVT(J) = JP
  90. 30      CONTINUE
  91.         PU = PU-1
  92. 40      CONTINUE
  93. 50    CONTINUE
  94. C
  95. C ***      COMPUTE THE NORMS OF THE FREE COLUMNS
  96. 60    CONTINUE
  97.       IF (PU.LT.PL) GO TO 80
  98.       DO 70 J = PL, PU
  99.         QRAUXR(J) = WNRM2 (N, XR(1,J), XI(1,J), 1)
  100.         QRAUXI(J) = 0.0D0
  101.         WORKR(J) = QRAUXR(J)
  102.         WORKI(J) = QRAUXI(J)
  103. 70    CONTINUE
  104. C
  105. C ***      PERFORM THE HOUSEHOLDER REDUCTION OF X
  106. 80    CONTINUE
  107.       LUP = MIN0 (N, P)
  108.       DO 210 L = 1, LUP
  109.         IF (L.LT.PL .OR. L.GE.PU) GO TO 120
  110. C
  111. C ***      LOCATE COLUMN OF LARGEST NORM AND BRING IT INTO THE PIVOT POSITION
  112.         MAXNRM = 0.0D0
  113.         MAXJ = L
  114.         DO 100 J = L, PU
  115.           IF (QRAUXR(J).LE.MAXNRM) GO TO 90
  116.           MAXNRM = QRAUXR(J)
  117.           MAXJ = J
  118. 90        CONTINUE
  119. 100     CONTINUE
  120.         IF (MAXJ.EQ.L) GO TO 110
  121.         CALL WSWAP (N, XR(1,L), XI(1,L), 1, XR(1,MAXJ), XI(1,MAXJ), 1)
  122.         QRAUXR(MAXJ) = QRAUXR(L)
  123.         QRAUXI(MAXJ) = QRAUXI(L)
  124.         WORKR(MAXJ) = WORKR(L)
  125.         WORKI(MAXJ) = WORKI(L)
  126.         JP = JPVT(MAXJ)
  127.         JPVT(MAXJ) = JPVT(L)
  128.         JPVT(L) = JP
  129. 110     CONTINUE
  130. 120     CONTINUE
  131.         QRAUXR(L) = 0.0D0
  132.         QRAUXI(L) = 0.0D0
  133.         IF (L.EQ.N) GO TO 200
  134. C
  135. C ***      COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L
  136.         NRMXLR = WNRM2 (N-L+1, XR(L,L), XI(L,L), 1)
  137.         NRMXLI = 0.0D0
  138.         IF (CABS1 (NRMXLR, NRMXLI).EQ.0.0D0) GO TO 190
  139.         IF (CABS1 (XR(L,L), XI(L,L)).EQ.0.0D0) GO TO 130
  140.         CALL WSIGN (NRMXLR, NRMXLI, XR(L,L), XI(L,L), NRMXLR, NRMXLI)
  141. 130     CONTINUE
  142.         CALL WDIV (1.0D0, 0.0D0, NRMXLR, NRMXLI, TR, TI)
  143.         CALL WSCAL (N-L+1, TR, TI, XR(L,L), XI(L,L), 1)
  144.         XR(L,L) = FLOP (1.0D0+XR(L,L))
  145. C
  146. C ***      APPLY TRANSFORMATION TO THE REMAINING COLUMNS, UPDATING THE NORMS
  147.         LP1 = L+1
  148.         IF (P.LT.LP1) GO TO 180
  149.         DO 170 J = LP1, P
  150.           TR = -WDOTCR (N-L+1, XR(L,L), XI(L,L), 1, XR(L,J), XI(L,J), 1)
  151.           TI = -WDOTCI (N-L+1, XR(L,L), XI(L,L), 1, XR(L,J), XI(L,J), 1)
  152.           CALL WDIV (TR, TI, XR(L,L), XI(L,L), TR, TI)
  153.           CALL WAXPY (N-L+1, TR, TI, XR(L,L), XI(L,L),
  154.      .                1, XR(L,J), XI(L,J), 1)
  155.           IF (J.LT.PL .OR. J.GT.PU) GO TO 160
  156.           IF (CABS1 (QRAUXR(J), QRAUXI(J)).EQ.0.0D0) GO TO 160
  157.           TT = 1.0D0-(PYTHAG (XR(L,J), XI(L,J))/QRAUXR(J))**2
  158.           TT = DMAX1 (TT, 0.0D0)
  159.           TR = FLOP (TT)
  160.           TT = FLOP (1.0D0+0.05D0*TT*(QRAUXR(J)/WORKR(J))**2)
  161.           IF (TT.EQ.1.0D0) GO TO 140
  162.           QRAUXR(J) = QRAUXR(J)*DSQRT (TR)
  163.           QRAUXI(J) = QRAUXI(J)*DSQRT (TR)
  164.           GO TO 150
  165. C
  166. 140       CONTINUE
  167.           QRAUXR(J) = WNRM2 (N-L, XR(L+1,J), XI(L+1,J), 1)
  168.           QRAUXI(J) = 0.0D0
  169.           WORKR(J) = QRAUXR(J)
  170.           WORKI(J) = QRAUXI(J)
  171. 150       CONTINUE
  172. 160       CONTINUE
  173. 170     CONTINUE
  174. C
  175. C ***      SAVE THE TRANSFORMATION
  176. 180     CONTINUE
  177.         QRAUXR(L) = XR(L,L)
  178.         QRAUXI(L) = XI(L,L)
  179.         XR(L,L) = -NRMXLR
  180.         XI(L,L) = -NRMXLI
  181. 190     CONTINUE
  182. 200     CONTINUE
  183. 210   CONTINUE
  184. C
  185.       RETURN
  186.       END
  187.