Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/ctrmv.f
5176 views
1
SUBROUTINE CTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
2
* .. Scalar Arguments ..
3
INTEGER INCX, LDA, N
4
CHARACTER*1 DIAG, TRANS, UPLO
5
* .. Array Arguments ..
6
COMPLEX A( LDA, * ), X( * )
7
* ..
8
*
9
* Purpose
10
* =======
11
*
12
* CTRMV performs one of the matrix-vector operations
13
*
14
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
15
*
16
* where x is an n element vector and A is an n by n unit, or non-unit,
17
* upper or lower triangular matrix.
18
*
19
* Parameters
20
* ==========
21
*
22
* UPLO - CHARACTER*1.
23
* On entry, UPLO specifies whether the matrix is an upper or
24
* lower triangular matrix as follows:
25
*
26
* UPLO = 'U' or 'u' A is an upper triangular matrix.
27
*
28
* UPLO = 'L' or 'l' A is a lower triangular matrix.
29
*
30
* Unchanged on exit.
31
*
32
* TRANS - CHARACTER*1.
33
* On entry, TRANS specifies the operation to be performed as
34
* follows:
35
*
36
* TRANS = 'N' or 'n' x := A*x.
37
*
38
* TRANS = 'T' or 't' x := A'*x.
39
*
40
* TRANS = 'C' or 'c' x := conjg( A' )*x.
41
*
42
* Unchanged on exit.
43
*
44
* DIAG - CHARACTER*1.
45
* On entry, DIAG specifies whether or not A is unit
46
* triangular as follows:
47
*
48
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
49
*
50
* DIAG = 'N' or 'n' A is not assumed to be unit
51
* triangular.
52
*
53
* Unchanged on exit.
54
*
55
* N - INTEGER.
56
* On entry, N specifies the order of the matrix A.
57
* N must be at least zero.
58
* Unchanged on exit.
59
*
60
* A - COMPLEX array of DIMENSION ( LDA, n ).
61
* Before entry with UPLO = 'U' or 'u', the leading n by n
62
* upper triangular part of the array A must contain the upper
63
* triangular matrix and the strictly lower triangular part of
64
* A is not referenced.
65
* Before entry with UPLO = 'L' or 'l', the leading n by n
66
* lower triangular part of the array A must contain the lower
67
* triangular matrix and the strictly upper triangular part of
68
* A is not referenced.
69
* Note that when DIAG = 'U' or 'u', the diagonal elements of
70
* A are not referenced either, but are assumed to be unity.
71
* Unchanged on exit.
72
*
73
* LDA - INTEGER.
74
* On entry, LDA specifies the first dimension of A as declared
75
* in the calling (sub) program. LDA must be at least
76
* max( 1, n ).
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 vector x. On exit, X is overwritten with the
83
* tranformed 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, 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, MAX
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( LDA.LT.MAX( 1, N ) )THEN
133
INFO = 6
134
ELSE IF( INCX.EQ.0 )THEN
135
INFO = 8
136
END IF
137
IF( INFO.NE.0 )THEN
138
CALL XERBLA( 'CTRMV ', INFO )
139
RETURN
140
END IF
141
*
142
* Quick return if possible.
143
*
144
IF( N.EQ.0 )
145
$ RETURN
146
*
147
NOCONJ = LSAME( TRANS, 'T' )
148
NOUNIT = LSAME( DIAG , 'N' )
149
*
150
* Set up the start point in X if the increment is not unity. This
151
* will be ( N - 1 )*INCX too small for descending loops.
152
*
153
IF( INCX.LE.0 )THEN
154
KX = 1 - ( N - 1 )*INCX
155
ELSE IF( INCX.NE.1 )THEN
156
KX = 1
157
END IF
158
*
159
* Start the operations. In this version the elements of A are
160
* accessed sequentially with one pass through A.
161
*
162
IF( LSAME( TRANS, 'N' ) )THEN
163
*
164
* Form x := A*x.
165
*
166
IF( LSAME( UPLO, 'U' ) )THEN
167
IF( INCX.EQ.1 )THEN
168
DO 20, J = 1, N
169
IF( X( J ).NE.ZERO )THEN
170
TEMP = X( J )
171
DO 10, I = 1, J - 1
172
X( I ) = X( I ) + TEMP*A( I, J )
173
10 CONTINUE
174
IF( NOUNIT )
175
$ X( J ) = X( J )*A( J, J )
176
END IF
177
20 CONTINUE
178
ELSE
179
JX = KX
180
DO 40, J = 1, N
181
IF( X( JX ).NE.ZERO )THEN
182
TEMP = X( JX )
183
IX = KX
184
DO 30, I = 1, J - 1
185
X( IX ) = X( IX ) + TEMP*A( I, J )
186
IX = IX + INCX
187
30 CONTINUE
188
IF( NOUNIT )
189
$ X( JX ) = X( JX )*A( J, J )
190
END IF
191
JX = JX + INCX
192
40 CONTINUE
193
END IF
194
ELSE
195
IF( INCX.EQ.1 )THEN
196
DO 60, J = N, 1, -1
197
IF( X( J ).NE.ZERO )THEN
198
TEMP = X( J )
199
DO 50, I = N, J + 1, -1
200
X( I ) = X( I ) + TEMP*A( I, J )
201
50 CONTINUE
202
IF( NOUNIT )
203
$ X( J ) = X( J )*A( J, J )
204
END IF
205
60 CONTINUE
206
ELSE
207
KX = KX + ( N - 1 )*INCX
208
JX = KX
209
DO 80, J = N, 1, -1
210
IF( X( JX ).NE.ZERO )THEN
211
TEMP = X( JX )
212
IX = KX
213
DO 70, I = N, J + 1, -1
214
X( IX ) = X( IX ) + TEMP*A( I, J )
215
IX = IX - INCX
216
70 CONTINUE
217
IF( NOUNIT )
218
$ X( JX ) = X( JX )*A( J, J )
219
END IF
220
JX = JX - INCX
221
80 CONTINUE
222
END IF
223
END IF
224
ELSE
225
*
226
* Form x := A'*x or x := conjg( A' )*x.
227
*
228
IF( LSAME( UPLO, 'U' ) )THEN
229
IF( INCX.EQ.1 )THEN
230
DO 110, J = N, 1, -1
231
TEMP = X( J )
232
IF( NOCONJ )THEN
233
IF( NOUNIT )
234
$ TEMP = TEMP*A( J, J )
235
DO 90, I = J - 1, 1, -1
236
TEMP = TEMP + A( I, J )*X( I )
237
90 CONTINUE
238
ELSE
239
IF( NOUNIT )
240
$ TEMP = TEMP*CONJG( A( J, J ) )
241
DO 100, I = J - 1, 1, -1
242
TEMP = TEMP + CONJG( A( I, J ) )*X( I )
243
100 CONTINUE
244
END IF
245
X( J ) = TEMP
246
110 CONTINUE
247
ELSE
248
JX = KX + ( N - 1 )*INCX
249
DO 140, J = N, 1, -1
250
TEMP = X( JX )
251
IX = JX
252
IF( NOCONJ )THEN
253
IF( NOUNIT )
254
$ TEMP = TEMP*A( J, J )
255
DO 120, I = J - 1, 1, -1
256
IX = IX - INCX
257
TEMP = TEMP + A( I, J )*X( IX )
258
120 CONTINUE
259
ELSE
260
IF( NOUNIT )
261
$ TEMP = TEMP*CONJG( A( J, J ) )
262
DO 130, I = J - 1, 1, -1
263
IX = IX - INCX
264
TEMP = TEMP + CONJG( A( I, J ) )*X( IX )
265
130 CONTINUE
266
END IF
267
X( JX ) = TEMP
268
JX = JX - INCX
269
140 CONTINUE
270
END IF
271
ELSE
272
IF( INCX.EQ.1 )THEN
273
DO 170, J = 1, N
274
TEMP = X( J )
275
IF( NOCONJ )THEN
276
IF( NOUNIT )
277
$ TEMP = TEMP*A( J, J )
278
DO 150, I = J + 1, N
279
TEMP = TEMP + A( I, J )*X( I )
280
150 CONTINUE
281
ELSE
282
IF( NOUNIT )
283
$ TEMP = TEMP*CONJG( A( J, J ) )
284
DO 160, I = J + 1, N
285
TEMP = TEMP + CONJG( A( I, J ) )*X( I )
286
160 CONTINUE
287
END IF
288
X( J ) = TEMP
289
170 CONTINUE
290
ELSE
291
JX = KX
292
DO 200, J = 1, N
293
TEMP = X( JX )
294
IX = JX
295
IF( NOCONJ )THEN
296
IF( NOUNIT )
297
$ TEMP = TEMP*A( J, J )
298
DO 180, I = J + 1, N
299
IX = IX + INCX
300
TEMP = TEMP + A( I, J )*X( IX )
301
180 CONTINUE
302
ELSE
303
IF( NOUNIT )
304
$ TEMP = TEMP*CONJG( A( J, J ) )
305
DO 190, I = J + 1, N
306
IX = IX + INCX
307
TEMP = TEMP + CONJG( A( I, J ) )*X( IX )
308
190 CONTINUE
309
END IF
310
X( JX ) = TEMP
311
JX = JX + INCX
312
200 CONTINUE
313
END IF
314
END IF
315
END IF
316
*
317
RETURN
318
*
319
* End of CTRMV .
320
*
321
END
322
323