← Back to branch summary

~ubuntu-branches/ubuntu/trusty/atlas/trusty

« back to all changes in this revision

Viewing changes to src/blas/f77reference/dtbmv.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-07-27 14:26:05 UTC
  • mfrom: (18.1.8 sid)
  • Revision ID: package-import@ubuntu.com-20130727142605-5rh3p972h1whdo99
Tags: 3.10.1-2
* Allow the generic package to build on machines with CPU throttling
  enabled. Otherwise the package FTBFS on some buildds (e.g. biber).
  Implementation is done by reactivating the "-Si cputhrchk 0" flag
  (cpu-throtthling-check.diff), and using it in debian/rules.
* Add architectural defaults for armel and mips.
* armhf.diff: do not enforce 32-registers FPU for Fortran

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/blas/f77reference/dtbmv.f
 
1
      SUBROUTINE DTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
 
2
*     .. Scalar Arguments ..
 
3
      INTEGER INCX,K,LDA,N
 
4
      CHARACTER DIAG,TRANS,UPLO
 
5
*     ..
 
6
*     .. Array Arguments ..
 
7
      DOUBLE PRECISION A(LDA,*),X(*)
 
8
*     ..
 
9
*
 
10
*  Purpose
 
11
*  =======
 
12
*
 
13
*  DTBMV  performs one of the matrix-vector operations
 
14
*
 
15
*     x := A*x,   or   x := A'*x,
 
16
*
 
17
*  where x is an n element vector and  A is an n by n unit, or non-unit,
 
18
*  upper or lower triangular band matrix, with ( k + 1 ) diagonals.
 
19
*
 
20
*  Arguments
 
21
*  ==========
 
22
*
 
23
*  UPLO   - CHARACTER*1.
 
24
*           On entry, UPLO specifies whether the matrix is an upper or
 
25
*           lower triangular matrix as follows:
 
26
*
 
27
*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
 
28
*
 
29
*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
 
30
*
 
31
*           Unchanged on exit.
 
32
*
 
33
*  TRANS  - CHARACTER*1.
 
34
*           On entry, TRANS specifies the operation to be performed as
 
35
*           follows:
 
36
*
 
37
*              TRANS = 'N' or 'n'   x := A*x.
 
38
*
 
39
*              TRANS = 'T' or 't'   x := A'*x.
 
40
*
 
41
*              TRANS = 'C' or 'c'   x := A'*x.
 
42
*
 
43
*           Unchanged on exit.
 
44
*
 
45
*  DIAG   - CHARACTER*1.
 
46
*           On entry, DIAG specifies whether or not A is unit
 
47
*           triangular as follows:
 
48
*
 
49
*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 
50
*
 
51
*              DIAG = 'N' or 'n'   A is not assumed to be unit
 
52
*                                  triangular.
 
53
*
 
54
*           Unchanged on exit.
 
55
*
 
56
*  N      - INTEGER.
 
57
*           On entry, N specifies the order of the matrix A.
 
58
*           N must be at least zero.
 
59
*           Unchanged on exit.
 
60
*
 
61
*  K      - INTEGER.
 
62
*           On entry with UPLO = 'U' or 'u', K specifies the number of
 
63
*           super-diagonals of the matrix A.
 
64
*           On entry with UPLO = 'L' or 'l', K specifies the number of
 
65
*           sub-diagonals of the matrix A.
 
66
*           K must satisfy  0 .le. K.
 
67
*           Unchanged on exit.
 
68
*
 
69
*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
 
70
*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
 
71
*           by n part of the array A must contain the upper triangular
 
72
*           band part of the matrix of coefficients, supplied column by
 
73
*           column, with the leading diagonal of the matrix in row
 
74
*           ( k + 1 ) of the array, the first super-diagonal starting at
 
75
*           position 2 in row k, and so on. The top left k by k triangle
 
76
*           of the array A is not referenced.
 
77
*           The following program segment will transfer an upper
 
78
*           triangular band matrix from conventional full matrix storage
 
79
*           to band storage:
 
80
*
 
81
*                 DO 20, J = 1, N
 
82
*                    M = K + 1 - J
 
83
*                    DO 10, I = MAX( 1, J - K ), J
 
84
*                       A( M + I, J ) = matrix( I, J )
 
85
*              10    CONTINUE
 
86
*              20 CONTINUE
 
87
*
 
88
*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
 
89
*           by n part of the array A must contain the lower triangular
 
90
*           band part of the matrix of coefficients, supplied column by
 
91
*           column, with the leading diagonal of the matrix in row 1 of
 
92
*           the array, the first sub-diagonal starting at position 1 in
 
93
*           row 2, and so on. The bottom right k by k triangle of the
 
94
*           array A is not referenced.
 
95
*           The following program segment will transfer a lower
 
96
*           triangular band matrix from conventional full matrix storage
 
97
*           to band storage:
 
98
*
 
99
*                 DO 20, J = 1, N
 
