← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/chifa.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/chifa.f
 
1
*DECK CHIFA
 
2
      SUBROUTINE CHIFA (A, LDA, N, KPVT, INFO)
 
3
C***BEGIN PROLOGUE  CHIFA
 
4
C***PURPOSE  Factor a complex Hermitian matrix by elimination
 
5
C            (symmetric pivoting).
 
6
C***LIBRARY   SLATEC (LINPACK)
 
7
C***CATEGORY  D2D1A
 
8
C***TYPE      COMPLEX (SSIFA-S, DSIFA-D, CHIFA-C, CSIFA-C)
 
9
C***KEYWORDS  HERMITIAN, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
 
10
C***AUTHOR  Bunch, J., (UCSD)
 
11
C***DESCRIPTION
 
12
C
 
13
C     CHIFA factors a complex Hermitian matrix by elimination
 
14
C     with symmetric pivoting.
 
15
C
 
16
C     To solve  A*X = B , follow CHIFA by CHISL.
 
17
C     To compute  INVERSE(A)*C , follow CHIFA by CHISL.
 
18
C     To compute  DETERMINANT(A) , follow CHIFA by CHIDI.
 
19
C     To compute  INERTIA(A) , follow CHIFA by CHIDI.
 
20
C     To compute  INVERSE(A) , follow CHIFA by CHIDI.
 
21
C
 
22
C     On Entry
 
23
C
 
24
C        A       COMPLEX(LDA,N)
 
25
C                the Hermitian matrix to be factored.
 
26
C                Only the diagonal and upper triangle are used.
 
27
C
 
28
C        LDA     INTEGER
 
29
C                the leading dimension of the array  A .
 
30
C
 
31
C        N       INTEGER
 
32
C                the order of the matrix  A .
 
33
C
 
34
C     On Return
 
35
C
 
36
C        A       a block diagonal matrix and the multipliers which
 
37
C                were used to obtain it.
 
38
C                The factorization can be written  A = U*D*CTRANS(U)
 
39
C                where  U  is a product of permutation and unit
 
40
C                upper triangular matrices , CTRANS(U) is the
 
41
C                conjugate transpose of  U , and  D  is block diagonal
 
42
C                with 1 by 1 and 2 by 2 blocks.
 
43
C
 
44
C        KVPT    INTEGER(N)
 
45
C                an integer vector of pivot indices.
 
46
C
 
47
C        INFO    INTEGER
 
48
C                = 0  normal value.
 
49
C                = K  if the K-th pivot block is singular.  This is
 
50
C                     not an error condition for this subroutine,
 
51
C                     but it does indicate that CHISL or CHIDI may
 
52
C                     divide by zero if called.
 
53
C
 
54
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
55
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
56
C***ROUTINES CALLED  CAXPY, CSWAP, ICAMAX
 
57
C***REVISION HISTORY  (YYMMDD)
 
58
C   780814  DATE WRITTEN
 
59
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
60
C   890831  Modified array declarations.  (WRB)
 
61
C   891107  Modified routine equivalence list.  (WRB)
 
62
C   891107  REVISION DATE from Version 3.2
 
63
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
64
C   900326  Removed duplicate information from DESCRIPTION section.
 
65
C           (WRB)
 
66
C   920501  Reformatted the REFERENCES section.  (WRB)
 
67
C***END PROLOGUE  CHIFA
 
68
      INTEGER LDA,N,KPVT(*),INFO
 
69
      COMPLEX A(LDA,*)
 
70
C
 
71
      COMPLEX AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T
 
72
      REAL ABSAKK,ALPHA,COLMAX,ROWMAX
 
73
      INTEGER IMAX,IMAXP1,J,JJ,JMAX,K,KM1,KM2,KSTEP,ICAMAX
 
74
      LOGICAL SWAP
 
75
      COMPLEX ZDUM
 
76
      REAL CABS1
 
77
      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
 
78
C***FIRST EXECUTABLE STATEMENT  CHIFA
 
79
C
 
80
C     INITIALIZE
 
81
C
 
82
C     ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
 
83
C
 
84
      ALPHA = (1.0E0 + SQRT(17.0E0))/8.0E0
 
85
C
 
86
      INFO = 0
 
87
C
 
88
C     MAIN LOOP ON K, WHICH GOES FROM N TO 1.
 
89
C
 
90
      K = N
 
91
   10 CONTINUE
 
92
C
 
93
C        LEAVE THE LOOP IF K=0 OR K=1.
 
94
C
 
95
         IF (K .EQ. 0) GO TO 200
 
96
         IF (K .GT. 1) GO TO 20
 
97
            KPVT(1) = 1
 
98
            IF (CABS1(A(1,1)) .EQ. 0.0E0) INFO = 1
 
99
            GO TO 200
 
100
   20    CONTINUE
 
101
C
 
102
C        THIS SECTION OF CODE DETERMINES THE KIND OF
 
103
C        ELIMINATION TO BE PERFORMED.  WHEN IT IS COMPLETED,
 
104
C        KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND
 
105
C        SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS
 
106
C        REQUIRED.
 
107
C
 
108
         KM1 = K - 1
 
109
         ABSAKK = CABS1(A(K,K))
 
110
C
 
111
C        DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
 
112
C        COLUMN K.
 
113
C
 
114
         IMAX = ICAMAX(K-1,A(1,K),1)
 
115
         COLMAX = CABS1(A(IMAX,K))
 
116
         IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30
 
117
            KSTEP = 1
 
118
            SWAP = .FALSE.
 
119
         GO TO 90
 
120
   30    CONTINUE
 
121
C
 
122
C           DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
 
123
C           ROW IMAX.
 
124
C
 
