← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dchdc.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/dchdc.f
 
1
*DECK DCHDC
 
2
      SUBROUTINE DCHDC (A, LDA, P, WORK, JPVT, JOB, INFO)
 
3
C***BEGIN PROLOGUE  DCHDC
 
4
C***PURPOSE  Compute the Cholesky decomposition of a positive definite
 
5
C            matrix.  A pivoting option allows the user to estimate the
 
6
C            condition number of a positive definite matrix or determine
 
7
C            the rank of a positive semidefinite matrix.
 
8
C***LIBRARY   SLATEC (LINPACK)
 
9
C***CATEGORY  D2B1B
 
10
C***TYPE      DOUBLE PRECISION (SCHDC-S, DCHDC-D, CCHDC-C)
 
11
C***KEYWORDS  CHOLESKY DECOMPOSITION, LINEAR ALGEBRA, LINPACK, MATRIX,
 
12
C             POSITIVE DEFINITE
 
13
C***AUTHOR  Dongarra, J., (ANL)
 
14
C           Stewart, G. W., (U. of Maryland)
 
15
C***DESCRIPTION
 
16
C
 
17
C     DCHDC computes the Cholesky decomposition of a positive definite
 
18
C     matrix.  A pivoting option allows the user to estimate the
 
19
C     condition of a positive definite matrix or determine the rank
 
20
C     of a positive semidefinite matrix.
 
21
C
 
22
C     On Entry
 
23
C
 
24
C         A      DOUBLE PRECISION(LDA,P).
 
25
C                A contains the matrix whose decomposition is to
 
26
C                be computed.  Only the upper half of A need be stored.
 
27
C                The lower part of the array A is not referenced.
 
28
C
 
29
C         LDA    INTEGER.
 
30
C                LDA is the leading dimension of the array A.
 
31
C
 
32
C         P      INTEGER.
 
33
C                P is the order of the matrix.
 
34
C
 
35
C         WORK   DOUBLE PRECISION.
 
36
C                WORK is a work array.
 
37
C
 
38
C         JPVT   INTEGER(P).
 
39
C                JPVT contains integers that control the selection
 
40
C                of the pivot elements, if pivoting has been requested.
 
41
C                Each diagonal element A(K,K)
 
42
C                is placed in one of three classes according to the
 
43
C                value of JPVT(K).
 
44
C
 
45
C                   If JPVT(K) .GT. 0, then X(K) is an initial
 
46
C                                      element.
 
47
C
 
48
C                   If JPVT(K) .EQ. 0, then X(K) is a free element.
 
49
C
 
50
C                   If JPVT(K) .LT. 0, then X(K) is a final element.
 
51
C
 
52
C                Before the decomposition is computed, initial elements
 
53
C                are moved by symmetric row and column interchanges to
 
54
C                the beginning of the array A and final
 
55
C                elements to the end.  Both initial and final elements
 
56
C                are frozen in place during the computation and only
 
57
C                free elements are moved.  At the K-th stage of the
 
58
C                reduction, if A(K,K) is occupied by a free element
 
59
C                it is interchanged with the largest free element
 
60
C                A(L,L) with L .GE. K.  JPVT is not referenced if
 
61
C                JOB .EQ. 0.
 
62
C
 
63
C        JOB     INTEGER.
 
64
C                JOB is an integer that initiates column pivoting.
 
65
C                If JOB .EQ. 0, no pivoting is done.
 
66
C                If JOB .NE. 0, pivoting is done.
 
67
C
 
68
C     On Return
 
69
C
 
70
C         A      A contains in its upper half the Cholesky factor
 
71
C                of the matrix A as it has been permuted by pivoting.
 
72
C
 
73
C         JPVT   JPVT(J) contains the index of the diagonal element
 
74
C                of a that was moved into the J-th position,
 
75
C                provided pivoting was requested.
 
76
C
 
77
C         INFO   contains the index of the last positive diagonal
 
78
C                element of the Cholesky factor.
 
79
C
 
80
C     For positive definite matrices INFO = P is the normal return.
 
81
C     For pivoting with positive semidefinite matrices INFO will
 
