233 SUBROUTINE zlatrs3( UPLO, TRANS, DIAG, NORMIN, N, NRHS, A, LDA,
234 $ X, LDX, SCALE, CNORM, WORK, LWORK, INFO )
238 CHARACTER DIAG, TRANS, NORMIN, UPLO
239 INTEGER INFO, LDA, LWORK, LDX, N, NRHS
242 COMPLEX*16 A( LDA, * ), X( LDX, * )
243 DOUBLE PRECISION CNORM( * ), SCALE( * ), WORK( * )
249 DOUBLE PRECISION ZERO, ONE
250 parameter( zero = 0.0d+0, one = 1.0d+0 )
251 COMPLEX*16 CZERO, CONE
252 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
253 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
254 INTEGER NBMAX, NBMIN, NBRHS, NRHSMIN
255 parameter( nrhsmin = 2, nbrhs = 32 )
256 parameter( nbmin = 8, nbmax = 64 )
259 DOUBLE PRECISION W( NBMAX ), XNRM( NBRHS )
262 LOGICAL LQUERY, NOTRAN, NOUNIT, UPPER
263 INTEGER AWRK, I, IFIRST, IINC, ILAST, II, I1, I2, J,
264 $ jfirst, jinc, jlast, j1, j2, k, kk, k1, k2,
265 $ lanrm, lds, lscale, nb, nba, nbx, rhs, lwmin
266 DOUBLE PRECISION ANRM, BIGNUM, BNRM, RSCAL, SCAL, SCALOC,
267 $ scamin, smlnum, tmax
272 DOUBLE PRECISION DLAMCH, ZLANGE, DLARMM
273 EXTERNAL ilaenv, lsame, dlamch, zlange,
280 INTRINSIC abs, max, min
285 upper = lsame( uplo,
'U' )
286 notran = lsame( trans,
'N' )
287 nounit = lsame( diag,
'N' )
288 lquery = ( lwork.EQ.-1 )
292 nb = max( nbmin, ilaenv( 1,
'ZLATRS',
'', n, n, -1, -1 ) )
293 nb = min( nbmax, nb )
294 nba = max( 1, (n + nb - 1) / nb )
295 nbx = max( 1, (nrhs + nbrhs - 1) / nbrhs )
306 lscale = nba * max( nba, min( nrhs, nbrhs ) )
317 IF( min( n, nrhs ).EQ.0 )
THEN
320 lwmin = lscale + lanrm
326 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
328 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
329 $ lsame( trans,
'C' ) )
THEN
331 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag,
'U' ) )
THEN
333 ELSE IF( .NOT.lsame( normin,
'Y' ) .AND. .NOT.
334 $ lsame( normin,
'N' ) )
THEN
336 ELSE IF( n.LT.0 )
THEN
338 ELSE IF( nrhs.LT.0 )
THEN
340 ELSE IF( lda.LT.max( 1, n ) )
THEN
342 ELSE IF( ldx.LT.max( 1, n ) )
THEN
344 ELSE IF( .NOT.lquery .AND. lwork.LT.lwmin )
THEN
348 CALL xerbla(
'ZLATRS3', -info )
350 ELSE IF( lquery )
THEN
362 IF( min( n, nrhs ).EQ.0 )
367 bignum = dlamch(
'Overflow' )
368 smlnum = dlamch(
'Safe Minimum' )
372 IF( nrhs.LT.nrhsmin )
THEN
373 CALL zlatrs( uplo, trans, diag, normin, n, a, lda, x( 1, 1),
374 $ scale( 1 ), cnorm, info )
376 CALL zlatrs( uplo, trans, diag,
'Y', n, a, lda, x( 1,
378 $ scale( k ), cnorm, info )
389 j2 = min( j*nb, n ) + 1
399 i2 = min( i*nb, n ) + 1
404 anrm = zlange(
'I', i2-i1, j2-j1, a( i1, j1 ), lda,
406 work( awrk + i+(j-1)*nba ) = anrm
408 anrm = zlange(
'1', i2-i1, j2-j1, a( i1, j1 ), lda,
410 work( awrk + j+(i-1) * nba ) = anrm
412 tmax = max( tmax, anrm )
416 IF( .NOT. tmax.LE.dlamch(
'Overflow') )
THEN
426 CALL zlatrs( uplo, trans, diag,
'N', n, a, lda, x( 1,
428 $ scale( k ), cnorm, info )
444 k2 = min( k*nbrhs, nrhs ) + 1
450 work( i+kk*lds ) = one
483 DO j = jfirst, jlast, jinc
488 j2 = min( j*nb, n ) + 1
495 CALL zlatrs( uplo, trans, diag,
'N', j2-j1,
496 $ a( j1, j1 ), lda, x( j1, rhs ),
497 $ scaloc, cnorm, info )
499 CALL zlatrs( uplo, trans, diag,
'Y', j2-j1,
500 $ a( j1, j1 ), lda, x( j1, rhs ),
501 $ scaloc, cnorm, info )
506 xnrm( kk ) = zlange(
'I', j2-j1, 1, x( j1, rhs ),
509 IF( scaloc .EQ. zero )
THEN
523 work( ii+kk*lds ) = one
526 ELSE IF( scaloc*work( j+kk*lds ) .EQ. zero )
THEN
533 scal = work( j+kk*lds ) / smlnum
534 scaloc = scaloc * scal
535 work( j+kk*lds ) = smlnum
540 IF( xnrm( kk )*rscal .LE. bignum )
THEN
541 xnrm( kk ) = xnrm( kk ) * rscal
542 CALL zdscal( j2-j1, rscal, x( j1, rhs ), 1 )
556 work( ii+kk*lds ) = one
561 scaloc = scaloc * work( j+kk*lds )
562 work( j+kk*lds ) = scaloc
589 DO i = ifirst, ilast, iinc
594 i2 = min( i*nb, n ) + 1
605 scamin = min( work( i+kk*lds), work( j+kk*lds ) )
610 bnrm = zlange(
'I', i2-i1, 1, x( i1, rhs ), ldx,
612 bnrm = bnrm*( scamin / work( i+kk*lds ) )
613 xnrm( kk ) = xnrm( kk )*( scamin / work( j+kk*lds) )
614 anrm = work( awrk + i+(j-1)*nba )
615 scaloc = dlarmm( anrm, xnrm( kk ), bnrm )
620 scal = ( scamin / work( i+kk*lds) )*scaloc
621 IF( scal.NE.one )
THEN
622 CALL zdscal( i2-i1, scal, x( i1, rhs ), 1 )
623 work( i+kk*lds ) = scamin*scaloc
626 scal = ( scamin / work( j+kk*lds ) )*scaloc
627 IF( scal.NE.one )
THEN
628 CALL zdscal( j2-j1, scal, x( j1, rhs ), 1 )
629 work( j+kk*lds ) = scamin*scaloc
637 CALL zgemm(
'N',
'N', i2-i1, k2-k1, j2-j1, -cone,
638 $ a( i1, j1 ), lda, x( j1, k1 ), ldx,
639 $ cone, x( i1, k1 ), ldx )
640 ELSE IF( lsame( trans,
'T' ) )
THEN
644 CALL zgemm(
'T',
'N', i2-i1, k2-k1, j2-j1, -cone,
645 $ a( j1, i1 ), lda, x( j1, k1 ), ldx,
646 $ cone, x( i1, k1 ), ldx )
651 CALL zgemm(
'C',
'N', i2-i1, k2-k1, j2-j1, -cone,
652 $ a( j1, i1 ), lda, x( j1, k1 ), ldx,
653 $ cone, x( i1, k1 ), ldx )
664 scale( rhs ) = min( scale( rhs ), work( i+kk*lds ) )
672 IF( scale( rhs ).NE.one .AND. scale( rhs ).NE. zero )
THEN
674 i1 = (i - 1) * nb + 1
675 i2 = min( i * nb, n ) + 1
676 scal = scale( rhs ) / work( i+kk*lds )
678 $
CALL zdscal( i2-i1, scal, x( i1, rhs ), 1 )