← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/dgbfa.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/dgbfa.f
 
1
*DECK DGBFA
 
2
      SUBROUTINE DGBFA (ABD, LDA, N, ML, MU, IPVT, INFO)
 
3
C***BEGIN PROLOGUE  DGBFA
 
4
C***PURPOSE  Factor a band matrix using Gaussian elimination.
 
5
C***LIBRARY   SLATEC (LINPACK)
 
6
C***CATEGORY  D2A2
 
7
C***TYPE      DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
 
8
C***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
 
9
C***AUTHOR  Moler, C. B., (U. of New Mexico)
 
10
C***DESCRIPTION
 
11
C
 
12
C     DGBFA factors a double precision band matrix by elimination.
 
13
C
 
14
C     DGBFA is usually called by DGBCO, but it can be called
 
15
C     directly with a saving in time if  RCOND  is not needed.
 
16
C
 
17
C     On Entry
 
18
C
 
19
C        ABD     DOUBLE PRECISION(LDA, N)
 
20
C                contains the matrix in band storage.  The columns
 
21
C                of the matrix are stored in the columns of  ABD  and
 
22
C                the diagonals of the matrix are stored in rows
 
23
C                ML+1 through 2*ML+MU+1 of  ABD .
 
24
C                See the comments below for details.
 
25
C
 
26
C        LDA     INTEGER
 
27
C                the leading dimension of the array  ABD .
 
28
C                LDA must be .GE. 2*ML + MU + 1 .
 
29
C
 
30
C        N       INTEGER
 
31
C                the order of the original matrix.
 
32
C
 
33
C        ML      INTEGER
 
34
C                number of diagonals below the main diagonal.
 
35
C                0 .LE. ML .LT.  N .
 
36
C
 
37
C        MU      INTEGER
 
38
C                number of diagonals above the main diagonal.
 
39
C                0 .LE. MU .LT.  N .
 
40
C                More efficient if  ML .LE. MU .
 
41
C     On Return
 
42
C
 
43
C        ABD     an upper triangular matrix in band storage and
 
44
C                the multipliers which were used to obtain it.
 
45
C                The factorization can be written  A = L*U  where
 
46
C                L  is a product of permutation and unit lower
 
47
C                triangular matrices and  U  is upper triangular.
 
48
C
 
49
C        IPVT    INTEGER(N)
 
50
C                an integer vector of pivot indices.
 
51
C
 
52
C        INFO    INTEGER
 
53
C                = 0  normal value.
 
54
C                = K  if  U(K,K) .EQ. 0.0 .  This is not an error
 
55
C                     condition for this subroutine, but it does
 
56
C                     indicate that DGBSL will divide by zero if
 
57
C                     called.  Use  RCOND  in DGBCO for a reliable
 
58
C                     indication of singularity.
 
59
C
 
60
C     Band Storage
 
61
C
 
62
C           If  A  is a band matrix, the following program segment
 
63
C           will set up the input.
 
64
C
 
65
C                   ML = (band width below the diagonal)
 
66
C                   MU = (band width above the diagonal)
 
67
C                   M = ML + MU + 1
 
68
C                   DO 20 J = 1, N
 
69
C                      I1 = MAX(1, J-MU)
 
70
C                      I2 = MIN(N, J+ML)
 
71
C                      DO 10 I = I1, I2
 
72
C                         K = I - J + M
 
73
C                         ABD(K,J) = A(I,J)
 
74
C                10    CONTINUE
 
75
C                20 CONTINUE
 
76
C
 
77
C           This uses rows  ML+1  through  2*ML+MU+1  of  ABD .
 
78
C           In addition, the first  ML  rows in  ABD  are used for
 
79
C           elements generated during the triangularization.
 
80
C           The total number of rows needed in  ABD  is  2*ML+MU+1 .
 
81
C           The  ML+MU by ML+MU  upper left triangle and the
 
82
C           ML by ML  lower right triangle are not referenced.
 
83
C
 
84
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
85
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
86
C***ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
 
87
C***REVISION HISTORY  (YYMMDD)
 
