← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/source_f/slatec/qrfac.f

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/source_f/slatec/qrfac.f
 
1
*DECK QRFAC
 
2
      SUBROUTINE QRFAC (M, N, A, LDA, PIVOT, IPVT, LIPVT, SIGMA, ACNORM,
 
3
     +   WA)
 
4
C***BEGIN PROLOGUE  QRFAC
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to SNLS1, SNLS1E, SNSQ and SNSQE
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      SINGLE PRECISION (QRFAC-S, DQRFAC-D)
 
9
C***AUTHOR  (UNKNOWN)
 
10
C***DESCRIPTION
 
11
C
 
12
C     This subroutine uses Householder transformations with column
 
13
C     pivoting (optional) to compute a QR factorization of the
 
14
C     M by N matrix A. That is, QRFAC determines an orthogonal
 
15
C     matrix Q, a permutation matrix P, and an upper trapezoidal
 
16
C     matrix R with diagonal elements of nonincreasing magnitude,
 
17
C     such that A*P = Q*R. The Householder transformation for
 
18
C     column K, K = 1,2,...,MIN(M,N), is of the form
 
19
C
 
20
C                           T
 
21
C           I - (1/U(K))*U*U
 
22
C
 
23
C     where U has zeros in the first K-1 positions. The form of
 
24
C     this transformation and the method of pivoting first
 
25
C     appeared in the corresponding LINPACK subroutine.
 
26
C
 
27
C     The subroutine statement is
 
28
C
 
29
C       SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
 
30
C
 
31
C     where
 
32
C
 
33
C       M is a positive integer input variable set to the number
 
34
C         of rows of A.
 
35
C
 
36
C       N is a positive integer input variable set to the number
 
37
C         of columns of A.
 
38
C
 
39
C       A is an M by N array. On input A contains the matrix for
 
40
C         which the QR factorization is to be computed. On output
 
41
C         the strict upper trapezoidal part of A contains the strict
 
42
C         upper trapezoidal part of R, and the lower trapezoidal
 
43
C         part of A contains a factored form of Q (the non-trivial
 
44
C         elements of the U vectors described above).
 
45
C
 
46
C       LDA is a positive integer input variable not less than M
 
47
C         which specifies the leading dimension of the array A.
 
48
C
 
49
C       PIVOT is a logical input variable. If pivot is set .TRUE.,
 
50
C         then column pivoting is enforced. If pivot is set .FALSE.,
 
51
C         then no column pivoting is done.
 
52
C
 
53
C       IPVT is an integer output array of length LIPVT. IPVT
 
54
C         defines the permutation matrix P such that A*P = Q*R.
 
55
C         Column J of P is column IPVT(J) of the identity matrix.
 
56
C         If pivot is .FALSE., IPVT is not referenced.
 
57
C
 
58
C       LIPVT is a positive integer input variable. If PIVOT is
 
59
C             .FALSE., then LIPVT may be as small as 1. If PIVOT is
 
60
C             .TRUE., then LIPVT must be at least N.
 
61
C
 
62
C       SIGMA is an output array of length N which contains the
 
63
C         diagonal elements of R.
 
64
C
 
65
C       ACNORM is an output array of length N which contains the
 
66
C         norms of the corresponding columns of the input matrix A.
 
67
C         If this information is not needed, then ACNORM can coincide
 
68
C         with SIGMA.
 
69
C
 
70
C       WA is a work array of length N. If pivot is .FALSE., then WA
 
71
C         can coincide with SIGMA.
 
72
C
 
73
C***SEE ALSO  SNLS1, SNLS1E, SNSQ, SNSQE
 
74
C***ROUTINES CALLED  ENORM, R1MACH
 
75
C***REVISION HISTORY  (YYMMDD)
 
76
C   800301  DATE WRITTEN
 
77
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
78
C   890831  Modified array declarations.  (WRB)
 
79
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
80
C   900326  Removed duplicate information from DESCRIPTION section.
 
