Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2789 added Finbarr O'Sullivan smoothing spline code

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