← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/sspco.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/sspco.f
 
1
*DECK SSPCO
 
2
      SUBROUTINE SSPCO (AP, N, KPVT, RCOND, Z)
 
3
C***BEGIN PROLOGUE  SSPCO
 
4
C***PURPOSE  Factor a real symmetric matrix stored in packed form
 
5
C            by elimination with symmetric pivoting and estimate the
 
6
C            condition number of the matrix.
 
7
C***LIBRARY   SLATEC (LINPACK)
 
8
C***CATEGORY  D2B1A
 
9
C***TYPE      SINGLE PRECISION (SSPCO-S, DSPCO-D, CHPCO-C, CSPCO-C)
 
10
C***KEYWORDS  CONDITION NUMBER, LINEAR ALGEBRA, LINPACK,
 
11
C             MATRIX FACTORIZATION, PACKED, SYMMETRIC
 
12
C***AUTHOR  Moler, C. B., (U. of New Mexico)
 
13
C***DESCRIPTION
 
14
C
 
15
C     SSPCO factors a real symmetric matrix stored in packed
 
16
C     form by elimination with symmetric pivoting and estimates
 
17
C     the condition of the matrix.
 
18
C
 
19
C     If  RCOND  is not needed, SSPFA is slightly faster.
 
20
C     To solve  A*X = B , follow SSPCO by SSPSL.
 
21
C     To compute  INVERSE(A)*C , follow SSPCO by SSPSL.
 
22
C     To compute  INVERSE(A) , follow SSPCO by SSPDI.
 
23
C     To compute  DETERMINANT(A) , follow SSPCO by SSPDI.
 
24
C     To compute  INERTIA(A), follow SSPCO by SSPDI.
 
25
C
 
26
C     On Entry
 
27
C
 
28
C        AP      REAL (N*(N+1)/2)
 
29
C                the packed form of a symmetric matrix  A .  The
 
30
C                columns of the upper triangle are stored sequentially
 
31
C                in a one-dimensional array of length  N*(N+1)/2 .
 
32
C                See comments below for details.
 
33
C
 
34
C        N       INTEGER
 
35
C                the order of the matrix  A .
 
36
C
 
37
C     Output
 
38
C
 
39
C        AP      a block diagonal matrix and the multipliers which
 
40
C                were used to obtain it stored in packed form.
 
41
C                The factorization can be written  A = U*D*TRANS(U)
 
42
C                where  U  is a product of permutation and unit
 
43
C                upper triangular matrices , TRANS(U) is the
 
44
C                transpose of  U , and  D  is block diagonal
 
45
C                with 1 by 1 and 2 by 2 blocks.
 
46
C
 
47
C        KPVT    INTEGER(N)
 
48
C                an integer vector of pivot indices.
 
49
C
 
50
C        RCOND   REAL
 
51
C                an estimate of the reciprocal condition of  A .
 
52
C                For the system  A*X = B , relative perturbations
 
53
C                in  A  and  B  of size  EPSILON  may cause
 
54
C                relative perturbations in  X  of size  EPSILON/RCOND .
 
55
C                If  RCOND  is so small that the logical expression
 
56
C                           1.0 + RCOND .EQ. 1.0
 
57
C                is true, then  A  may be singular to working
 
58
C                precision.  In particular,  RCOND  is zero  if
 
59
C                exact singularity is detected or the estimate
 
60
C                underflows.
 
61
C
 
62
C        Z       REAL(N)
 
63
C                a work vector whose contents are usually unimportant.
 
64
C                If  A  is close to a singular matrix, then  Z  is
 
65
C                an approximate null vector in the sense that
 
66
C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
 
67
C
 
68
C     Packed Storage
 
69
C
 
70
C          The following program segment will pack the upper
 
71
C          triangle of a symmetric matrix.
 
72
C
 
73
C                K = 0
 
74
C                DO 20 J = 1, N
 
75
C                   DO 10 I = 1, J
 
76
C                      K = K + 1
 
77
C                      AP(K) = A(I,J)
 
78
C             10    CONTINUE
 
79
C             20 CONTINUE
 
80
C
 
81
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
82
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
83
C***ROUTINES CALLED  SASUM, SAXPY, SDOT, SSCAL, SSPFA
 
84
C***REVISION HISTORY  (YYMMDD)
 
85
C   780814  DATE WRITTEN
 
86
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
87
C   890831  Modified array declarations.  (WRB)
 
88
C   891107  Modified routine equivalence list.  (WRB)
 
89
C   891107  REVISION DATE from Version 3.2
 
90
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
91
C   900326  Removed duplicate information from DESCRIPTION section.
 
92
C           (WRB)
 
93
C   920501  Reformatted the REFERENCES section.  (WRB)
 
94
C***END PROLOGUE  SSPCO
 
95
      INTEGER N,KPVT(*)
 
96
      REAL AP(*),Z(*)
 
97
      REAL RCOND
 
98
C
 
99
      REAL AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T
 
100
      REAL ANORM,S,SASUM,YNORM
 
101
      INTEGER I,IJ,IK,IKM1,IKP1,INFO,J,JM1,J1
 
102
      INTEGER K,KK,KM1K,KM1KM1,KP,KPS,KS
 
103
C
 
104
C     FIND NORM OF A USING ONLY UPPER HALF
 
105
C
 
106
C***FIRST EXECUTABLE STATEMENT  SSPCO
 
107
      J1 = 1
 
108
      DO 30 J = 1, N
 
109
         Z(J) = SASUM(J,AP(J1),1)
 
110
         IJ = J1
 
111
         J1 = J1 + J
 
112
         JM1 = J - 1
 
113
         IF (JM1 .LT. 1) GO TO 20
 
114
         DO 10 I = 1, JM1
 
115
            Z(I) = Z(I) + ABS(AP(IJ))
 
116
            IJ = IJ + 1
 
117
   10    CONTINUE
 
118
   20    CONTINUE
 
119
   30 CONTINUE
 
120
      ANORM = 0.0E0
 
121
      DO 40 J = 1, N
 
122
         ANORM = MAX(ANORM,Z(J))
 
123
   40 CONTINUE
 
124
C
 
125
C     FACTOR
 
126
C
 
127
      CALL SSPFA(AP,N,KPVT,INFO)
 
128
C
 
129
C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
 
130
C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
 
131
C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
 
132
C     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
 
133
C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
 
134
C
 
135
C     SOLVE U*D*W = E
 
136
C
 
137
      EK = 1.0E0
 
138
      DO 50 J = 1, N
 
139
         Z(J) = 0.0E0
 
140
   50 CONTINUE
 
141
      K = N
 
142
      IK = (N*(N - 1))/2
 
143
   60 IF (K .EQ. 0) GO TO 120
 
144
         KK = IK + K
 
145
         IKM1 = IK - (K - 1)
 
146
         KS = 1
 
147
         IF (KPVT(K) .LT. 0) KS = 2
 
148
         KP = ABS(KPVT(K))
 
149
         KPS = K + 1 - KS
 
150
         IF (KP .EQ. KPS) GO TO 70
 
151
            T = Z(KPS)
 
152
            Z(KPS) = Z(KP)
 
153
            Z(KP) = T
 
154
   70    CONTINUE
 
155
         IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,Z(K))
 
156
         Z(K) = Z(K) + EK
 
157
         CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
 
158
         IF (KS .EQ. 1) GO TO 80
 
159
            IF (Z(K-1) .NE. 0.0E0) EK = SIGN(EK,Z(K-1))
 
160
            Z(K-1) = Z(K-1) + EK
 
161
            CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
 
162
   80    CONTINUE
 
163
         IF (KS .EQ. 2) GO TO 100
 
164
            IF (ABS(Z(K)) .LE. ABS(AP(KK))) GO TO 90
 
165
               S = ABS(AP(KK))/ABS(Z(K))
 
166
               CALL SSCAL(N,S,Z,1)
 
167
               EK = S*EK
 
168
   90       CONTINUE
 
169
            IF (AP(KK) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
 
170
            IF (AP(KK) .EQ. 0.0E0) Z(K) = 1.0E0
 
171
         GO TO 110
 
172
  100    CONTINUE
 
173
            KM1K = IK + K - 1
 
174
            KM1KM1 = IKM1 + K - 1
 
175
            AK = AP(KK)/AP(KM1K)
 
176
            AKM1 = AP(KM1KM1)/AP(KM1K)
 
177
            BK = Z(K)/AP(KM1K)
 
178
            BKM1 = Z(K-1)/AP(KM1K)
 
179
            DENOM = AK*AKM1 - 1.0E0
 
180
            Z(K) = (AKM1*BK - BKM1)/DENOM
 
181
            Z(K-1) = (AK*BKM1 - BK)/DENOM
 
182
  110    CONTINUE
 
183
         K = K - KS
 
184
         IK = IK - K
 
185
         IF (KS .EQ. 2) IK = IK - (K + 1)
 
186
      GO TO 60
 
187
  120 CONTINUE
 
188
      S = 1.0E0/SASUM(N,Z,1)
 
189
      CALL SSCAL(N,S,Z,1)
 
190
C
 
191
C     SOLVE TRANS(U)*Y = W
 
192
C
 
193
      K = 1
 
194
      IK = 0
 
195
  130 IF (K .GT. N) GO TO 160
 
196
         KS = 1
 
197
         IF (KPVT(K) .LT. 0) KS = 2
 
198
         IF (K .EQ. 1) GO TO 150
 
199
            Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1)
 
