136
137
138
139
140
141
142 CHARACTER UPLO
143 INTEGER INFO, LDA, LDB, N, NRHS
144
145
146 INTEGER IPIV( * )
147 REAL A( LDA, * ), B( LDB, * )
148
149
150
151
152
153 REAL ONE
154 parameter( one = 1.0e+0 )
155
156
157 LOGICAL UPPER
158 INTEGER J, K, KP
159 REAL AK, AKM1, AKM1K, BK, BKM1, DENOM
160
161
162 LOGICAL LSAME
164
165
167
168
169 INTRINSIC max
170
171
172
173 info = 0
174 upper =
lsame( uplo,
'U' )
175 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
176 info = -1
177 ELSE IF( n.LT.0 ) THEN
178 info = -2
179 ELSE IF( nrhs.LT.0 ) THEN
180 info = -3
181 ELSE IF( lda.LT.max( 1, n ) ) THEN
182 info = -5
183 ELSE IF( ldb.LT.max( 1, n ) ) THEN
184 info = -8
185 END IF
186 IF( info.NE.0 ) THEN
187 CALL xerbla(
'SSYTRS_ROOK', -info )
188 RETURN
189 END IF
190
191
192
193 IF( n.EQ.0 .OR. nrhs.EQ.0 )
194 $ RETURN
195
196 IF( upper ) THEN
197
198
199
200
201
202
203
204
205 k = n
206 10 CONTINUE
207
208
209
210 IF( k.LT.1 )
211 $ GO TO 30
212
213 IF( ipiv( k ).GT.0 ) THEN
214
215
216
217
218
219 kp = ipiv( k )
220 IF( kp.NE.k )
221 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
222
223
224
225
226 CALL sger( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
227 $ b( 1, 1 ), ldb )
228
229
230
231 CALL sscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
232 k = k - 1
233 ELSE
234
235
236
237
238
239 kp = -ipiv( k )
240 IF( kp.NE.k )
241 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
242
243 kp = -ipiv( k-1 )
244 IF( kp.NE.k-1 )
245 $
CALL sswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
246
247
248
249
250 IF( k.GT.2 ) THEN
251 CALL sger( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ),
252 $ ldb, b( 1, 1 ), ldb )
253 CALL sger( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 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 - one
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 sgemv(
'Transpose', k-1, nrhs, -one, b,
297 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
298
299
300
301 kp = ipiv( k )
302 IF( kp.NE.k )
303 $
CALL sswap( 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 sgemv(
'Transpose', k-1, nrhs, -one, b,
314 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
315 CALL sgemv(
'Transpose', k-1, nrhs, -one, b,
316 $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
317 END IF
318
319
320
321 kp = -ipiv( k )
322 IF( kp.NE.k )
323 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
324
325 kp = -ipiv( k+1 )
326 IF( kp.NE.k+1 )
327 $
CALL sswap( 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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
361
362
363
364
365 IF( k.LT.n )
366 $
CALL sger( n-k, nrhs, -one, a( k+1, k ), 1, b( k, 1 ),
367 $ ldb, b( k+1, 1 ), ldb )
368
369
370
371 CALL sscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
372 k = k + 1
373 ELSE
374
375
376
377
378
379 kp = -ipiv( k )
380 IF( kp.NE.k )
381 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
382
383 kp = -ipiv( k+1 )
384 IF( kp.NE.k+1 )
385 $
CALL sswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
386
387
388
389
390 IF( k.LT.n-1 ) THEN
391 CALL sger( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
392 $ ldb, b( k+2, 1 ), ldb )
393 CALL sger( n-k-1, nrhs, -one, a( k+2, k+1 ), 1,
394 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
395 END IF
396
397
398
399 akm1k = a( k+1, k )
400 akm1 = a( k, k ) / akm1k
401 ak = a( k+1, k+1 ) / akm1k
402 denom = akm1*ak - one
403 DO 70 j = 1, nrhs
404 bkm1 = b( k, j ) / akm1k
405 bk = b( k+1, j ) / akm1k
406 b( k, j ) = ( ak*bkm1-bk ) / denom
407 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
408 70 CONTINUE
409 k = k + 2
410 END IF
411
412 GO TO 60
413 80 CONTINUE
414
415
416
417
418
419
420 k = n
421 90 CONTINUE
422
423
424
425 IF( k.LT.1 )
426 $ GO TO 100
427
428 IF( ipiv( k ).GT.0 ) THEN
429
430
431
432
433
434
435 IF( k.LT.n )
436 $
CALL sgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
437 $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
438
439
440
441 kp = ipiv( k )
442 IF( kp.NE.k )
443 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
444 k = k - 1
445 ELSE
446
447
448
449
450
451
452 IF( k.LT.n ) THEN
453 CALL sgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
454 $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
455 CALL sgemv(
'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
456 $ ldb, a( k+1, k-1 ), 1, one, b( k-1, 1 ),
457 $ ldb )
458 END IF
459
460
461
462 kp = -ipiv( k )
463 IF( kp.NE.k )
464 $
CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
465
466 kp = -ipiv( k-1 )
467 IF( kp.NE.k-1 )
468 $
CALL sswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
469
470 k = k - 2
471 END IF
472
473 GO TO 90
474 100 CONTINUE
475 END IF
476
477 RETURN
478
479
480
subroutine xerbla(srname, info)
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
subroutine sger(m, n, alpha, x, incx, y, incy, a, lda)
SGER
logical function lsame(ca, cb)
LSAME
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine sswap(n, sx, incx, sy, incy)
SSWAP