267 SUBROUTINE dlals0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
268 $ perm, givptr, givcol, ldgcol, givnum, ldgnum,
269 $ poles, difl, difr, z, k, c, s, work, info )
277 INTEGER givptr, icompq, info, k, ldb, ldbx, ldgcol,
278 $ ldgnum, nl, nr, nrhs, sqre
279 DOUBLE PRECISION c, s
282 INTEGER givcol( ldgcol, * ), perm( * )
283 DOUBLE PRECISION b( ldb, * ), bx( ldbx, * ), difl( * ),
284 $ difr( ldgnum, * ), givnum( ldgnum, * ),
285 $ poles( ldgnum, * ), work( * ), z( * )
291 DOUBLE PRECISION one, zero, negone
292 parameter( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
295 INTEGER i,
j, m, n, nlp1
296 DOUBLE PRECISION diflj, difrj, dj, dsigj, dsigjp, temp
315 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) )
THEN
317 ELSE IF( nl.LT.1 )
THEN
319 ELSE IF( nr.LT.1 )
THEN
321 ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) )
THEN
329 ELSE IF( ldb.LT.n )
THEN
331 ELSE IF( ldbx.LT.n )
THEN
333 ELSE IF( givptr.LT.0 )
THEN
335 ELSE IF( ldgcol.LT.n )
THEN
337 ELSE IF( ldgnum.LT.n )
THEN
339 ELSE IF( k.LT.1 )
THEN
343 CALL
xerbla(
'DLALS0', -info )
350 IF( icompq.EQ.0 )
THEN
357 CALL
drot( nrhs,
b( givcol( i, 2 ), 1 ), ldb,
358 $
b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
364 CALL
dcopy( nrhs,
b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
366 CALL
dcopy( nrhs,
b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
373 CALL
dcopy( nrhs, bx, ldbx,
b, ldb )
374 IF( z( 1 ).LT.zero )
THEN
375 CALL
dscal( nrhs, negone,
b, ldb )
381 dsigj = -poles(
j, 2 )
383 difrj = -difr(
j, 1 )
384 dsigjp = -poles(
j+1, 2 )
386 IF( ( z(
j ).EQ.zero ) .OR. ( poles(
j, 2 ).EQ.zero ) )
390 work(
j ) = -poles(
j, 2 )*z(
j ) / diflj /
391 $ ( poles(
j, 2 )+dj )
394 IF( ( z( i ).EQ.zero ) .OR.
395 $ ( poles( i, 2 ).EQ.zero ) )
THEN
398 work( i ) = poles( i, 2 )*z( i ) /
399 $ (
dlamc3( poles( i, 2 ), dsigj )-
400 $ diflj ) / ( poles( i, 2 )+dj )
404 IF( ( z( i ).EQ.zero ) .OR.
405 $ ( poles( i, 2 ).EQ.zero ) )
THEN
408 work( i ) = poles( i, 2 )*z( i ) /
409 $ (
dlamc3( poles( i, 2 ), dsigjp )+
410 $ difrj ) / ( poles( i, 2 )+dj )
414 temp =
dnrm2( k, work, 1 )
415 CALL
dgemv(
'T', k, nrhs, one, bx, ldbx, work, 1, zero,
417 CALL
dlascl(
'G', 0, 0, temp, one, 1, nrhs,
b(
j, 1 ),
424 IF( k.LT.max( m, n ) )
425 $ CALL
dlacpy(
'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
435 CALL
dcopy( nrhs,
b, ldb, bx, ldbx )
438 dsigj = poles(
j, 2 )
439 IF( z(
j ).EQ.zero )
THEN
442 work(
j ) = -z(
j ) / difl(
j ) /
443 $ ( dsigj+poles(
j, 1 ) ) / difr(
j, 2 )
446 IF( z(
j ).EQ.zero )
THEN
449 work( i ) = z(
j ) / (
dlamc3( dsigj, -poles( i+1,
450 $ 2 ) )-difr( i, 1 ) ) /
451 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
455 IF( z(
j ).EQ.zero )
THEN
458 work( i ) = z(
j ) / (
dlamc3( dsigj, -poles( i,
459 $ 2 ) )-difl( i ) ) /
460 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
463 CALL
dgemv(
'T', k, nrhs, one,
b, ldb, work, 1, zero,
472 CALL
dcopy( nrhs,
b( m, 1 ), ldb, bx( m, 1 ), ldbx )
473 CALL
drot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
475 IF( k.LT.max( m, n ) )
476 $ CALL
dlacpy(
'A', n-k, nrhs,
b( k+1, 1 ), ldb, bx( k+1, 1 ),
481 CALL
dcopy( nrhs, bx( 1, 1 ), ldbx,
b( nlp1, 1 ), ldb )
483 CALL
dcopy( nrhs, bx( m, 1 ), ldbx,
b( m, 1 ), ldb )
486 CALL
dcopy( nrhs, bx( i, 1 ), ldbx,
b( perm( i ), 1 ), ldb )
491 DO 100 i = givptr, 1, -1
492 CALL
drot( nrhs,
b( givcol( i, 2 ), 1 ), ldb,
493 $
b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function dlamc3(A, B)
DLAMC3
subroutine dlals0(ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO)
DLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dnrm2(N, X, INCX)
DNRM2
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV