source: branches/GRISLIv3/SOURCES/BLAS/sgbmv.f @ 474

Last change on this file since 474 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk

File size: 9.1 KB
Line 
1      SUBROUTINE SGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
2     $                   BETA, Y, INCY )
3*     .. Scalar Arguments ..
4      REAL               ALPHA, BETA
5      INTEGER            INCX, INCY, KL, KU, LDA, M, N
6      CHARACTER*1        TRANS
7*     .. Array Arguments ..
8      REAL               A( LDA, * ), X( * ), Y( * )
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  SGBMV  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  - REAL            .
57*           On entry, ALPHA specifies the scalar alpha.
58*           Unchanged on exit.
59*
60*  A      - REAL             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      - REAL             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   - REAL            .
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      - REAL             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      REAL               ONE         , ZERO
129      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
130*     .. Local Scalars ..
131      REAL               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( 'SGBMV ', 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 SGBMV .
299*
300      END
Note: See TracBrowser for help on using the repository browser.