← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/blas/dgbmv.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

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