← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dchex.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/dchex.f
 
1
*DECK DCHEX
 
2
      SUBROUTINE DCHEX (R, LDR, P, K, L, Z, LDZ, NZ, C, S, JOB)
 
3
C***BEGIN PROLOGUE  DCHEX
 
4
C***PURPOSE  Update the Cholesky factorization  A=TRANS(R)*R  of a
 
5
C            positive definite matrix A of order P under diagonal
 
6
C            permutations of the form  TRANS(E)*A*E, where E is a
 
7
C            permutation matrix.
 
8
C***LIBRARY   SLATEC (LINPACK)
 
9
C***CATEGORY  D7B
 
10
C***TYPE      DOUBLE PRECISION (SCHEX-S, DCHEX-D, CCHEX-C)
 
11
C***KEYWORDS  CHOLESKY DECOMPOSITION, EXCHANGE, LINEAR ALGEBRA, LINPACK,
 
12
C             MATRIX, POSITIVE DEFINITE
 
13
C***AUTHOR  Stewart, G. W., (U. of Maryland)
 
14
C***DESCRIPTION
 
15
C
 
16
C     DCHEX updates the Cholesky factorization
 
17
C
 
18
C                   A = TRANS(R)*R
 
19
C
 
20
C     of a positive definite matrix A of order P under diagonal
 
21
C     permutations of the form
 
22
C
 
23
C                   TRANS(E)*A*E
 
24
C
 
25
C     where E is a permutation matrix.  Specifically, given
 
26
C     an upper triangular matrix R and a permutation matrix
 
27
C     E (which is specified by K, L, and JOB), DCHEX determines
 
28
C     an orthogonal matrix U such that
 
29
C
 
30
C                           U*R*E = RR,
 
31
C
 
32
C     where RR is upper triangular.  At the users option, the
 
33
C     transformation U will be multiplied into the array Z.
 
34
C     If A = TRANS(X)*X, so that R is the triangular part of the
 
35
C     QR factorization of X, then RR is the triangular part of the
 
36
C     QR factorization of X*E, i.e. X with its columns permuted.
 
37
C     For a less terse description of what DCHEX does and how
 
38
C     it may be applied, see the LINPACK guide.
 
39
C
 
40
C     The matrix Q is determined as the product U(L-K)*...*U(1)
 
41
C     of plane rotations of the form
 
42
C
 
43
C                           (    C(I)       S(I) )
 
44
C                           (                    ) ,
 
45
C                           (    -S(I)      C(I) )
 
46
C
 
47
C     where C(I) is double precision.  The rows these rotations operate
 
48
C     on are described below.
 
49
C
 
50
C     There are two types of permutations, which are determined
 
51
C     by the value of JOB.
 
52
C
 
53
C     1. Right circular shift (JOB = 1).
 
54
C
 
55
C         The columns are rearranged in the following order.
 
56
C
 
57
C                1,...,K-1,L,K,K+1,...,L-1,L+1,...,P.
 
58
C
 
59
C         U is the product of L-K rotations U(I), where U(I)
 
60
C         acts in the (L-I,L-I+1)-plane.
 
61
C
 
62
C     2. Left circular shift (JOB = 2).
 
63
C         The columns are rearranged in the following order
 
64
C
 
65
C                1,...,K-1,K+1,K+2,...,L,K,L+1,...,P.
 
66
C
 
67
C         U is the product of L-K rotations U(I), where U(I)
 
68
C         acts in the (K+I-1,K+I)-plane.
 
69
C
 
70
C     On Entry
 
71
C
 
72
C         R      DOUBLE PRECISION(LDR,P), where LDR .GE. P.
 
73
C                R contains the upper triangular factor
 
74
C                that is to be updated.  Elements of R
 
75
C                below the diagonal are not referenced.
 
76
C
 
77
C         LDR    INTEGER.
 
78
C                LDR is the leading dimension of the array R.
 
79
C
 
80
C         P      INTEGER.
 
81
C                P is the order of the matrix R.
 
82
C
 
83
C         K      INTEGER.
 
84
C                K is the first column to be permuted.
 
85
C
 
86
C         L      INTEGER.
 
87
C                L is the last column to be permuted.
 
88
C                L must be strictly greater than K.
 
89
C
 
90
C         Z      DOUBLE PRECISION(LDZ,N)Z), where LDZ .GE. P.
 
91
C                Z is an array of NZ P-vectors into which the
 
92
C                transformation U is multiplied.  Z is
 
93
C                not referenced if NZ = 0.
 
94
C
 
95
C         LDZ    INTEGER.
 
96
C                LDZ is the leading dimension of the array Z.
 
97
C
 
98
C         NZ     INTEGER.
 
99
C                NZ is the number of columns of the matrix Z.
 
100
C
 
101
C         JOB    INTEGER.
 
102
C                JOB determines the type of permutation.
 
103
C                       JOB = 1  right circular shift.
 
104
C                       JOB = 2  left circular shift.
 
105
C
 
106
C     On Return
 
107
C
 
108
C         R      contains the updated factor.
 
109
C
 
110
C         Z      contains the updated matrix Z.
 
111
C
 
112
C         C      DOUBLE PRECISION(P).
 
113
C                C contains the cosines of the transforming rotations.
 
114
C
 
115
C         S      DOUBLE PRECISION(P).
 
116
C                S contains the sines of the transforming rotations.
 
117
C
 
118
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
119
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
120
C***ROUTINES CALLED  DROTG
 
121
C***REVISION HISTORY  (YYMMDD)
 
122
C   780814  DATE WRITTEN
 
123
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
124
C   890831  Modified array declarations.  (WRB)
 
125
C   890831  REVISION DATE from Version 3.2
 
