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