100
*                    M = 1 - J
 
101
*                    DO 10, I = J, MIN( N, J + K )
 
102
*                       A( M + I, J ) = matrix( I, J )
 
103
*              10    CONTINUE
 
104
*              20 CONTINUE
 
105
*
 
106
*           Note that when DIAG = 'U' or 'u' the elements of the array A
 
107
*           corresponding to the diagonal elements of the matrix are not
 
108
*           referenced, but are assumed to be unity.
 
109
*           Unchanged on exit.
 
110
*
 
111
*  LDA    - INTEGER.
 
112
*           On entry, LDA specifies the first dimension of A as declared
 
113
*           in the calling (sub) program. LDA must be at least
 
114
*           ( k + 1 ).
 
115
*           Unchanged on exit.
 
116
*
 
117
*  X      - DOUBLE PRECISION array of dimension at least
 
118
*           ( 1 + ( n - 1 )*abs( INCX ) ).
 
119
*           Before entry, the incremented array X must contain the n
 
120
*           element vector x. On exit, X is overwritten with the
 
121
*           tranformed vector x.
 
122
*
 
123
*  INCX   - INTEGER.
 
124
*           On entry, INCX specifies the increment for the elements of
 
125
*           X. INCX must not be zero.
 
126
*           Unchanged on exit.
 
127
*
 
128
*
 
129
*  Level 2 Blas routine.
 
130
*
 
131
*  -- Written on 22-October-1986.
 
132
*     Jack Dongarra, Argonne National Lab.
 
133
*     Jeremy Du Croz, Nag Central Office.
 
134
*     Sven Hammarling, Nag Central Office.
 
135
*     Richard Hanson, Sandia National Labs.
 
136
*
 
137
*
 
138
*     .. Parameters ..
 
139
      DOUBLE PRECISION ZERO
 
140
      PARAMETER (ZERO=0.0D+0)
 
141
*     ..
 
142
*     .. Local Scalars ..
 
143
      DOUBLE PRECISION TEMP
 
144
      INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
 
145
      LOGICAL NOUNIT
 
146
*     ..
 
147
*     .. External Functions ..
 
148
      LOGICAL LSAME
 
149
      EXTERNAL LSAME
 
150
*     ..
 
151
*     .. External Subroutines ..
 
152
      EXTERNAL XERBLA
 
153
*     ..
 
154
*     .. Intrinsic Functions ..
 
155
      INTRINSIC MAX,MIN
 
156
*     ..
 
157
*
 
158
*     Test the input parameters.
 
159
*
 
160
      INFO = 0
 
161
      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
 
162
          INFO = 1
 
163
      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
 
164
     +         .NOT.LSAME(TRANS,'C')) THEN
 
165
          INFO = 2
 
166
      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
 
167
          INFO = 3
 
168
      ELSE IF (N.LT.0) THEN
 
169
          INFO = 4
 
170
      ELSE IF (K.LT.0) THEN
 
171
          INFO = 5
 
172
      ELSE IF (LDA.LT. (K+1)) THEN
 
173
          INFO = 7
 
174
      ELSE IF (INCX.EQ.0) THEN
 
175
          INFO = 9
 
176
      END IF
 
177
      IF (INFO.NE.0) THEN
 
178
          CALL XERBLA('DTBMV ',INFO)
 
179
          RETURN
 
180
      END IF
 
181
*
 
182
*     Quick return if possible.
 
183
*
 
184
      IF (N.EQ.0) RETURN
 
185
*
 
186
      NOUNIT = LSAME(DIAG,'N')
 
187
*
 
188
*     Set up the start point in X if the increment is not unity. This
 
189
*     will be  ( N - 1 )*INCX   too small for descending loops.
 
190
*
 
191
      IF (INCX.LE.0) THEN
 
192
          KX = 1 - (N-1)*INCX
 
193
      ELSE IF (INCX.NE.1) THEN
 
194
          KX = 1
 
195
      END IF
 
196
*
 
197
*     Start the operations. In this version the elements of A are
 
198
*     accessed sequentially with one pass through A.
 
199
*
 
200
      IF (LSAME(TRANS,'N')) THEN
 
201
*
 
202
*         Form  x := A*x.
 
203
*
 
204
          IF (LSAME(UPLO,'U')) THEN
 
205
              KPLUS1 = K + 1
 
206
              IF (INCX.EQ.1) THEN
 
207
                  DO 20 J = 1,N
 
208
                      IF (X(J).NE.ZERO) THEN
 
209
                          TEMP = X(J)
 
210
                          L = KPLUS1 - J
 
211
                          DO 10 I = MAX(1,J-K),J - 1
 
212
                              X(I) = X(I) + TEMP*A(L+I,J)
 
213
   10                     CONTINUE
 
214
                          IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
 
215
                      END IF
 
216
   20             CONTINUE
 
217
              ELSE
 
218
                  JX = KX
 
219
                  DO 40 J = 1,N
 
220
                      IF (X(JX).NE.ZERO) THEN
 
221
                          TEMP = X(JX)
 
222
                          IX = KX
 
223
                          L = KPLUS1 - J
 
224
                          DO 30 I = MAX(1,J-K),J - 1
 
225
                              X(IX) = X(IX) + TEMP*A(L+I,J)
 
226
                              IX = IX + INCX
 
227
   30                     CONTINUE
 
228
                          IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
 
229
                      END IF
 
230
                      JX = JX + INCX
 
231
                      IF (J.GT.K) KX = KX + INCX
 
232
   40             CONTINUE
 
233
              END IF
 
234
          ELSE
 
235
              IF (INCX.EQ.1) THEN
 
236
                  DO 60 J = N,1,-1
 
237
                      IF (X(J).NE.ZERO) THEN
 
238
                          TEMP = X(J)
 
239
                          L = 1 - J
 
240
                          DO 50 I = MIN(N,J+K),J + 1,-1
 
241
                              X(I) = X(I) + TEMP*A(L+I,J)
 
242
   50                     CONTINUE
 
243
                          IF (NOUNIT) X(J) = X(J)*A(1,J)
 
244
                      END IF
 
245
   60             CONTINUE
 
246
              ELSE
 
247
                  KX = KX + (N-1)*INCX
 
248
                  JX = KX
 
249
                  DO 80 J = N,1,-1
 
250
                      IF (X(JX).NE.ZERO) THEN
 
251
                          TEMP = X(JX)
 
252
                          IX = KX
 
253
                          L = 1 - J
 
254
                          DO 70 I = MIN(N,J+K),J + 1,-1
 
255
                              X(IX) = X(IX) + TEMP*A(L+I,J)
 
256
                              IX = IX - INCX
 
257
   70                     CONTINUE
 
258
                          IF (NOUNIT) X(JX) = X(JX)*A(1,J)
 
259
                      END IF
 
260
                      JX = JX - INCX
 
261
                      IF ((N-J).GE.K) KX = KX - INCX
 
262
   80             CONTINUE
 
263
              END IF
 
264
          END IF
 
265
      ELSE
 
266
*
 
267
*        Form  x := A'*x.
 
268
*
 
269
          IF (LSAME(UPLO,'U')) THEN
 
270
              KPLUS1 = K + 1
 
271
              IF (INCX.EQ.1) THEN
 
272
                  DO 100 J = N,1,-1
 
273
                      TEMP = X(J)
 
274
                      L = KPLUS1 - J
 
275
                      IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
 
276
                      DO 90 I = J - 1,MAX(1,J-K),-1
 
277
                          TEMP = TEMP + A(L+I,J)*X(I)
 
278
   90                 CONTINUE
 
279
                      X(J) = TEMP
 
280
  100             CONTINUE
 
281
              ELSE
 
282
                  KX = KX + (N-1)*INCX
 
283
                  JX = KX
 
284
                  DO 120 J = N,1,-1
 
285
                      TEMP = X(JX)
 
286
                      KX = KX - INCX
 
287
                      IX = KX
 
288
                      L = KPLUS1 - J
 
289
                      IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
 
290
                      DO 110 I = J - 1,MAX(1,J-K),-1
 
291
                          TEMP = TEMP + A(L+I,J)*X(IX)
 
292
                          IX = IX - INCX
 
293
  110                 CONTINUE
 
294
                      X(JX) = TEMP
 
295
                      JX = JX - INCX
 
296
  120             CONTINUE
 
297
              END IF
 
298
          ELSE
 
299
              IF (INCX.EQ.1) THEN
 
300
                  DO 140 J = 1,N
 
301
                      TEMP = X(J)
 
302
                      L = 1 - J
 
303
                      IF (NOUNIT) TEMP = TEMP*A(1,J)
 
304
                      DO 130 I = J + 1,MIN(N,J+K)
 
305
                          TEMP = TEMP + A(L+I,J)*X(I)
 
306
  130                 CONTINUE
 
307
                      X(J) = TEMP
 
308
  140             CONTINUE
 
309
              ELSE
 
310
                  JX = KX
 
311
                  DO 160 J = 1,N
 
312
                      TEMP = X(JX)
 
313
                      KX = KX + INCX
 
314
                      IX = KX
 
315
                      L = 1 - J
 
316
                      IF (NOUNIT) TEMP = TEMP*A(1,J)
 
317
                      DO 150 I = J + 1,MIN(N,J+K)
 
318
                          TEMP = TEMP + A(L+I,J)*X(IX)
 
319
                          IX = IX + INCX
 
320
  150                 CONTINUE
 
321
                      X(JX) = TEMP
 
322
                      JX = JX + INCX
 
323
  160             CONTINUE
 
324
              END IF
 
325
          END IF
 
326
      END IF
 
327
*
 
328
      RETURN
 
329
*
 
330
*     End of DTBMV .
 
331
*
 
332
      END