126
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
127
C   900326  Removed duplicate information from DESCRIPTION section.
 
128
C           (WRB)
 
129
C   920501  Reformatted the REFERENCES section.  (WRB)
 
130
C***END PROLOGUE  DCHEX
 
131
      INTEGER LDR,P,K,L,LDZ,NZ,JOB
 
132
      DOUBLE PRECISION R(LDR,*),Z(LDZ,*),S(*)
 
133
      DOUBLE PRECISION C(*)
 
134
C
 
135
      INTEGER I,II,IL,IU,J,JJ,KM1,KP1,LMK,LM1
 
136
      DOUBLE PRECISION T
 
137
C
 
138
C     INITIALIZE
 
139
C
 
140
C***FIRST EXECUTABLE STATEMENT  DCHEX
 
141
      KM1 = K - 1
 
142
      KP1 = K + 1
 
143
      LMK = L - K
 
144
      LM1 = L - 1
 
145
C
 
146
C     PERFORM THE APPROPRIATE TASK.
 
147
C
 
148
      GO TO (10,130), JOB
 
149
C
 
150
C     RIGHT CIRCULAR SHIFT.
 
151
C
 
152
   10 CONTINUE
 
153
C
 
154
C        REORDER THE COLUMNS.
 
155
C
 
156
         DO 20 I = 1, L
 
157
            II = L - I + 1
 
158
            S(I) = R(II,L)
 
159
   20    CONTINUE
 
160
         DO 40 JJ = K, LM1
 
161
            J = LM1 - JJ + K
 
162
            DO 30 I = 1, J
 
163
               R(I,J+1) = R(I,J)
 
164
   30       CONTINUE
 
165
            R(J+1,J+1) = 0.0D0
 
166
   40    CONTINUE
 
167
         IF (K .EQ. 1) GO TO 60
 
168
            DO 50 I = 1, KM1
 
169
               II = L - I + 1
 
170
               R(I,K) = S(II)
 
171
   50       CONTINUE
 
172
   60    CONTINUE
 
173
C
 
174
C        CALCULATE THE ROTATIONS.
 
175
C
 
176
         T = S(1)
 
177
         DO 70 I = 1, LMK
 
178
            CALL DROTG(S(I+1),T,C(I),S(I))
 
179
            T = S(I+1)
 
180
   70    CONTINUE
 
181
         R(K,K) = T
 
182
         DO 90 J = KP1, P
 
183
            IL = MAX(1,L-J+1)
 
184
            DO 80 II = IL, LMK
 
185
               I = L - II
 
186
               T = C(II)*R(I,J) + S(II)*R(I+1,J)
 
187
               R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J)
 
188
               R(I,J) = T
 
189
   80       CONTINUE
 
190
   90    CONTINUE
 
191
C
 
192
C        IF REQUIRED, APPLY THE TRANSFORMATIONS TO Z.
 
193
C
 
194
         IF (NZ .LT. 1) GO TO 120
 
195
         DO 110 J = 1, NZ
 
196
            DO 100 II = 1, LMK
 
197
               I = L - II
 
198
               T = C(II)*Z(I,J) + S(II)*Z(I+1,J)
 
199
               Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J)
 
200
               Z(I,J) = T
 
201
  100       CONTINUE
 
202
  110    CONTINUE
 
203
  120    CONTINUE
 
204
      GO TO 260
 
205
C
 
206
C     LEFT CIRCULAR SHIFT
 
207
C
 
208
  130 CONTINUE
 
209
C
 
210
C        REORDER THE COLUMNS
 
211
C
 
212
         DO 140 I = 1, K
 
213
            II = LMK + I
 
214
            S(II) = R(I,K)
 
215
  140    CONTINUE
 
216
         DO 160 J = K, LM1
 
217
            DO 150 I = 1, J
 
218
               R(I,J) = R(I,J+1)
 
219
  150       CONTINUE
 
220
            JJ = J - KM1
 
221
            S(JJ) = R(J+1,J+1)
 
222
  160    CONTINUE
 
223
         DO 170 I = 1, K
 
224
            II = LMK + I
 
225
            R(I,L) = S(II)
 
226
  170    CONTINUE
 
227
         DO 180 I = KP1, L
 
228
            R(I,L) = 0.0D0
 
229
  180    CONTINUE
 
230
C
 
231
C        REDUCTION LOOP.
 
232
C
 
233
         DO 220 J = K, P
 
234
            IF (J .EQ. K) GO TO 200
 
235
C
 
236
C              APPLY THE ROTATIONS.
 
237
C
 
238
               IU = MIN(J-1,L-1)
 
239
               DO 190 I = K, IU
 
240
                  II = I - K + 1
 
241
                  T = C(II)*R(I,J) + S(II)*R(I+1,J)
 
242
                  R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J)
 
243
                  R(I,J) = T
 
244
  190          CONTINUE
 
245
  200       CONTINUE
 
246
            IF (J .GE. L) GO TO 210
 
247
               JJ = J - K + 1
 
248
               T = S(JJ)
 
249
               CALL DROTG(R(J,J),T,C(JJ),S(JJ))
 
250
  210       CONTINUE
 
251
  220    CONTINUE
 
252
C
 
253
C        APPLY THE ROTATIONS TO Z.
 
254
C
 
255
         IF (NZ .LT. 1) GO TO 250
 
256
         DO 240 J = 1, NZ
 
257
            DO 230 I = K, LM1
 
258
               II = I - KM1
 
259
               T = C(II)*Z(I,J) + S(II)*Z(I+1,J)
 
260
               Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J)
 
261
               Z(I,J) = T
 
262
  230       CONTINUE
 
263
  240    CONTINUE
 
264
  250    CONTINUE
 
265
  260 CONTINUE
 
266
      RETURN
 
267
      END