source: branches/iLoveclim/SOURCES/BLAS/stbsv.f @ 254

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

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