~ubuntu-branches/debian/wheezy/arpack/wheezy

« back to all changes in this revision

Viewing changes to ARPACK/BLAS/ctbsv.f

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-12-10 20:32:45 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20111210203245-g0fo30pqvuo92fqh
Tags: 3.0-1
* Switch to arpack-ng since upstream is dead.
* New upstream release
* Daniel Leidert removed from the uploaders (Closes: #651351)

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
      SUBROUTINE CTBSV ( 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
 
      COMPLEX            A( LDA, * ), X( * )
7
 
*     ..
8
 
*
9
 
*  Purpose
10
 
*  =======
11
 
*
12
 
*  CTBSV  solves one of the systems of equations
13
 
*
14
 
*     A*x = b,   or   A'*x = b,   or   conjg( 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'   conjg( 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      - COMPLEX          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      - COMPLEX          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
 
      COMPLEX            ZERO
143
 
      PARAMETER        ( ZERO = ( 0.0E+0, 0.0E+0 ) )
144
 
*     .. Local Scalars ..
145
 
      COMPLEX            TEMP
146
 
      INTEGER            I, INFO, IX, J, JX, KPLUS1, KX, L
147
 
      LOGICAL            NOCONJ, NOUNIT
148
 
*     .. External Functions ..
149
 
      LOGICAL            LSAME
150
 
      EXTERNAL           LSAME
151
 
*     .. External Subroutines ..
152
 
      EXTERNAL           XERBLA
153
 
*     .. Intrinsic Functions ..
154
 
      INTRINSIC          CONJG, 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( 'CTBSV ', INFO )
182
 
         RETURN
183
 
      END IF
184
 
*
185
 
*     Quick return if possible.
186
 
*
187
 
      IF( N.EQ.0 )
188
 
     $   RETURN
189
 
*
190
 
      NOCONJ = LSAME( TRANS, 'T' )
191
 
      NOUNIT = LSAME( DIAG , 'N' )
192
 
*
193
 
*     Set up the start point in X if the increment is not unity. This
194
 
*     will be  ( N - 1 )*INCX  too small for descending loops.
195
 
*
196
 
      IF( INCX.LE.0 )THEN
197
 
         KX = 1 - ( N - 1 )*INCX
198
 
      ELSE IF( INCX.NE.1 )THEN
199
 
         KX = 1
200
 
      END IF
201
 
*
202
 
*     Start the operations. In this version the elements of A are
203
 
*     accessed by sequentially with one pass through A.
204
 
*
205
 
      IF( LSAME( TRANS, 'N' ) )THEN
206
 
*
207
 
*        Form  x := inv( A )*x.
208
 
*
209
 
         IF( LSAME( UPLO, 'U' ) )THEN
210
 
            KPLUS1 = K + 1
211
 
            IF( INCX.EQ.1 )THEN
212
 
               DO 20, J = N, 1, -1
213
 
                  IF( X( J ).NE.ZERO )THEN
214
 
                     L = KPLUS1 - J
215
 
                     IF( NOUNIT )
216
 
     $                  X( J ) = X( J )/A( KPLUS1, J )
217
 
                     TEMP = X( J )
218
 
                     DO 10, I = J - 1, MAX( 1, J - K ), -1
219
 
                        X( I ) = X( I ) - TEMP*A( L + I, J )
220
 
   10                CONTINUE
221
 
                  END IF
222
 
   20          CONTINUE
223
 
            ELSE
224
 
               KX = KX + ( N - 1 )*INCX
225
 
               JX = KX
226
 
               DO 40, J = N, 1, -1
227
 
                  KX = KX - INCX
228
 
                  IF( X( JX ).NE.ZERO )THEN
229
 
                     IX = KX
230
 
                     L  = KPLUS1 - J
231
 
                     IF( NOUNIT )
232
 
     $                  X( JX ) = X( JX )/A( KPLUS1, J )
233
 
                     TEMP = X( JX )
234
 
                     DO 30, I = J - 1, MAX( 1, J - K ), -1
235
 
                        X( IX ) = X( IX ) - TEMP*A( L + I, J )
236
 
                        IX      = IX      - INCX
237
 
   30                CONTINUE
238
 
                  END IF
239
 
                  JX = JX - INCX
240
 
   40          CONTINUE
241
 
            END IF
242
 
         ELSE
243
 
            IF( INCX.EQ.1 )THEN
244
 
               DO 60, J = 1, N
245
 
                  IF( X( J ).NE.ZERO )THEN
246
 
                     L = 1 - J
247
 
                     IF( NOUNIT )
248
 
     $                  X( J ) = X( J )/A( 1, J )
249
 
                     TEMP = X( J )
250
 
                     DO 50, I = J + 1, MIN( N, J + K )
251
 
                        X( I ) = X( I ) - TEMP*A( L + I, J )
252
 
   50                CONTINUE
253
 
                  END IF
254
 
   60          CONTINUE
255
 
            ELSE
256
 
               JX = KX
257
 
               DO 80, J = 1, N
258
 
                  KX = KX + INCX
259
 
                  IF( X( JX ).NE.ZERO )THEN
260
 
                     IX = KX
261
 
                     L  = 1  - J
262
 
                     IF( NOUNIT )
263
 
     $                  X( JX ) = X( JX )/A( 1, J )
264
 
                     TEMP = X( JX )
265
 
                     DO 70, I = J + 1, MIN( N, J + K )
266
 
                        X( IX ) = X( IX ) - TEMP*A( L + I, J )
267
 
                        IX      = IX      + INCX
268
 
   70                CONTINUE
269
 
                  END IF
270
 
                  JX = JX + INCX
271
 
   80          CONTINUE
272
 
            END IF
273
 
         END IF
274
 
      ELSE
275
 
*
276
 
*        Form  x := inv( A' )*x  or  x := inv( conjg( A') )*x.
277
 
*
278
 
         IF( LSAME( UPLO, 'U' ) )THEN
279
 
            KPLUS1 = K + 1
280
 
            IF( INCX.EQ.1 )THEN
281
 
               DO 110, J = 1, N
282
 
                  TEMP = X( J )
283
 
                  L    = KPLUS1 - J
284
 
                  IF( NOCONJ )THEN
285
 
                     DO 90, I = MAX( 1, J - K ), J - 1
286
 
                        TEMP = TEMP - A( L + I, J )*X( I )
287
 
   90                CONTINUE
288
 
                     IF( NOUNIT )
289
 
     $                  TEMP = TEMP/A( KPLUS1, J )
290
 
                  ELSE
291
 
                     DO 100, I = MAX( 1, J - K ), J - 1
292
 
                        TEMP = TEMP - CONJG( A( L + I, J ) )*X( I )
293
 
  100                CONTINUE
294
 
                     IF( NOUNIT )
295
 
     $                  TEMP = TEMP/CONJG( A( KPLUS1, J ) )
296
 
                  END IF
297
 
                  X( J ) = TEMP
298
 
  110          CONTINUE
299
 
            ELSE
300
 
               JX = KX
301
 
               DO 140, J = 1, N
302
 
                  TEMP = X( JX )
303
 
                  IX   = KX
304
 
                  L    = KPLUS1  - J
305
 
                  IF( NOCONJ )THEN
306
 
                     DO 120, I = MAX( 1, J - K ), J - 1
307
 
                        TEMP = TEMP - A( L + I, J )*X( IX )
308
 
                        IX   = IX   + INCX
309
 
  120                CONTINUE
310
 
                     IF( NOUNIT )
311
 
     $                  TEMP = TEMP/A( KPLUS1, J )
312
 
                  ELSE
313
 
                     DO 130, I = MAX( 1, J - K ), J - 1
314
 
                        TEMP = TEMP - CONJG( A( L + I, J ) )*X( IX )
315
 
                        IX   = IX   + INCX
316
 
  130                CONTINUE
317
 
                     IF( NOUNIT )
318
 
     $                  TEMP = TEMP/CONJG( A( KPLUS1, J ) )
319
 
                  END IF
320
 
                  X( JX ) = TEMP
321
 
                  JX      = JX   + INCX
322
 
                  IF( J.GT.K )
323
 
     $               KX = KX + INCX
324
 
  140          CONTINUE
325
 
            END IF
326
 
         ELSE
327
 
            IF( INCX.EQ.1 )THEN
328
 
               DO 170, J = N, 1, -1
329
 
                  TEMP = X( J )
330
 
                  L    = 1      - J
331
 
                  IF( NOCONJ )THEN
332
 
                     DO 150, I = MIN( N, J + K ), J + 1, -1
333
 
                        TEMP = TEMP - A( L + I, J )*X( I )
334
 
  150                CONTINUE
335
 
                     IF( NOUNIT )
336
 
     $                  TEMP = TEMP/A( 1, J )
337
 
                  ELSE
338
 
                     DO 160, I = MIN( N, J + K ), J + 1, -1
339
 
                        TEMP = TEMP - CONJG( A( L + I, J ) )*X( I )
340
 
  160                CONTINUE
341
 
                     IF( NOUNIT )
342
 
     $                  TEMP = TEMP/CONJG( A( 1, J ) )
343
 
                  END IF
344
 
                  X( J ) = TEMP
345
 
  170          CONTINUE
346
 
            ELSE
347
 
               KX = KX + ( N - 1 )*INCX
348
 
               JX = KX
349
 
               DO 200, J = N, 1, -1
350
 
                  TEMP = X( JX )
351
 
                  IX   = KX
352
 
                  L    = 1       - J
353
 
                  IF( NOCONJ )THEN
354
 
                     DO 180, I = MIN( N, J + K ), J + 1, -1
355
 
                        TEMP = TEMP - A( L + I, J )*X( IX )
356
 
                        IX   = IX   - INCX
357
 
  180                CONTINUE
358
 
                     IF( NOUNIT )
359
 
     $                  TEMP = TEMP/A( 1, J )
360
 
                  ELSE
361
 
                     DO 190, I = MIN( N, J + K ), J + 1, -1
362
 
                        TEMP = TEMP - CONJG( A( L + I, J ) )*X( IX )
363
 
                        IX   = IX   - INCX
364
 
  190                CONTINUE
365
 
                     IF( NOUNIT )
366
 
     $                  TEMP = TEMP/CONJG( A( 1, J ) )
367
 
                  END IF
368
 
                  X( JX ) = TEMP
369
 
                  JX      = JX   - INCX
370
 
                  IF( ( N - J ).GE.K )
371
 
     $               KX = KX - INCX
372
 
  200          CONTINUE
373
 
            END IF
374
 
         END IF
375
 
      END IF
376
 
*
377
 
      RETURN
378
 
*
379
 
*     End of CTBSV .
380
 
*
381
 
      END