200
            IKP1 = IK + K
 
201
            IF (KS .EQ. 2)
 
202
     1         Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1)
 
203
            KP = ABS(KPVT(K))
 
204
            IF (KP .EQ. K) GO TO 140
 
205
               T = Z(K)
 
206
               Z(K) = Z(KP)
 
207
               Z(KP) = T
 
208
  140       CONTINUE
 
209
  150    CONTINUE
 
210
         IK = IK + K
 
211
         IF (KS .EQ. 2) IK = IK + (K + 1)
 
212
         K = K + KS
 
213
      GO TO 130
 
214
  160 CONTINUE
 
215
      S = 1.0E0/SASUM(N,Z,1)
 
216
      CALL SSCAL(N,S,Z,1)
 
217
C
 
218
      YNORM = 1.0E0
 
219
C
 
220
C     SOLVE U*D*V = Y
 
221
C
 
222
      K = N
 
223
      IK = N*(N - 1)/2
 
224
  170 IF (K .EQ. 0) GO TO 230
 
225
         KK = IK + K
 
226
         IKM1 = IK - (K - 1)
 
227
         KS = 1
 
228
         IF (KPVT(K) .LT. 0) KS = 2
 
229
         IF (K .EQ. KS) GO TO 190
 
230
            KP = ABS(KPVT(K))
 
231
            KPS = K + 1 - KS
 
232
            IF (KP .EQ. KPS) GO TO 180
 
233
               T = Z(KPS)
 
234
               Z(KPS) = Z(KP)
 
235
               Z(KP) = T
 
236
  180       CONTINUE
 
237
            CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
 
238
            IF (KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
 
239
  190    CONTINUE
 
240
         IF (KS .EQ. 2) GO TO 210
 
241
            IF (ABS(Z(K)) .LE. ABS(AP(KK))) GO TO 200
 
242
               S = ABS(AP(KK))/ABS(Z(K))
 
243
               CALL SSCAL(N,S,Z,1)
 
244
               YNORM = S*YNORM
 
245
  200       CONTINUE
 
246
            IF (AP(KK) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
 
247
            IF (AP(KK) .EQ. 0.0E0) Z(K) = 1.0E0
 
248
         GO TO 220
 
249
  210    CONTINUE
 
250
            KM1K = IK + K - 1
 
251
            KM1KM1 = IKM1 + K - 1
 
252
            AK = AP(KK)/AP(KM1K)
 
253
            AKM1 = AP(KM1KM1)/AP(KM1K)
 
254
            BK = Z(K)/AP(KM1K)
 
255
            BKM1 = Z(K-1)/AP(KM1K)
 
256
            DENOM = AK*AKM1 - 1.0E0
 
257
            Z(K) = (AKM1*BK - BKM1)/DENOM
 
258
            Z(K-1) = (AK*BKM1 - BK)/DENOM
 
259
  220    CONTINUE
 
260
         K = K - KS
 
261
         IK = IK - K
 
262
         IF (KS .EQ. 2) IK = IK - (K + 1)
 
263
      GO TO 170
 
264
  230 CONTINUE
 
265
      S = 1.0E0/SASUM(N,Z,1)
 
266
      CALL SSCAL(N,S,Z,1)
 
267
      YNORM = S*YNORM
 
268
C
 
269
C     SOLVE TRANS(U)*Z = V
 
270
C
 
271
      K = 1
 
272
      IK = 0
 
273
  240 IF (K .GT. N) GO TO 270
 
274
         KS = 1
 
275
         IF (KPVT(K) .LT. 0) KS = 2
 
276
         IF (K .EQ. 1) GO TO 260
 
277
            Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1)
 
278
            IKP1 = IK + K
 
279
            IF (KS .EQ. 2)
 
280
     1         Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1)
 
281
            KP = ABS(KPVT(K))
 
282
            IF (KP .EQ. K) GO TO 250
 
283
               T = Z(K)
 
284
               Z(K) = Z(KP)
 
285
               Z(KP) = T
 
286
  250       CONTINUE
 
287
  260    CONTINUE
 
288
         IK = IK + K
 
289
         IF (KS .EQ. 2) IK = IK + (K + 1)
 
290
         K = K + KS
 
291
      GO TO 240
 
292
  270 CONTINUE
 
293
C     MAKE ZNORM = 1.0
 
294
      S = 1.0E0/SASUM(N,Z,1)
 
295
      CALL SSCAL(N,S,Z,1)
 
296
      YNORM = S*YNORM
 
297
C
 
298
      IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
 
299
      IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
 
300
      RETURN
 
301
      END