81
C           (WRB)
 
82
C   900328  Added TYPE section.  (WRB)
 
83
C***END PROLOGUE  QRFAC
 
84
      INTEGER M,N,LDA,LIPVT
 
85
      INTEGER IPVT(*)
 
86
      LOGICAL PIVOT
 
87
      REAL A(LDA,*),SIGMA(*),ACNORM(*),WA(*)
 
88
      INTEGER I,J,JP1,K,KMAX,MINMN
 
89
      REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
 
90
      REAL R1MACH,ENORM
 
91
      SAVE ONE, P05, ZERO
 
92
      DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/
 
93
C***FIRST EXECUTABLE STATEMENT  QRFAC
 
94
      EPSMCH = R1MACH(4)
 
95
C
 
96
C     COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
 
97
C
 
98
      DO 10 J = 1, N
 
99
         ACNORM(J) = ENORM(M,A(1,J))
 
100
         SIGMA(J) = ACNORM(J)
 
101
         WA(J) = SIGMA(J)
 
102
         IF (PIVOT) IPVT(J) = J
 
103
   10    CONTINUE
 
104
C
 
105
C     REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
 
106
C
 
107
      MINMN = MIN(M,N)
 
108
      DO 110 J = 1, MINMN
 
109
         IF (.NOT.PIVOT) GO TO 40
 
110
C
 
111
C        BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
 
112
C
 
113
         KMAX = J
 
114
         DO 20 K = J, N
 
115
            IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
 
116
   20       CONTINUE
 
117
         IF (KMAX .EQ. J) GO TO 40
 
118
         DO 30 I = 1, M
 
119
            TEMP = A(I,J)
 
120
            A(I,J) = A(I,KMAX)
 
121
            A(I,KMAX) = TEMP
 
122
   30       CONTINUE
 
123
         SIGMA(KMAX) = SIGMA(J)
 
124
         WA(KMAX) = WA(J)
 
125
         K = IPVT(J)
 
126
         IPVT(J) = IPVT(KMAX)
 
127
         IPVT(KMAX) = K
 
128
   40    CONTINUE
 
129
C
 
130
C        COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
 
131
C        J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
 
132
C
 
133
         AJNORM = ENORM(M-J+1,A(J,J))
 
134
         IF (AJNORM .EQ. ZERO) GO TO 100
 
135
         IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
 
136
         DO 50 I = J, M
 
137
            A(I,J) = A(I,J)/AJNORM
 
138
   50       CONTINUE
 
139
         A(J,J) = A(J,J) + ONE
 
140
C
 
141
C        APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
 
142
C        AND UPDATE THE NORMS.
 
143
C
 
144
         JP1 = J + 1
 
145
         IF (N .LT. JP1) GO TO 100
 
146
         DO 90 K = JP1, N
 
147
            SUM = ZERO
 
148
            DO 60 I = J, M
 
149
               SUM = SUM + A(I,J)*A(I,K)
 
150
   60          CONTINUE
 
151
            TEMP = SUM/A(J,J)
 
152
            DO 70 I = J, M
 
153
               A(I,K) = A(I,K) - TEMP*A(I,J)
 
154
   70          CONTINUE
 
155
            IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
 
156
            TEMP = A(J,K)/SIGMA(K)
 
157
            SIGMA(K) = SIGMA(K)*SQRT(MAX(ZERO,ONE-TEMP**2))
 
158
            IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
 
159
            SIGMA(K) = ENORM(M-J,A(JP1,K))
 
160
            WA(K) = SIGMA(K)
 
161
   80       CONTINUE
 
162
   90       CONTINUE
 
163
  100    CONTINUE
 
164
         SIGMA(J) = -AJNORM
 
165
  110    CONTINUE
 
166
      RETURN
 
167
C
 
168
C     LAST CARD OF SUBROUTINE QRFAC.
 
169
C
 
170
      END