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