120
121
122
123
124
125
126 CHARACTER UPLO
127 INTEGER INFO, LDA, LDB, N, NRHS
128
129
130 INTEGER IPIV( * )
131 COMPLEX*16 A( LDA, * ), B( LDB, * )
132
133
134
135
136
137 COMPLEX*16 ONE
138 parameter( one = ( 1.0d+0, 0.0d+0 ) )
139
140
141 LOGICAL UPPER
142 INTEGER J, K, KP
143 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
144
145
146 LOGICAL LSAME
148
149
151
152
153 INTRINSIC max
154
155
156
157 info = 0
158 upper =
lsame( uplo,
'U' )
159 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
160 info = -1
161 ELSE IF( n.LT.0 ) THEN
162 info = -2
163 ELSE IF( nrhs.LT.0 ) THEN
164 info = -3
165 ELSE IF( lda.LT.max( 1, n ) ) THEN
166 info = -5
167 ELSE IF( ldb.LT.max( 1, n ) ) THEN
168 info = -8
169 END IF
170 IF( info.NE.0 ) THEN
171 CALL xerbla(
'ZSYTRS', -info )
172 RETURN
173 END IF
174
175
176
177 IF( n.EQ.0 .OR. nrhs.EQ.0 )
178 $ RETURN
179
180 IF( upper ) THEN
181
182
183
184
185
186
187
188
189 k = n
190 10 CONTINUE
191
192
193
194 IF( k.LT.1 )
195 $ GO TO 30
196
197 IF( ipiv( k ).GT.0 ) THEN
198
199
200
201
202
203 kp = ipiv( k )
204 IF( kp.NE.k )
205 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
206
207
208
209
210 CALL zgeru( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
211 $ b( 1, 1 ), ldb )
212
213
214
215 CALL zscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
216 k = k - 1
217 ELSE
218
219
220
221
222
223 kp = -ipiv( k )
224 IF( kp.NE.k-1 )
225 $
CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
226
227
228
229
230 CALL zgeru( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
231 $ b( 1, 1 ), ldb )
232 CALL zgeru( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 1 ),
233 $ ldb, b( 1, 1 ), ldb )
234
235
236
237 akm1k = a( k-1, k )
238 akm1 = a( k-1, k-1 ) / akm1k
239 ak = a( k, k ) / akm1k
240 denom = akm1*ak - one
241 DO 20 j = 1, nrhs
242 bkm1 = b( k-1, j ) / akm1k
243 bk = b( k, j ) / akm1k
244 b( k-1, j ) = ( ak*bkm1-bk ) / denom
245 b( k, j ) = ( akm1*bk-bkm1 ) / denom
246 20 CONTINUE
247 k = k - 2
248 END IF
249
250 GO TO 10
251 30 CONTINUE
252
253
254
255
256
257
258 k = 1
259 40 CONTINUE
260
261
262
263 IF( k.GT.n )
264 $ GO TO 50
265
266 IF( ipiv( k ).GT.0 ) THEN
267
268
269
270
271
272
273 CALL zgemv(
'Transpose', k-1, nrhs, -one, b, ldb, a( 1, k ),
274 $ 1, one, b( k, 1 ), ldb )
275
276
277
278 kp = ipiv( k )
279 IF( kp.NE.k )
280 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
281 k = k + 1
282 ELSE
283
284
285
286
287
288
289 CALL zgemv(
'Transpose', k-1, nrhs, -one, b, ldb, a( 1, k ),
290 $ 1, one, b( k, 1 ), ldb )
291 CALL zgemv(
'Transpose', k-1, nrhs, -one, b, ldb,
292 $ a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
293
294
295
296 kp = -ipiv( k )
297 IF( kp.NE.k )
298 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
299 k = k + 2
300 END IF
301
302 GO TO 40
303 50 CONTINUE
304
305 ELSE
306
307
308
309
310
311
312
313
314 k = 1
315 60 CONTINUE
316
317
318
319 IF( k.GT.n )
320 $ GO TO 80
321
322 IF( ipiv( k ).GT.0 ) THEN
323
324
325
326
327
328 kp = ipiv( k )
329 IF( kp.NE.k )
330 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
331
332
333
334
335 IF( k.LT.n )
336 $
CALL zgeru( n-k, nrhs, -one, a( k+1, k ), 1, b( k, 1 ),
337 $ ldb, b( k+1, 1 ), ldb )
338
339
340
341 CALL zscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
342 k = k + 1
343 ELSE
344
345
346
347
348
349 kp = -ipiv( k )
350 IF( kp.NE.k+1 )
351 $
CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
352
353
354
355
356 IF( k.LT.n-1 ) THEN
357 CALL zgeru( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
358 $ ldb, b( k+2, 1 ), ldb )
359 CALL zgeru( n-k-1, nrhs, -one, a( k+2, k+1 ), 1,
360 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
361 END IF
362
363
364
365 akm1k = a( k+1, k )
366 akm1 = a( k, k ) / akm1k
367 ak = a( k+1, k+1 ) / akm1k
368 denom = akm1*ak - one
369 DO 70 j = 1, nrhs
370 bkm1 = b( k, j ) / akm1k
371 bk = b( k+1, j ) / akm1k
372 b( k, j ) = ( ak*bkm1-bk ) / denom
373 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
374 70 CONTINUE
375 k = k + 2
376 END IF
377
378 GO TO 60
379 80 CONTINUE
380
381
382
383
384
385
386 k = n
387 90 CONTINUE
388
389
390
391 IF( k.LT.1 )
392 $ GO TO 100
393
394 IF( ipiv( k ).GT.0 ) THEN
395
396
397
398
399
400
401 IF( k.LT.n )
402 $
CALL zgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
403 $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
404
405
406
407 kp = ipiv( k )
408 IF( kp.NE.k )
409 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
410 k = k - 1
411 ELSE
412
413
414
415
416
417
418 IF( k.LT.n ) THEN
419 CALL zgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
420 $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
421 CALL zgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
422 $ ldb, a( k+1, k-1 ), 1, one, b( k-1, 1 ),
423 $ ldb )
424 END IF
425
426
427
428 kp = -ipiv( k )
429 IF( kp.NE.k )
430 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
431 k = k - 2
432 END IF
433
434 GO TO 90
435 100 CONTINUE
436 END IF
437
438 RETURN
439
440
441
subroutine xerbla(srname, info)
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
logical function lsame(ca, cb)
LSAME
subroutine zscal(n, za, zx, incx)
ZSCAL
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP