LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dtgsy2()

subroutine dtgsy2 ( character  trans,
integer  ijob,
integer  m,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( ldc, * )  c,
integer  ldc,
double precision, dimension( ldd, * )  d,
integer  ldd,
double precision, dimension( lde, * )  e,
integer  lde,
double precision, dimension( ldf, * )  f,
integer  ldf,
double precision  scale,
double precision  rdsum,
double precision  rdscal,
integer, dimension( * )  iwork,
integer  pq,
integer  info 
)

DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Download DTGSY2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DTGSY2 solves the generalized Sylvester equation:

             A * R - L * B = scale * C                (1)
             D * R - L * E = scale * F,

 using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
 (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
 N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
 must be in generalized Schur canonical form, i.e. A, B are upper
 quasi triangular and D, E are upper triangular. The solution (R, L)
 overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
 chosen to avoid overflow.

 In matrix notation solving equation (1) corresponds to solve
 Z*x = scale*b, where Z is defined as

        Z = [ kron(In, A)  -kron(B**T, Im) ]             (2)
            [ kron(In, D)  -kron(E**T, Im) ],

 Ik is the identity matrix of size k and X**T is the transpose of X.
 kron(X, Y) is the Kronecker product between the matrices X and Y.
 In the process of solving (1), we solve a number of such systems
 where Dim(In), Dim(In) = 1 or 2.

 If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
 which is equivalent to solve for R and L in

             A**T * R  + D**T * L   = scale * C           (3)
             R  * B**T + L  * E**T  = scale * -F

 This case is used to compute an estimate of Dif[(A, D), (B, E)] =
 sigma_min(Z) using reverse communication with DLACON.

 DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
 of an upper bound on the separation between to matrix pairs. Then
 the input (A, D), (B, E) are sub-pencils of the matrix pair in
 DTGSYL. See DTGSYL for details.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': solve the generalized Sylvester equation (1).
          = 'T': solve the 'transposed' system (3).
[in]IJOB
          IJOB is INTEGER
          Specifies what kind of functionality to be performed.
          = 0: solve (1) only.
          = 1: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (look ahead strategy is used).
          = 2: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (DGECON on sub-systems is used.)
          Not referenced if TRANS = 'T'.
[in]M
          M is INTEGER
          On entry, M specifies the order of A and D, and the row
          dimension of C, F, R and L.
[in]N
          N is INTEGER
          On entry, N specifies the order of B and E, and the column
          dimension of C, F, R and L.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA, M)
          On entry, A contains an upper quasi triangular matrix.
[in]LDA
          LDA is INTEGER
          The leading dimension of the matrix A. LDA >= max(1, M).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB, N)
          On entry, B contains an upper quasi triangular matrix.
[in]LDB
          LDB is INTEGER
          The leading dimension of the matrix B. LDB >= max(1, N).
[in,out]C
          C is DOUBLE PRECISION array, dimension (LDC, N)
          On entry, C contains the right-hand-side of the first matrix
          equation in (1).
          On exit, if IJOB = 0, C has been overwritten by the
          solution R.
[in]LDC
          LDC is INTEGER
          The leading dimension of the matrix C. LDC >= max(1, M).
[in]D
          D is DOUBLE PRECISION array, dimension (LDD, M)
          On entry, D contains an upper triangular matrix.
[in]LDD
          LDD is INTEGER
          The leading dimension of the matrix D. LDD >= max(1, M).
[in]E
          E is DOUBLE PRECISION array, dimension (LDE, N)
          On entry, E contains an upper triangular matrix.
[in]LDE
          LDE is INTEGER
          The leading dimension of the matrix E. LDE >= max(1, N).
[in,out]F
          F is DOUBLE PRECISION array, dimension (LDF, N)
          On entry, F contains the right-hand-side of the second matrix
          equation in (1).
          On exit, if IJOB = 0, F has been overwritten by the
          solution L.
[in]LDF
          LDF is INTEGER
          The leading dimension of the matrix F. LDF >= max(1, M).
