Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/sgbmv.f @ 15678

Last change on this file since 15678 was 15457, checked in by gkronber, 7 years ago

#2789 added Finbarr O'Sullivan smoothing spline code

File size: 8.8 KB
Line 
1      SUBROUTINE SGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2*     .. Scalar Arguments ..
3      REAL ALPHA,BETA
4      INTEGER INCX,INCY,KL,KU,LDA,M,N
5      CHARACTER TRANS
6*     ..
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*  Arguments
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*     ..
131*     .. Local Scalars ..
132      REAL TEMP
133      INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
134*     ..
135*     .. External Functions ..
136      LOGICAL LSAME
137      EXTERNAL LSAME
138*     ..
139*     .. External Subroutines ..
140      EXTERNAL XERBLA
141*     ..
142*     .. Intrinsic Functions ..
143      INTRINSIC MAX,MIN
144*     ..
145*
146*     Test the input parameters.
147*
148      INFO = 0
149      IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
150     +    .NOT.LSAME(TRANS,'C')) THEN
151          INFO = 1
152      ELSE IF (M.LT.0) THEN
153          INFO = 2
154      ELSE IF (N.LT.0) THEN
155          INFO = 3
156      ELSE IF (KL.LT.0) THEN
157          INFO = 4
158      ELSE IF (KU.LT.0) THEN
159          INFO = 5
160      ELSE IF (LDA.LT. (KL+KU+1)) THEN
161          INFO = 8
162      ELSE IF (INCX.EQ.0) THEN
163          INFO = 10
164      ELSE IF (INCY.EQ.0) THEN
165          INFO = 13
166      END IF
167      IF (INFO.NE.0) THEN
168          CALL XERBLA('SGBMV ',INFO)
169          RETURN
170      END IF
171*
172*     Quick return if possible.
173*
174      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
175     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) 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) RETURN
230      KUP1 = KU + 1
231      IF (LSAME(TRANS,'N')) THEN
232*
233*        Form  y := alpha*A*x + y.
234*
235          JX = KX
236          IF (INCY.EQ.1) THEN
237              DO 60 J = 1,N
238                  IF (X(JX).NE.ZERO) THEN
239                      TEMP = ALPHA*X(JX)
240                      K = KUP1 - J
241                      DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
242                          Y(I) = Y(I) + TEMP*A(K+I,J)
243   50                 CONTINUE
244                  END IF
245                  JX = JX + INCX
246   60         CONTINUE
247          ELSE
248              DO 80 J = 1,N
249                  IF (X(JX).NE.ZERO) THEN
250                      TEMP = ALPHA*X(JX)
251                      IY = KY
252                      K = KUP1 - J
253                      DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
254                          Y(IY) = Y(IY) + TEMP*A(K+I,J)
255                          IY = IY + INCY
256   70                 CONTINUE
257                  END IF
258                  JX = JX + INCX
259                  IF (J.GT.KU) KY = KY + INCY
260   80         CONTINUE
261          END IF
262      ELSE
263*
264*        Form  y := alpha*A'*x + y.
265*
266          JY = KY
267          IF (INCX.EQ.1) THEN
268              DO 100 J = 1,N
269                  TEMP = ZERO
270                  K = KUP1 - J
271                  DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
272                      TEMP = TEMP + A(K+I,J)*X(I)
273   90             CONTINUE
274                  Y(JY) = Y(JY) + ALPHA*TEMP
275                  JY = JY + INCY
276  100         CONTINUE
277          ELSE
278              DO 120 J = 1,N
279                  TEMP = ZERO
280                  IX = KX
281                  K = KUP1 - J
282                  DO 110 I = MAX(1,J-KU),MIN(M,J+KL)
283                      TEMP = TEMP + A(K+I,J)*X(IX)
284                      IX = IX + INCX
285  110             CONTINUE
286                  Y(JY) = Y(JY) + ALPHA*TEMP
287                  JY = JY + INCY
288                  IF (J.GT.KU) KX = KX + INCX
289  120         CONTINUE
290          END IF
291      END IF
292*
293      RETURN
294*
295*     End of SGBMV .
296*
297      END
Note: See TracBrowser for help on using the repository browser.