88
C   780814  DATE WRITTEN
 
89
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
90
C   890831  Modified array declarations.  (WRB)
 
91
C   890831  REVISION DATE from Version 3.2
 
92
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
93
C   900326  Removed duplicate information from DESCRIPTION section.
 
94
C           (WRB)
 
95
C   920501  Reformatted the REFERENCES section.  (WRB)
 
96
C***END PROLOGUE  DGBFA
 
97
      INTEGER LDA,N,ML,MU,IPVT(*),INFO
 
98
      DOUBLE PRECISION ABD(LDA,*)
 
99
C
 
100
      DOUBLE PRECISION T
 
101
      INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
 
102
C
 
103
C***FIRST EXECUTABLE STATEMENT  DGBFA
 
104
      M = ML + MU + 1
 
105
      INFO = 0
 
106
C
 
107
C     ZERO INITIAL FILL-IN COLUMNS
 
108
C
 
109
      J0 = MU + 2
 
110
      J1 = MIN(N,M) - 1
 
111
      IF (J1 .LT. J0) GO TO 30
 
112
      DO 20 JZ = J0, J1
 
113
         I0 = M + 1 - JZ
 
114
         DO 10 I = I0, ML
 
115
            ABD(I,JZ) = 0.0D0
 
116
   10    CONTINUE
 
117
   20 CONTINUE
 
118
   30 CONTINUE
 
119
      JZ = J1
 
120
      JU = 0
 
121
C
 
122
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
 
123
C
 
124
      NM1 = N - 1
 
125
      IF (NM1 .LT. 1) GO TO 130
 
126
      DO 120 K = 1, NM1
 
127
         KP1 = K + 1
 
128
C
 
129
C        ZERO NEXT FILL-IN COLUMN
 
130
C
 
131
         JZ = JZ + 1
 
132
         IF (JZ .GT. N) GO TO 50
 
133
         IF (ML .LT. 1) GO TO 50
 
134
            DO 40 I = 1, ML
 
135
               ABD(I,JZ) = 0.0D0
 
136
   40       CONTINUE
 
137
   50    CONTINUE
 
138
C
 
139
C        FIND L = PIVOT INDEX
 
140
C
 
141
         LM = MIN(ML,N-K)
 
142
         L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
 
143
         IPVT(K) = L + K - M
 
144
C
 
145
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
 
146
C
 
147
         IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
 
148
C
 
149
C           INTERCHANGE IF NECESSARY
 
150
C
 
151
            IF (L .EQ. M) GO TO 60
 
152
               T = ABD(L,K)
 
153
               ABD(L,K) = ABD(M,K)
 
154
               ABD(M,K) = T
 
155
   60       CONTINUE
 
156
C
 
157
C           COMPUTE MULTIPLIERS
 
158
C
 
159
            T = -1.0D0/ABD(M,K)
 
160
            CALL DSCAL(LM,T,ABD(M+1,K),1)
 
161
C
 
162
C           ROW ELIMINATION WITH COLUMN INDEXING
 
163
C
 
164
            JU = MIN(MAX(JU,MU+IPVT(K)),N)
 
165
            MM = M
 
166
            IF (JU .LT. KP1) GO TO 90
 
167
            DO 80 J = KP1, JU
 
168
               L = L - 1
 
169
               MM = MM - 1
 
170
               T = ABD(L,J)
 
171
               IF (L .EQ. MM) GO TO 70
 
172
                  ABD(L,J) = ABD(MM,J)
 
173
                  ABD(MM,J) = T
 
174
   70          CONTINUE
 
175
               CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
 
176
   80       CONTINUE
 
177
   90       CONTINUE
 
178
         GO TO 110
 
179
  100    CONTINUE
 
180
            INFO = K
 
181
  110    CONTINUE
 
182
  120 CONTINUE
 
183
  130 CONTINUE
 
184
      IPVT(N) = N
 
185
      IF (ABD(M,N) .EQ. 0.0D0) INFO = N
 
186
      RETURN
 
187
      END