82
C     in general be less than P.  However, INFO may be greater than
 
83
C     the rank of A, since rounding error can cause an otherwise zero
 
84
C     element to be positive.  Indefinite systems will always cause
 
85
C     INFO to be less than P.
 
86
C
 
87
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
88
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
89
C***ROUTINES CALLED  DAXPY, DSWAP
 
90
C***REVISION HISTORY  (YYMMDD)
 
91
C   790319  DATE WRITTEN
 
92
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
93
C   890831  Modified array declarations.  (WRB)
 
94
C   890831  REVISION DATE from Version 3.2
 
95
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
96
C   900326  Removed duplicate information from DESCRIPTION section.
 
97
C           (WRB)
 
98
C   920501  Reformatted the REFERENCES section.  (WRB)
 
99
C***END PROLOGUE  DCHDC
 
100
      INTEGER LDA,P,JPVT(*),JOB,INFO
 
101
      DOUBLE PRECISION A(LDA,*),WORK(*)
 
102
C
 
103
      INTEGER PU,PL,PLP1,J,JP,JT,K,KB,KM1,KP1,L,MAXL
 
104
      DOUBLE PRECISION TEMP
 
105
      DOUBLE PRECISION MAXDIA
 
106
      LOGICAL SWAPK,NEGK
 
107
C***FIRST EXECUTABLE STATEMENT  DCHDC
 
108
      PL = 1
 
109
      PU = 0
 
110
      INFO = P
 
111
      IF (JOB .EQ. 0) GO TO 160
 
112
C
 
113
C        PIVOTING HAS BEEN REQUESTED. REARRANGE THE
 
114
C        THE ELEMENTS ACCORDING TO JPVT.
 
115
C
 
116
         DO 70 K = 1, P
 
117
            SWAPK = JPVT(K) .GT. 0
 
118
            NEGK = JPVT(K) .LT. 0
 
119
            JPVT(K) = K
 
120
            IF (NEGK) JPVT(K) = -JPVT(K)
 
121
            IF (.NOT.SWAPK) GO TO 60
 
122
               IF (K .EQ. PL) GO TO 50
 
123
                  CALL DSWAP(PL-1,A(1,K),1,A(1,PL),1)
 
124
                  TEMP = A(K,K)
 
125
                  A(K,K) = A(PL,PL)
 
126
                  A(PL,PL) = TEMP
 
127
                  PLP1 = PL + 1
 
128
                  IF (P .LT. PLP1) GO TO 40
 
129
                  DO 30 J = PLP1, P
 
130
                     IF (J .GE. K) GO TO 10
 
131
                        TEMP = A(PL,J)
 
132
                        A(PL,J) = A(J,K)
 
133
                        A(J,K) = TEMP
 
134
                     GO TO 20
 
135
   10                CONTINUE
 
136
                     IF (J .EQ. K) GO TO 20
 
137
                        TEMP = A(K,J)
 
138
                        A(K,J) = A(PL,J)
 
139
                        A(PL,J) = TEMP
 
140
   20                CONTINUE
 
141
   30             CONTINUE
 
142
   40             CONTINUE
 
143
                  JPVT(K) = JPVT(PL)
 
144
                  JPVT(PL) = K
 
145
   50          CONTINUE
 
146
               PL = PL + 1
 
147
   60       CONTINUE
 
148
   70    CONTINUE
 
149
         PU = P
 
150
         IF (P .LT. PL) GO TO 150
 
151
         DO 140 KB = PL, P
 
152
            K = P - KB + PL
 
153
            IF (JPVT(K) .GE. 0) GO TO 130
 
154
               JPVT(K) = -JPVT(K)
 
155
               IF (PU .EQ. K) GO TO 120
 
156
                  CALL DSWAP(K-1,A(1,K),1,A(1,PU),1)
 
157
                  TEMP = A(K,K)
 
158
                  A(K,K) = A(PU,PU)
 
159
                  A(PU,PU) = TEMP
 
160
                  KP1 = K + 1
 
161
                  IF (P .LT. KP1) GO TO 110
 
162
                  DO 100 J = KP1, P
 
163
                     IF (J .GE. PU) GO TO 80
 
164
                        TEMP = A(K,J)
 
165
                        A(K,J) = A(J,PU)
 
166
                        A(J,PU) = TEMP
 
167
                     GO TO 90
 
168
   80                CONTINUE
 
169
                     IF (J .EQ. PU) GO TO 90
 
170
                        TEMP = A(K,J)
 
171
                        A(K,J) = A(PU,J)
 
172
                        A(PU,J) = TEMP
 
173
   90                CONTINUE
 
174
  100             CONTINUE
 
175
  110             CONTINUE
 
176
                  JT = JPVT(K)
 
177
                  JPVT(K) = JPVT(PU)
 
178
                  JPVT(PU) = JT
 
179
  120          CONTINUE
 
180
               PU = PU - 1
 
181
  130       CONTINUE
 
182
  140    CONTINUE
 
183
  150    CONTINUE
 
184
  160 CONTINUE
 
185
      DO 270 K = 1, P
 
186
C
 
187
C        REDUCTION LOOP.
 
188
C
 
189
         MAXDIA = A(K,K)
 
190
         KP1 = K + 1
 
191
         MAXL = K
 
192
C
 
193
C        DETERMINE THE PIVOT ELEMENT.
 
194
C
 
195
         IF (K .LT. PL .OR. K .GE. PU) GO TO 190
 
196
            DO 180 L = KP1, PU
 
197
               IF (A(L,L) .LE. MAXDIA) GO TO 170
 
198
                  MAXDIA = A(L,L)
 
199
                  MAXL = L
 
200
  170          CONTINUE
 
201
  180       CONTINUE
 
202
  190    CONTINUE
 
203
C
 
204
C        QUIT IF THE PIVOT ELEMENT IS NOT POSITIVE.
 
205
C
 
206
         IF (MAXDIA .GT. 0.0D0) GO TO 200
 
207
            INFO = K - 1
 
208
            GO TO 280
 
209
  200    CONTINUE
 
210
         IF (K .EQ. MAXL) GO TO 210
 
211
C
 
212
C           START THE PIVOTING AND UPDATE JPVT.
 
213
C
 
214
            KM1 = K - 1
 
215
            CALL DSWAP(KM1,A(1,K),1,A(1,MAXL),1)
 
216
            A(MAXL,MAXL) = A(K,K)
 
217
            A(K,K) = MAXDIA
 
218
            JP = JPVT(MAXL)
 
219
            JPVT(MAXL) = JPVT(K)
 
220
            JPVT(K) = JP
 
221
  210    CONTINUE
 
222
C
 
223
C        REDUCTION STEP. PIVOTING IS CONTAINED ACROSS THE ROWS.
 
224
C
 
225
         WORK(K) = SQRT(A(K,K))
 
226
         A(K,K) = WORK(K)
 
227
         IF (P .LT. KP1) GO TO 260
 
228
         DO 250 J = KP1, P
 
229
            IF (K .EQ. MAXL) GO TO 240
 
230
               IF (J .GE. MAXL) GO TO 220
 
231
                  TEMP = A(K,J)
 
232
                  A(K,J) = A(J,MAXL)
 
233
                  A(J,MAXL) = TEMP
 
234
               GO TO 230
 
235
  220          CONTINUE
 
236
               IF (J .EQ. MAXL) GO TO 230
 
237
                  TEMP = A(K,J)
 
238
                  A(K,J) = A(MAXL,J)
 
239
                  A(MAXL,J) = TEMP
 
240
  230          CONTINUE
 
241
  240       CONTINUE
 
242
            A(K,J) = A(K,J)/WORK(K)
 
243
            WORK(J) = A(K,J)
 
244
            TEMP = -A(K,J)
 
245
            CALL DAXPY(J-K,TEMP,WORK(KP1),1,A(KP1,J),1)
 
246
  250    CONTINUE
 
247
  260    CONTINUE
 
248
  270 CONTINUE
 
249
  280 CONTINUE
 
250
      RETURN
 
251
      END