125
            ROWMAX = 0.0E0
 
126
            IMAXP1 = IMAX + 1
 
127
            DO 40 J = IMAXP1, K
 
128
               ROWMAX = MAX(ROWMAX,CABS1(A(IMAX,J)))
 
129
   40       CONTINUE
 
130
            IF (IMAX .EQ. 1) GO TO 50
 
131
               JMAX = ICAMAX(IMAX-1,A(1,IMAX),1)
 
132
               ROWMAX = MAX(ROWMAX,CABS1(A(JMAX,IMAX)))
 
133
   50       CONTINUE
 
134
            IF (CABS1(A(IMAX,IMAX)) .LT. ALPHA*ROWMAX) GO TO 60
 
135
               KSTEP = 1
 
136
               SWAP = .TRUE.
 
137
            GO TO 80
 
138
   60       CONTINUE
 
139
            IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70
 
140
               KSTEP = 1
 
141
               SWAP = .FALSE.
 
142
            GO TO 80
 
143
   70       CONTINUE
 
144
               KSTEP = 2
 
145
               SWAP = IMAX .NE. KM1
 
146
   80       CONTINUE
 
147
   90    CONTINUE
 
148
         IF (MAX(ABSAKK,COLMAX) .NE. 0.0E0) GO TO 100
 
149
C
 
150
C           COLUMN K IS ZERO.  SET INFO AND ITERATE THE LOOP.
 
151
C
 
152
            KPVT(K) = K
 
153
            INFO = K
 
154
         GO TO 190
 
155
  100    CONTINUE
 
156
         IF (KSTEP .EQ. 2) GO TO 140
 
157
C
 
158
C           1 X 1 PIVOT BLOCK.
 
159
C
 
160
            IF (.NOT.SWAP) GO TO 120
 
161
C
 
162
C              PERFORM AN INTERCHANGE.
 
163
C
 
164
               CALL CSWAP(IMAX,A(1,IMAX),1,A(1,K),1)
 
165
               DO 110 JJ = IMAX, K
 
166
                  J = K + IMAX - JJ
 
167
                  T = CONJG(A(J,K))
 
168
                  A(J,K) = CONJG(A(IMAX,J))
 
169
                  A(IMAX,J) = T
 
170
  110          CONTINUE
 
171
  120       CONTINUE
 
172
C
 
173
C           PERFORM THE ELIMINATION.
 
174
C
 
175
            DO 130 JJ = 1, KM1
 
176
               J = K - JJ
 
177
               MULK = -A(J,K)/A(K,K)
 
178
               T = CONJG(MULK)
 
179
               CALL CAXPY(J,T,A(1,K),1,A(1,J),1)
 
180
               A(J,J) = CMPLX(REAL(A(J,J)),0.0E0)
 
181
               A(J,K) = MULK
 
182
  130       CONTINUE
 
183
C
 
184
C           SET THE PIVOT ARRAY.
 
185
C
 
186
            KPVT(K) = K
 
187
            IF (SWAP) KPVT(K) = IMAX
 
188
         GO TO 190
 
189
  140    CONTINUE
 
190
C
 
191
C           2 X 2 PIVOT BLOCK.
 
192
C
 
193
            IF (.NOT.SWAP) GO TO 160
 
194
C
 
195
C              PERFORM AN INTERCHANGE.
 
196
C
 
197
               CALL CSWAP(IMAX,A(1,IMAX),1,A(1,K-1),1)
 
198
               DO 150 JJ = IMAX, KM1
 
199
                  J = KM1 + IMAX - JJ
 
200
                  T = CONJG(A(J,K-1))
 
201
                  A(J,K-1) = CONJG(A(IMAX,J))
 
202
                  A(IMAX,J) = T
 
203
  150          CONTINUE
 
204
               T = A(K-1,K)
 
205
               A(K-1,K) = A(IMAX,K)
 
206
               A(IMAX,K) = T
 
207
  160       CONTINUE
 
208
C
 
209
C           PERFORM THE ELIMINATION.
 
210
C
 
211
            KM2 = K - 2
 
212
            IF (KM2 .EQ. 0) GO TO 180
 
213
               AK = A(K,K)/A(K-1,K)
 
214
               AKM1 = A(K-1,K-1)/CONJG(A(K-1,K))
 
215
               DENOM = 1.0E0 - AK*AKM1
 
216
               DO 170 JJ = 1, KM2
 
217
                  J = KM1 - JJ
 
218
                  BK = A(J,K)/A(K-1,K)
 
219
                  BKM1 = A(J,K-1)/CONJG(A(K-1,K))
 
220
                  MULK = (AKM1*BK - BKM1)/DENOM
 
221
                  MULKM1 = (AK*BKM1 - BK)/DENOM
 
222
                  T = CONJG(MULK)
 
223
                  CALL CAXPY(J,T,A(1,K),1,A(1,J),1)
 
224
                  T = CONJG(MULKM1)
 
225
                  CALL CAXPY(J,T,A(1,K-1),1,A(1,J),1)
 
226
                  A(J,K) = MULK
 
227
                  A(J,K-1) = MULKM1
 
228
                  A(J,J) = CMPLX(REAL(A(J,J)),0.0E0)
 
229
  170          CONTINUE
 
230
  180       CONTINUE
 
231
C
 
232
C           SET THE PIVOT ARRAY.
 
233
C
 
234
            KPVT(K) = 1 - K
 
235
            IF (SWAP) KPVT(K) = -IMAX
 
236
            KPVT(K-1) = KPVT(K)
 
237
  190    CONTINUE
 
238
         K = K - KSTEP
 
239
      GO TO 10
 
240
  200 CONTINUE
 
241
      RETURN
 
242
      END