home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / fortran / linpklib.arc / SQRDC.FOR < prev    next >
Text File  |  1984-01-06  |  7KB  |  208 lines
  1.       SUBROUTINE SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
  2.       INTEGER LDX,N,P,JOB
  3.       INTEGER JPVT(1)
  4.       REAL X(LDX,1),QRAUX(1),WORK(1)
  5. C
  6. C     SQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR
  7. C     FACTORIZATION OF AN N BY P MATRIX X.  COLUMN PIVOTING
  8. C     BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE
  9. C     PERFORMED AT THE USERS OPTION.
  10. C
  11. C     ON ENTRY
  12. C
  13. C        X       REAL(LDX,P), WHERE LDX .GE. N.
  14. C                X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE
  15. C                COMPUTED.
  16. C
  17. C        LDX     INTEGER.
  18. C                LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  19. C
  20. C        N       INTEGER.
  21. C                N IS THE NUMBER OF ROWS OF THE MATRIX X.
  22. C
  23. C        P       INTEGER.
  24. C                P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  25. C
  26. C        JPVT    INTEGER(P).
  27. C                JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION
  28. C                OF THE PIVOT COLUMNS.  THE K-TH COLUMN X(K) OF X
  29. C                IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE
  30. C                VALUE OF JPVT(K).
  31. C
  32. C                   IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL
  33. C                                      COLUMN.
  34. C
  35. C                   IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN.
  36. C
  37. C                   IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN.
  38. C
  39. C                BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS
  40. C                ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL
  41. C                COLUMNS TO THE END.  BOTH INITIAL AND FINAL COLUMNS
  42. C                ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY
  43. C                FREE COLUMNS ARE MOVED.  AT THE K-TH STAGE OF THE
  44. C                REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN
  45. C                IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST
  46. C                REDUCED NORM.  JPVT IS NOT REFERENCED IF
  47. C                JOB .EQ. 0.
  48. C
  49. C        WORK    REAL(P).
  50. C                WORK IS A WORK ARRAY.  WORK IS NOT REFERENCED IF
  51. C                JOB .EQ. 0.
  52. C
  53. C        JOB     INTEGER.
  54. C                JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING.
  55. C                IF JOB .EQ. 0, NO PIVOTING IS DONE.
  56. C                IF JOB .NE. 0, PIVOTING IS DONE.
  57. C
  58. C     ON RETURN
  59. C
  60. C        X       X CONTAINS IN ITS UPPER TRIANGLE THE UPPER
  61. C                TRIANGULAR MATRIX R OF THE QR FACTORIZATION.
  62. C                BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM
  63. C                WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION
  64. C                CAN BE RECOVERED.  NOTE THAT IF PIVOTING HAS
  65. C                BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT
  66. C                OF THE ORIGINAL MATRIX X BUT THAT OF X
  67. C                WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT.
  68. C
  69. C        QRAUX   REAL(P).
  70. C                QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER
  71. C                THE ORTHOGONAL PART OF THE DECOMPOSITION.
  72. C
  73. C        JPVT    JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE
  74. C                ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO
  75. C                THE K-TH COLUMN, IF PIVOTING WAS REQUESTED.
  76. C
  77. C     LINPACK. THIS VERSION DATED 08/14/78 .
  78. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  79. C
  80. C     SQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  81. C
  82. C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2
  83. C     FORTRAN ABS,AMAX1,MIN0,SQRT
  84. C
  85. C     INTERNAL VARIABLES
  86. C
  87.       INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
  88.       REAL MAXNRM,SNRM2,TT
  89.       REAL SDOT,NRMXL,T
  90.       LOGICAL NEGJ,SWAPJ
  91. C
  92. C
  93.       PL = 1
  94.       PU = 0
  95.       IF (JOB .EQ. 0) GO TO 60
  96. C
  97. C        PIVOTING HAS BEEN REQUESTED.  REARRANGE THE COLUMNS
  98. C        ACCORDING TO JPVT.
  99. C
  100.          DO 20 J = 1, P
  101.             SWAPJ = JPVT(J) .GT. 0
  102.             NEGJ = JPVT(J) .LT. 0
  103.             JPVT(J) = J
  104.             IF (NEGJ) JPVT(J) = -J
  105.             IF (.NOT.SWAPJ) GO TO 10
  106.                IF (J .NE. PL) CALL SSWAP(N,X(1,PL),1,X(1,J),1)
  107.                JPVT(J) = JPVT(PL)
  108.                JPVT(PL) = J
  109.                PL = PL + 1
  110.    10       CONTINUE
  111.    20    CONTINUE
  112.          PU = P
  113.          DO 50 JJ = 1, P
  114.             J = P - JJ + 1
  115.             IF (JPVT(J) .GE. 0) GO TO 40
  116.                JPVT(J) = -JPVT(J)
  117.                IF (J .EQ. PU) GO TO 30
  118.                   CALL SSWAP(N,X(1,PU),1,X(1,J),1)
  119.                   JP = JPVT(PU)
  120.                   JPVT(PU) = JPVT(J)
  121.                   JPVT(J) = JP
  122.    30          CONTINUE
  123.                PU = PU - 1
  124.    40       CONTINUE
  125.    50    CONTINUE
  126.    60 CONTINUE
  127. C
  128. C     COMPUTE THE NORMS OF THE FREE COLUMNS.
  129. C
  130.       IF (PU .LT. PL) GO TO 80
  131.       DO 70 J = PL, PU
  132.          QRAUX(J) = SNRM2(N,X(1,J),1)
  133.          WORK(J) = QRAUX(J)
  134.    70 CONTINUE
  135.    80 CONTINUE
  136. C
  137. C     PERFORM THE HOUSEHOLDER REDUCTION OF X.
  138. C
  139.       LUP = MIN0(N,P)
  140.       DO 200 L = 1, LUP
  141.          IF (L .LT. PL .OR. L .GE. PU) GO TO 120
  142. C
  143. C           LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
  144. C           INTO THE PIVOT POSITION.
  145. C
  146.             MAXNRM = 0.0E0
  147.             MAXJ = L
  148.             DO 100 J = L, PU
  149.                IF (QRAUX(J) .LE. MAXNRM) GO TO 90
  150.                   MAXNRM = QRAUX(J)
  151.                   MAXJ = J
  152.    90          CONTINUE
  153.   100       CONTINUE
  154.             IF (MAXJ .EQ. L) GO TO 110
  155.                CALL SSWAP(N,X(1,L),1,X(1,MAXJ),1)
  156.                QRAUX(MAXJ) = QRAUX(L)
  157.                WORK(MAXJ) = WORK(L)
  158.                JP = JPVT(MAXJ)
  159.                JPVT(MAXJ) = JPVT(L)
  160.                JPVT(L) = JP
  161.   110       CONTINUE
  162.   120    CONTINUE
  163.          QRAUX(L) = 0.0E0
  164.          IF (L .EQ. N) GO TO 190
  165. C
  166. C           COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
  167. C
  168.             NRMXL = SNRM2(N-L+1,X(L,L),1)
  169.             IF (NRMXL .EQ. 0.0E0) GO TO 180
  170.                IF (X(L,L) .NE. 0.0E0) NRMXL = SIGN(NRMXL,X(L,L))
  171.                CALL SSCAL(N-L+1,1.0E0/NRMXL,X(L,L),1)
  172.                X(L,L) = 1.0E0 + X(L,L)
  173. C
  174. C              APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
  175. C              UPDATING THE NORMS.
  176. C
  177.                LP1 = L + 1
  178.                IF (P .LT. LP1) GO TO 170
  179.                DO 160 J = LP1, P
  180.                   T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
  181.                   CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
  182.                   IF (J .LT. PL .OR. J .GT. PU) GO TO 150
  183.                   IF (QRAUX(J) .EQ. 0.0E0) GO TO 150
  184.                      TT = 1.0E0 - (ABS(X(L,J))/QRAUX(J))**2
  185.                      TT = AMAX1(TT,0.0E0)
  186.                      T = TT
  187.                      TT = 1.0E0 + 0.05E0*TT*(QRAUX(J)/WORK(J))**2
  188.                      IF (TT .EQ. 1.0E0) GO TO 130
  189.                         QRAUX(J) = QRAUX(J)*SQRT(T)
  190.                      GO TO 140
  191.   130                CONTINUE
  192.                         QRAUX(J) = SNRM2(N-L,X(L+1,J),1)
  193.                         WORK(J) = QRAUX(J)
  194.   140                CONTINUE
  195.   150             CONTINUE
  196.   160          CONTINUE
  197.   170          CONTINUE
  198. C
  199. C              SAVE THE TRANSFORMATION.
  200. C
  201.                QRAUX(L) = X(L,L)
  202.                X(L,L) = -NRMXL
  203.   180       CONTINUE
  204.   190    CONTINUE
  205.   200 CONTINUE
  206.       RETURN
  207.       END
  208.