[out]SCALE
          SCALE is DOUBLE PRECISION
          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
          R and L (C and F on entry) will hold the solutions to a
          slightly perturbed system but the input matrices A, B, D and
          E have not been changed. If SCALE = 0, R and L will hold the
          solutions to the homogeneous system with C = F = 0. Normally,
          SCALE = 1.
[in,out]RDSUM
          RDSUM is DOUBLE PRECISION
          On entry, the sum of squares of computed contributions to
          the Dif-estimate under computation by DTGSYL, where the
          scaling factor RDSCAL (see below) has been factored out.
          On exit, the corresponding sum of squares updated with the
          contributions from the current sub-system.
          If TRANS = 'T' RDSUM is not touched.
          NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
[in,out]RDSCAL
          RDSCAL is DOUBLE PRECISION
          On entry, scaling factor used to prevent overflow in RDSUM.
          On exit, RDSCAL is updated w.r.t. the current contributions
          in RDSUM.
          If TRANS = 'T', RDSCAL is not touched.
          NOTE: RDSCAL only makes sense when DTGSY2 is called by
                DTGSYL.
[out]IWORK
          IWORK is INTEGER array, dimension (M+N+2)
[out]PQ
          PQ is INTEGER
          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
          8-by-8) solved by this routine.
[out]INFO
          INFO is INTEGER
          On exit, if INFO is set to
            =0: Successful exit
            <0: If INFO = -i, the i-th argument had an illegal value.
            >0: The matrix pairs (A, D) and (B, E) have common or very
                close eigenvalues.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 271 of file dtgsy2.f.

274*
275* -- LAPACK auxiliary routine --
276* -- LAPACK is a software package provided by Univ. of Tennessee, --
277* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
278*
279* .. Scalar Arguments ..
280 CHARACTER TRANS
281 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
282 $ PQ
283 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
284* ..
285* .. Array Arguments ..
286 INTEGER IWORK( * )
287 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
288 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
289* ..
290*
291* =====================================================================
292* Replaced various illegal calls to DCOPY by calls to DLASET.
293* Sven Hammarling, 27/5/02.
294*
295* .. Parameters ..
296 INTEGER LDZ
297 parameter( ldz = 8 )
298 DOUBLE PRECISION ZERO, ONE
299 parameter( zero = 0.0d+0, one = 1.0d+0 )
300* ..
301* .. Local Scalars ..
302 LOGICAL NOTRAN
303 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
304 $ K, MB, NB, P, Q, ZDIM
305 DOUBLE PRECISION ALPHA, SCALOC
306* ..
307* .. Local Arrays ..
308 INTEGER IPIV( LDZ ), JPIV( LDZ )
309 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
310* ..
311* .. External Functions ..
312 LOGICAL LSAME
313 EXTERNAL lsame
314* ..
315* .. External Subroutines ..
316 EXTERNAL daxpy, dcopy, dgemm, dgemv, dger, dgesc2,
317 $ dgetc2, dlaset, dlatdf, dscal, xerbla
318* ..
319* .. Intrinsic Functions ..
320 INTRINSIC max
321* ..
322* .. Executable Statements ..
323*
324* Decode and test input parameters
325*
326 info = 0
327 ierr = 0
328 notran = lsame( trans, 'N' )
329 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
330 info = -1
331 ELSE IF( notran ) THEN
332 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
333 info = -2
334 END IF
335 END IF
336 IF( info.EQ.0 ) THEN
337 IF( m.LE.0 ) THEN
338 info = -3
339 ELSE IF( n.LE.0 ) THEN
340 info = -4
341 ELSE IF( lda.LT.max( 1, m ) ) THEN
342 info = -6
343 ELSE IF( ldb.LT.max( 1, n ) ) THEN
344 info = -8
345 ELSE IF( ldc.LT.max( 1, m ) ) THEN
346 info = -10
347 ELSE IF( ldd.LT.max( 1, m ) ) THEN
348 info = -12
349 ELSE IF( lde.LT.max( 1, n ) ) THEN
350 info = -14
351 ELSE IF( ldf.LT.max( 1, m ) ) THEN
352 info = -16
353 END IF
354 END IF
355 IF( info.NE.0 ) THEN
356 CALL xerbla( 'DTGSY2', -info )
357 RETURN
358 END IF
359*
360* Determine block structure of A
361*
362 pq = 0
363 p = 0
364 i = 1
365 10 CONTINUE
366 IF( i.GT.m )
367 $ GO TO 20
368 p = p + 1
369 iwork( p ) = i
370 IF( i.EQ.m )
371 $ GO TO 20
372 IF( a( i+1, i ).NE.zero ) THEN
373 i = i + 2
374 ELSE
375 i = i + 1
376 END IF
377 GO TO 10
378 20 CONTINUE
379 iwork( p+1 ) = m + 1
380*
381* Determine block structure of B
382*
383 q = p + 1
384 j = 1
385 30 CONTINUE
386 IF( j.GT.n )
387 $ GO TO 40
388 q = q + 1
389 iwork( q ) = j
390 IF( j.EQ.n )
391 $ GO TO 40
392 IF( b( j+1, j ).NE.zero ) THEN
393 j = j + 2
394 ELSE
395 j = j + 1
396 END IF
397 GO TO 30
398 40 CONTINUE
399 iwork( q+1 ) = n + 1
400 pq = p*( q-p-1 )
401*
402 IF( notran ) THEN
403*
404* Solve (I, J) - subsystem
405* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
406* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
407* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
408*
409 scale = one
410 scaloc = one
411 DO 120 j = p + 2, q
412 js = iwork( j )
413 jsp1 = js + 1
414 je = iwork( j+1 ) - 1
415 nb = je - js + 1
416 DO 110 i = p, 1, -1
417*
418 is = iwork( i )
419 isp1 = is + 1
420 ie = iwork( i+1 ) - 1
421 mb = ie - is + 1
422 zdim = mb*nb*2
423*
424 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
425*
426* Build a 2-by-2 system Z * x = RHS
427*
428 z( 1, 1 ) = a( is, is )
429 z( 2, 1 ) = d( is, is )
430 z( 1, 2 ) = -b( js, js )
431 z( 2, 2 ) = -e( js, js )
432*
433* Set up right hand side(s)
434*
435 rhs( 1 ) = c( is, js )
436 rhs( 2 ) = f( is, js )
437*
438* Solve Z * x = RHS
439*
440 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441 IF( ierr.GT.0 )
442 $ info = ierr
443*
444 IF( ijob.EQ.0 ) THEN
445 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446 $ scaloc )
447 IF( scaloc.NE.one ) THEN
448 DO 50 k = 1, n
449 CALL dscal( m, scaloc, c( 1, k ), 1 )
450 CALL dscal( m, scaloc, f( 1, k ), 1 )
451 50 CONTINUE
452 scale = scale*scaloc
453 END IF
454 ELSE
455 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
456 $ rdscal, ipiv, jpiv )
457 END IF
458*
459* Unpack solution vector(s)
460*
461 c( is, js ) = rhs( 1 )
462 f( is, js ) = rhs( 2 )
463*
464* Substitute R(I, J) and L(I, J) into remaining
465* equation.
466*
467 IF( i.GT.1 ) THEN
468 alpha = -rhs( 1 )
469 CALL daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
470 $ 1 )
471 CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
472 $ 1 )
473 END IF
474 IF( j.LT.q ) THEN
475 CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
476 $ c( is, je+1 ), ldc )
477 CALL daxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
478 $ f( is, je+1 ), ldf )
479 END IF
480*
481 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
482*
483* Build a 4-by-4 system Z * x = RHS
484*
485 z( 1, 1 ) = a( is, is )
486 z( 2, 1 ) = zero
487 z( 3, 1 ) = d( is, is )
488 z( 4, 1 ) = zero
489*
490 z( 1, 2 ) = zero
491 z( 2, 2 ) = a( is, is )
492 z( 3, 2 ) = zero
493 z( 4, 2 ) = d( is, is )
494*
495 z( 1, 3 ) = -b( js, js )
496 z( 2, 3 ) = -b( js, jsp1 )
497 z( 3, 3 ) = -e( js, js )
498 z( 4, 3 ) = -e( js, jsp1 )
499*
500 z( 1, 4 ) = -b( jsp1, js )
501 z( 2, 4 ) = -b( jsp1, jsp1 )
502 z( 3, 4 ) = zero
503 z( 4, 4 ) = -e( jsp1, jsp1 )
504*
505* Set up right hand side(s)
506*
507 rhs( 1 ) = c( is, js )
508 rhs( 2 ) = c( is, jsp1 )
509 rhs( 3 ) = f( is, js )
510 rhs( 4 ) = f( is, jsp1 )
511*
512* Solve Z * x = RHS
513*
514 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
515 IF( ierr.GT.0 )
516 $ info = ierr
517*
518 IF( ijob.EQ.0 ) THEN
519 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
520 $ scaloc )
521 IF( scaloc.NE.one ) THEN
522 DO 60 k = 1, n
523 CALL dscal( m, scaloc, c( 1, k ), 1 )
524 CALL dscal( m, scaloc, f( 1, k ), 1 )
525 60 CONTINUE
526 scale = scale*scaloc
527 END IF
528 ELSE
529 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
530 $ rdscal, ipiv, jpiv )
531 END IF
532*
533* Unpack solution vector(s)
534*
535 c( is, js ) = rhs( 1 )
536 c( is, jsp1 ) = rhs( 2 )
537 f( is, js ) = rhs( 3 )
538 f( is, jsp1 ) = rhs( 4 )
539*
540* Substitute R(I, J) and L(I, J) into remaining
541* equation.
542*
543 IF( i.GT.1 ) THEN
544 CALL dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
545 $ 1, c( 1, js ), ldc )
546 CALL dger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
547 $ 1, f( 1, js ), ldf )
548 END IF
549 IF( j.LT.q ) THEN
550 CALL daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
551 $ c( is, je+1 ), ldc )
552 CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
553 $ f( is, je+1 ), ldf )
554 CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
555 $ c( is, je+1 ), ldc )
556 CALL daxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
557 $ f( is, je+1 ), ldf )
558 END IF
559*
560 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
561*
562* Build a 4-by-4 system Z * x = RHS
563*
564 z( 1, 1 ) = a( is, is )
565 z( 2, 1 ) = a( isp1, is )
566 z( 3, 1 ) = d( is, is )
567 z( 4, 1 ) = zero
568*
569 z( 1, 2 ) = a( is, isp1 )
570 z( 2, 2 ) = a( isp1, isp1 )
571 z( 3, 2 ) = d( is, isp1 )
572 z( 4, 2 ) = d( isp1, isp1 )
573*
574 z( 1, 3 ) = -b( js, js )
575 z( 2, 3 ) = zero
576 z( 3, 3 ) = -e( js, js )
577 z( 4, 3 ) = zero
578*
579 z( 1, 4 ) = zero
580 z( 2, 4 ) = -b( js, js )
581 z( 3, 4 ) = zero
582 z( 4, 4 ) = -e( js, js )
583*
584* Set up right hand side(s)
585*
586 rhs( 1 ) = c( is, js )
587 rhs( 2 ) = c( isp1, js )
588 rhs( 3 ) = f( is, js )
589 rhs( 4 ) = f( isp1, js )
590*
591* Solve Z * x = RHS
592*
593 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
594 IF( ierr.GT.0 )
595 $ info = ierr
596 IF( ijob.EQ.0 ) THEN
597 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
598 $ scaloc )
599 IF( scaloc.NE.one ) THEN
600 DO 70 k = 1, n
601 CALL dscal( m, scaloc, c( 1, k ), 1 )
602 CALL dscal( m, scaloc, f( 1, k ), 1 )
603 70 CONTINUE
604 scale = scale*scaloc
605 END IF
606 ELSE
607 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
608 $ rdscal, ipiv, jpiv )
609 END IF
610*
611* Unpack solution vector(s)
612*
613 c( is, js ) = rhs( 1 )
614 c( isp1, js ) = rhs( 2 )
615 f( is, js ) = rhs( 3 )
616 f( isp1, js ) = rhs( 4 )
617*
618* Substitute R(I, J) and L(I, J) into remaining
619* equation.
620*
621 IF( i.GT.1 ) THEN
622 CALL dgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
623 $ rhs( 1 ), 1, one, c( 1, js ), 1 )
624 CALL dgemv( 'N', is-1, mb, -one, d( 1, is ), ldd,
625 $ rhs( 1 ), 1, one, f( 1, js ), 1 )
626 END IF
627 IF( j.LT.q ) THEN
628 CALL dger( mb, n-je, one, rhs( 3 ), 1,
629 $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
630 CALL dger( mb, n-je, one, rhs( 3 ), 1,
631 $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
632 END IF
633*
634 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
635*
636* Build an 8-by-8 system Z * x = RHS
637*
638 CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
639*
640 z( 1, 1 ) = a( is, is )
641 z( 2, 1 ) = a( isp1, is )
642 z( 5, 1 ) = d( is, is )
643*
644 z( 1, 2 ) = a( is, isp1 )
645 z( 2, 2 ) = a( isp1, isp1 )
646 z( 5, 2 ) = d( is, isp1 )
647 z( 6, 2 ) = d( isp1, isp1 )
648*
649 z( 3, 3 ) = a( is, is )
650 z( 4, 3 ) = a( isp1, is )
651 z( 7, 3 ) = d( is, is )
652*
653 z( 3, 4 ) = a( is, isp1 )
654 z( 4, 4 ) = a( isp1, isp1 )
655 z( 7, 4 ) = d( is, isp1 )
656 z( 8, 4 ) = d( isp1, isp1 )
657*
658 z( 1, 5 ) = -b( js, js )
659 z( 3, 5 ) = -b( js, jsp1 )
660 z( 5, 5 ) = -e( js, js )
661 z( 7, 5 ) = -e( js, jsp1 )
662*
663 z( 2, 6 ) = -b( js, js )
664 z( 4, 6 ) = -b( js, jsp1 )
665 z( 6, 6 ) = -e( js, js )
666 z( 8, 6 ) = -e( js, jsp1 )
667*
668 z( 1, 7 ) = -b( jsp1, js )
669 z( 3, 7 ) = -b( jsp1, jsp1 )
670 z( 7, 7 ) = -e( jsp1, jsp1 )
671*
672 z( 2, 8 ) = -b( jsp1, js )
673 z( 4, 8 ) = -b( jsp1, jsp1 )
674 z( 8, 8 ) = -e( jsp1, jsp1 )
675*
676* Set up right hand side(s)
677*
678 k = 1
679 ii = mb*nb + 1
680 DO 80 jj = 0, nb - 1
681 CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
682 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
683 k = k + mb
684 ii = ii + mb
685 80 CONTINUE
686*
687* Solve Z * x = RHS
688*
689 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
690 IF( ierr.GT.0 )
691 $ info = ierr
692 IF( ijob.EQ.0 ) THEN
693 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
694 $ scaloc )
695 IF( scaloc.NE.one ) THEN
696 DO 90 k = 1, n
697 CALL dscal( m, scaloc, c( 1, k ), 1 )
698 CALL dscal( m, scaloc, f( 1, k ), 1 )
699 90 CONTINUE
700 scale = scale*scaloc
701 END IF
702 ELSE
703 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
704 $ rdscal, ipiv, jpiv )
705 END IF
706*
707* Unpack solution vector(s)
708*
709 k = 1
710 ii = mb*nb + 1
711 DO 100 jj = 0, nb - 1
712 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
713 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
714 k = k + mb
715 ii = ii + mb
716 100 CONTINUE
717*
718* Substitute R(I, J) and L(I, J) into remaining
719* equation.
720*
721 IF( i.GT.1 ) THEN
722 CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
723 $ a( 1, is ), lda, rhs( 1 ), mb, one,
724 $ c( 1, js ), ldc )
725 CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
726 $ d( 1, is ), ldd, rhs( 1 ), mb, one,
727 $ f( 1, js ), ldf )
728 END IF
729 IF( j.LT.q ) THEN
730 k = mb*nb + 1
731 CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
732 $ mb, b( js, je+1 ), ldb, one,
733 $ c( is, je+1 ), ldc )
734 CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
735 $ mb, e( js, je+1 ), lde, one,
736 $ f( is, je+1 ), ldf )
737 END IF
738*
739 END IF
740*
741 110 CONTINUE
742 120 CONTINUE
743 ELSE
744*
745* Solve (I, J) - subsystem
746* A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
747* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
748* for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
749*
750 scale = one
751 scaloc = one
752 DO 200 i = 1, p
753*
754 is = iwork( i )
755 isp1 = is + 1
756 ie = iwork( i+1 ) - 1
757 mb = ie - is + 1
758 DO 190 j = q, p + 2, -1
759*
760 js = iwork( j )
761 jsp1 = js + 1
762 je = iwork( j+1 ) - 1
763 nb = je - js + 1
764 zdim = mb*nb*2
765 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
766*
767* Build a 2-by-2 system Z**T * x = RHS
768*
769 z( 1, 1 ) = a( is, is )
770 z( 2, 1 ) = -b( js, js )
771 z( 1, 2 ) = d( is, is )
772 z( 2, 2 ) = -e( js, js )
773*
774* Set up right hand side(s)
775*
776 rhs( 1 ) = c( is, js )
777 rhs( 2 ) = f( is, js )
778*
779* Solve Z**T * x = RHS
780*
781 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
782 IF( ierr.GT.0 )
783 $ info = ierr
784*
785 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
786 IF( scaloc.NE.one ) THEN
787 DO 130 k = 1, n
788 CALL dscal( m, scaloc, c( 1, k ), 1 )
789 CALL dscal( m, scaloc, f( 1, k ), 1 )
790 130 CONTINUE
791 scale = scale*scaloc
792 END IF
793*
794* Unpack solution vector(s)
795*
796 c( is, js ) = rhs( 1 )
797 f( is, js ) = rhs( 2 )
798*
799* Substitute R(I, J) and L(I, J) into remaining
800* equation.
801*
802 IF( j.GT.p+2 ) THEN
803 alpha = rhs( 1 )
804 CALL daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
805 $ ldf )
806 alpha = rhs( 2 )
807 CALL daxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
808 $ ldf )
809 END IF
810 IF( i.LT.p ) THEN
811 alpha = -rhs( 1 )
812 CALL daxpy( m-ie, alpha, a( is, ie+1 ), lda,
813 $ c( ie+1, js ), 1 )
814 alpha = -rhs( 2 )
815 CALL daxpy( m-ie, alpha, d( is, ie+1 ), ldd,
816 $ c( ie+1, js ), 1 )
817 END IF
818*
819 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
820*
821* Build a 4-by-4 system Z**T * x = RHS
822*
823 z( 1, 1 ) = a( is, is )
824 z( 2, 1 ) = zero
825 z( 3, 1 ) = -b( js, js )
826 z( 4, 1 ) = -b( jsp1, js )
827*
828 z( 1, 2 ) = zero
829 z( 2, 2 ) = a( is, is )
830 z( 3, 2 ) = -b( js, jsp1 )
831 z( 4, 2 ) = -b( jsp1, jsp1 )
832*
833 z( 1, 3 ) = d( is, is )
834 z( 2, 3 ) = zero
835 z( 3, 3 ) = -e( js, js )
836 z( 4, 3 ) = zero
837*
838 z( 1, 4 ) = zero
839 z( 2, 4 ) = d( is, is )
840 z( 3, 4 ) = -e( js, jsp1 )
841 z( 4, 4 ) = -e( jsp1, jsp1 )
842*
843* Set up right hand side(s)
844*
845 rhs( 1 ) = c( is, js )
846 rhs( 2 ) = c( is, jsp1 )
847 rhs( 3 ) = f( is, js )
848 rhs( 4 ) = f( is, jsp1 )
849*
850* Solve Z**T * x = RHS
851*
852 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
853 IF( ierr.GT.0 )
854 $ info = ierr
855 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
856 IF( scaloc.NE.one ) THEN
857 DO 140 k = 1, n
858 CALL dscal( m, scaloc, c( 1, k ), 1 )
859 CALL dscal( m, scaloc, f( 1, k ), 1 )
860 140 CONTINUE
861 scale = scale*scaloc
862 END IF
863*
864* Unpack solution vector(s)
865*
866 c( is, js ) = rhs( 1 )
867 c( is, jsp1 ) = rhs( 2 )
868 f( is, js ) = rhs( 3 )
869 f( is, jsp1 ) = rhs( 4 )
870*
871* Substitute R(I, J) and L(I, J) into remaining
872* equation.
873*
874 IF( j.GT.p+2 ) THEN
875 CALL daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
876 $ f( is, 1 ), ldf )
877 CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
878 $ f( is, 1 ), ldf )
879 CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
880 $ f( is, 1 ), ldf )
881 CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
882 $ f( is, 1 ), ldf )
883 END IF
884 IF( i.LT.p ) THEN
885 CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
886 $ rhs( 1 ), 1, c( ie+1, js ), ldc )
887 CALL dger( m-ie, nb, -one, d( is, ie+1 ), ldd,
888 $ rhs( 3 ), 1, c( ie+1, js ), ldc )
889 END IF
890*
891 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
892*
893* Build a 4-by-4 system Z**T * x = RHS
894*
895 z( 1, 1 ) = a( is, is )
896 z( 2, 1 ) = a( is, isp1 )
897 z( 3, 1 ) = -b( js, js )
898 z( 4, 1 ) = zero
899*
900 z( 1, 2 ) = a( isp1, is )
901 z( 2, 2 ) = a( isp1, isp1 )
902 z( 3, 2 ) = zero
903 z( 4, 2 ) = -b( js, js )
904*
905 z( 1, 3 ) = d( is, is )
906 z( 2, 3 ) = d( is, isp1 )
907 z( 3, 3 ) = -e( js, js )
908 z( 4, 3 ) = zero
909*
910 z( 1, 4 ) = zero
911 z( 2, 4 ) = d( isp1, isp1 )
912 z( 3, 4 ) = zero
913 z( 4, 4 ) = -e( js, js )
914*
915* Set up right hand side(s)
916*
917 rhs( 1 ) = c( is, js )
918 rhs( 2 ) = c( isp1, js )
919 rhs( 3 ) = f( is, js )
920 rhs( 4 ) = f( isp1, js )
921*
922* Solve Z**T * x = RHS
923*
924 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
925 IF( ierr.GT.0 )
926 $ info = ierr
927*
928 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
929 IF( scaloc.NE.one ) THEN
930 DO 150 k = 1, n
931 CALL dscal( m, scaloc, c( 1, k ), 1 )
932 CALL dscal( m, scaloc, f( 1, k ), 1 )
933 150 CONTINUE
934 scale = scale*scaloc
935 END IF
936*
937* Unpack solution vector(s)
938*
939 c( is, js ) = rhs( 1 )
940 c( isp1, js ) = rhs( 2 )
941 f( is, js ) = rhs( 3 )
942 f( isp1, js ) = rhs( 4 )
943*
944* Substitute R(I, J) and L(I, J) into remaining
945* equation.
946*
947 IF( j.GT.p+2 ) THEN
948 CALL dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
949 $ 1, f( is, 1 ), ldf )
950 CALL dger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
951 $ 1, f( is, 1 ), ldf )
952 END IF
953 IF( i.LT.p ) THEN
954 CALL dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
955 $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
956 $ 1 )
957 CALL dgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
958 $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
959 $ 1 )
960 END IF
961*
962 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
963*
964* Build an 8-by-8 system Z**T * x = RHS
965*
966 CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
967*
968 z( 1, 1 ) = a( is, is )
969 z( 2, 1 ) = a( is, isp1 )
970 z( 5, 1 ) = -b( js, js )
971 z( 7, 1 ) = -b( jsp1, js )
972*
973 z( 1, 2 ) = a( isp1, is )
974 z( 2, 2 ) = a( isp1, isp1 )
975 z( 6, 2 ) = -b( js, js )
976 z( 8, 2 ) = -b( jsp1, js )
977*
978 z( 3, 3 ) = a( is, is )
979 z( 4, 3 ) = a( is, isp1 )
980 z( 5, 3 ) = -b( js, jsp1 )
981 z( 7, 3 ) = -b( jsp1, jsp1 )
982*
983 z( 3, 4 ) = a( isp1, is )
984 z( 4, 4 ) = a( isp1, isp1 )
985 z( 6, 4 ) = -b( js, jsp1 )
986 z( 8, 4 ) = -b( jsp1, jsp1 )
987*
988 z( 1, 5 ) = d( is, is )
989 z( 2, 5 ) = d( is, isp1 )
990 z( 5, 5 ) = -e( js, js )
991*
992 z( 2, 6 ) = d( isp1, isp1 )
993 z( 6, 6 ) = -e( js, js )
994*
995 z( 3, 7 ) = d( is, is )
996 z( 4, 7 ) = d( is, isp1 )
997 z( 5, 7 ) = -e( js, jsp1 )
998 z( 7, 7 ) = -e( jsp1, jsp1 )
999*
1000 z( 4, 8 ) = d( isp1, isp1 )
1001 z( 6, 8 ) = -e( js, jsp1 )
1002 z( 8, 8 ) = -e( jsp1, jsp1 )
1003*
1004* Set up right hand side(s)
1005*
1006 k = 1
1007 ii = mb*nb + 1
1008 DO 160 jj = 0, nb - 1
1009 CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1010 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1011 k = k + mb
1012 ii = ii + mb
1013 160 CONTINUE
1014*
1015*
1016* Solve Z**T * x = RHS
1017*
1018 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1019 IF( ierr.GT.0 )
1020 $ info = ierr
1021*
1022 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1023 IF( scaloc.NE.one ) THEN
1024 DO 170 k = 1, n
1025 CALL dscal( m, scaloc, c( 1, k ), 1 )
1026 CALL dscal( m, scaloc, f( 1, k ), 1 )
1027 170 CONTINUE
1028 scale = scale*scaloc
1029 END IF
1030*
1031* Unpack solution vector(s)
1032*
1033 k = 1
1034 ii = mb*nb + 1
1035 DO 180 jj = 0, nb - 1
1036 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1037 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1038 k = k + mb
1039 ii = ii + mb
1040 180 CONTINUE
1041*
1042* Substitute R(I, J) and L(I, J) into remaining
1043* equation.
1044*
1045 IF( j.GT.p+2 ) THEN
1046 CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1047 $ c( is, js ), ldc, b( 1, js ), ldb, one,
1048 $ f( is, 1 ), ldf )
1049 CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1050 $ f( is, js ), ldf, e( 1, js ), lde, one,
1051 $ f( is, 1 ), ldf )
1052 END IF
1053 IF( i.LT.p ) THEN
1054 CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1055 $ a( is, ie+1 ), lda, c( is, js ), ldc,
1056 $ one, c( ie+1, js ), ldc )
1057 CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1058 $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1059 $ one, c( ie+1, js ), ldc )
1060 END IF
1061*
1062 END IF
1063*
1064 190 CONTINUE
1065 200 CONTINUE
1066*
1067 END IF
1068 RETURN
1069*
1070* End of DTGSY2
1071*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dgesc2
subroutine dgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition dgesc2.f:114
dgetc2
subroutine dgetc2(n, a, lda, ipiv, jpiv, info)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition dgetc2.f:111
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dlatdf
subroutine dlatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition dlatdf.f:171
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: