Annotation of rpl/lapack/lapack/dtgsy2.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b DTGSY2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTGSY2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsy2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsy2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsy2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
! 22: * LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
! 23: * IWORK, PQ, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER TRANS
! 27: * INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
! 28: * $ PQ
! 29: * DOUBLE PRECISION RDSCAL, RDSUM, SCALE
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * INTEGER IWORK( * )
! 33: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
! 34: * $ D( LDD, * ), E( LDE, * ), F( LDF, * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> DTGSY2 solves the generalized Sylvester equation:
! 44: *>
! 45: *> A * R - L * B = scale * C (1)
! 46: *> D * R - L * E = scale * F,
! 47: *>
! 48: *> using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
! 49: *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
! 50: *> N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
! 51: *> must be in generalized Schur canonical form, i.e. A, B are upper
! 52: *> quasi triangular and D, E are upper triangular. The solution (R, L)
! 53: *> overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
! 54: *> chosen to avoid overflow.
! 55: *>
! 56: *> In matrix notation solving equation (1) corresponds to solve
! 57: *> Z*x = scale*b, where Z is defined as
! 58: *>
! 59: *> Z = [ kron(In, A) -kron(B**T, Im) ] (2)
! 60: *> [ kron(In, D) -kron(E**T, Im) ],
! 61: *>
! 62: *> Ik is the identity matrix of size k and X**T is the transpose of X.
! 63: *> kron(X, Y) is the Kronecker product between the matrices X and Y.
! 64: *> In the process of solving (1), we solve a number of such systems
! 65: *> where Dim(In), Dim(In) = 1 or 2.
! 66: *>
! 67: *> If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
! 68: *> which is equivalent to solve for R and L in
! 69: *>
! 70: *> A**T * R + D**T * L = scale * C (3)
! 71: *> R * B**T + L * E**T = scale * -F
! 72: *>
! 73: *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
! 74: *> sigma_min(Z) using reverse communicaton with DLACON.
! 75: *>
! 76: *> DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
! 77: *> of an upper bound on the separation between to matrix pairs. Then
! 78: *> the input (A, D), (B, E) are sub-pencils of the matrix pair in
! 79: *> DTGSYL. See DTGSYL for details.
! 80: *> \endverbatim
! 81: *
! 82: * Arguments:
! 83: * ==========
! 84: *
! 85: *> \param[in] TRANS
! 86: *> \verbatim
! 87: *> TRANS is CHARACTER*1
! 88: *> = 'N', solve the generalized Sylvester equation (1).
! 89: *> = 'T': solve the 'transposed' system (3).
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] IJOB
! 93: *> \verbatim
! 94: *> IJOB is INTEGER
! 95: *> Specifies what kind of functionality to be performed.
! 96: *> = 0: solve (1) only.
! 97: *> = 1: A contribution from this subsystem to a Frobenius
! 98: *> norm-based estimate of the separation between two matrix
! 99: *> pairs is computed. (look ahead strategy is used).
! 100: *> = 2: A contribution from this subsystem to a Frobenius
! 101: *> norm-based estimate of the separation between two matrix
! 102: *> pairs is computed. (DGECON on sub-systems is used.)
! 103: *> Not referenced if TRANS = 'T'.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in] M
! 107: *> \verbatim
! 108: *> M is INTEGER
! 109: *> On entry, M specifies the order of A and D, and the row
! 110: *> dimension of C, F, R and L.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] N
! 114: *> \verbatim
! 115: *> N is INTEGER
! 116: *> On entry, N specifies the order of B and E, and the column
! 117: *> dimension of C, F, R and L.
! 118: *> \endverbatim
! 119: *>
! 120: *> \param[in] A
! 121: *> \verbatim
! 122: *> A is DOUBLE PRECISION array, dimension (LDA, M)
! 123: *> On entry, A contains an upper quasi triangular matrix.
! 124: *> \endverbatim
! 125: *>
! 126: *> \param[in] LDA
! 127: *> \verbatim
! 128: *> LDA is INTEGER
! 129: *> The leading dimension of the matrix A. LDA >= max(1, M).
! 130: *> \endverbatim
! 131: *>
! 132: *> \param[in] B
! 133: *> \verbatim
! 134: *> B is DOUBLE PRECISION array, dimension (LDB, N)
! 135: *> On entry, B contains an upper quasi triangular matrix.
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[in] LDB
! 139: *> \verbatim
! 140: *> LDB is INTEGER
! 141: *> The leading dimension of the matrix B. LDB >= max(1, N).
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[in,out] C
! 145: *> \verbatim
! 146: *> C is DOUBLE PRECISION array, dimension (LDC, N)
! 147: *> On entry, C contains the right-hand-side of the first matrix
! 148: *> equation in (1).
! 149: *> On exit, if IJOB = 0, C has been overwritten by the
! 150: *> solution R.
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[in] LDC
! 154: *> \verbatim
! 155: *> LDC is INTEGER
! 156: *> The leading dimension of the matrix C. LDC >= max(1, M).
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[in] D
! 160: *> \verbatim
! 161: *> D is DOUBLE PRECISION array, dimension (LDD, M)
! 162: *> On entry, D contains an upper triangular matrix.
! 163: *> \endverbatim
! 164: *>
! 165: *> \param[in] LDD
! 166: *> \verbatim
! 167: *> LDD is INTEGER
! 168: *> The leading dimension of the matrix D. LDD >= max(1, M).
! 169: *> \endverbatim
! 170: *>
! 171: *> \param[in] E
! 172: *> \verbatim
! 173: *> E is DOUBLE PRECISION array, dimension (LDE, N)
! 174: *> On entry, E contains an upper triangular matrix.
! 175: *> \endverbatim
! 176: *>
! 177: *> \param[in] LDE
! 178: *> \verbatim
! 179: *> LDE is INTEGER
! 180: *> The leading dimension of the matrix E. LDE >= max(1, N).
! 181: *> \endverbatim
! 182: *>
! 183: *> \param[in,out] F
! 184: *> \verbatim
! 185: *> F is DOUBLE PRECISION array, dimension (LDF, N)
! 186: *> On entry, F contains the right-hand-side of the second matrix
! 187: *> equation in (1).
! 188: *> On exit, if IJOB = 0, F has been overwritten by the
! 189: *> solution L.
! 190: *> \endverbatim
! 191: *>
! 192: *> \param[in] LDF
! 193: *> \verbatim
! 194: *> LDF is INTEGER
! 195: *> The leading dimension of the matrix F. LDF >= max(1, M).
! 196: *> \endverbatim
! 197: *>
! 198: *> \param[out] SCALE
! 199: *> \verbatim
! 200: *> SCALE is DOUBLE PRECISION
! 201: *> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
! 202: *> R and L (C and F on entry) will hold the solutions to a
! 203: *> slightly perturbed system but the input matrices A, B, D and
! 204: *> E have not been changed. If SCALE = 0, R and L will hold the
! 205: *> solutions to the homogeneous system with C = F = 0. Normally,
! 206: *> SCALE = 1.
! 207: *> \endverbatim
! 208: *>
! 209: *> \param[in,out] RDSUM
! 210: *> \verbatim
! 211: *> RDSUM is DOUBLE PRECISION
! 212: *> On entry, the sum of squares of computed contributions to
! 213: *> the Dif-estimate under computation by DTGSYL, where the
! 214: *> scaling factor RDSCAL (see below) has been factored out.
! 215: *> On exit, the corresponding sum of squares updated with the
! 216: *> contributions from the current sub-system.
! 217: *> If TRANS = 'T' RDSUM is not touched.
! 218: *> NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
! 219: *> \endverbatim
! 220: *>
! 221: *> \param[in,out] RDSCAL
! 222: *> \verbatim
! 223: *> RDSCAL is DOUBLE PRECISION
! 224: *> On entry, scaling factor used to prevent overflow in RDSUM.
! 225: *> On exit, RDSCAL is updated w.r.t. the current contributions
! 226: *> in RDSUM.
! 227: *> If TRANS = 'T', RDSCAL is not touched.
! 228: *> NOTE: RDSCAL only makes sense when DTGSY2 is called by
! 229: *> DTGSYL.
! 230: *> \endverbatim
! 231: *>
! 232: *> \param[out] IWORK
! 233: *> \verbatim
! 234: *> IWORK is INTEGER array, dimension (M+N+2)
! 235: *> \endverbatim
! 236: *>
! 237: *> \param[out] PQ
! 238: *> \verbatim
! 239: *> PQ is INTEGER
! 240: *> On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
! 241: *> 8-by-8) solved by this routine.
! 242: *> \endverbatim
! 243: *>
! 244: *> \param[out] INFO
! 245: *> \verbatim
! 246: *> INFO is INTEGER
! 247: *> On exit, if INFO is set to
! 248: *> =0: Successful exit
! 249: *> <0: If INFO = -i, the i-th argument had an illegal value.
! 250: *> >0: The matrix pairs (A, D) and (B, E) have common or very
! 251: *> close eigenvalues.
! 252: *> \endverbatim
! 253: *
! 254: * Authors:
! 255: * ========
! 256: *
! 257: *> \author Univ. of Tennessee
! 258: *> \author Univ. of California Berkeley
! 259: *> \author Univ. of Colorado Denver
! 260: *> \author NAG Ltd.
! 261: *
! 262: *> \date November 2011
! 263: *
! 264: *> \ingroup doubleSYauxiliary
! 265: *
! 266: *> \par Contributors:
! 267: * ==================
! 268: *>
! 269: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
! 270: *> Umea University, S-901 87 Umea, Sweden.
! 271: *
! 272: * =====================================================================
1.1 bertrand 273: SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
274: $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
275: $ IWORK, PQ, INFO )
276: *
1.10 ! bertrand 277: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 278: * -- LAPACK is a software package provided by Univ. of Tennessee, --
279: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 280: * November 2011
1.1 bertrand 281: *
282: * .. Scalar Arguments ..
283: CHARACTER TRANS
284: INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
285: $ PQ
286: DOUBLE PRECISION RDSCAL, RDSUM, SCALE
287: * ..
288: * .. Array Arguments ..
289: INTEGER IWORK( * )
290: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
291: $ D( LDD, * ), E( LDE, * ), F( LDF, * )
292: * ..
293: *
294: * =====================================================================
295: * Replaced various illegal calls to DCOPY by calls to DLASET.
296: * Sven Hammarling, 27/5/02.
297: *
298: * .. Parameters ..
299: INTEGER LDZ
300: PARAMETER ( LDZ = 8 )
301: DOUBLE PRECISION ZERO, ONE
302: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
303: * ..
304: * .. Local Scalars ..
305: LOGICAL NOTRAN
306: INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
307: $ K, MB, NB, P, Q, ZDIM
308: DOUBLE PRECISION ALPHA, SCALOC
309: * ..
310: * .. Local Arrays ..
311: INTEGER IPIV( LDZ ), JPIV( LDZ )
312: DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
313: * ..
314: * .. External Functions ..
315: LOGICAL LSAME
316: EXTERNAL LSAME
317: * ..
318: * .. External Subroutines ..
319: EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2,
320: $ DGETC2, DLASET, DLATDF, DSCAL, XERBLA
321: * ..
322: * .. Intrinsic Functions ..
323: INTRINSIC MAX
324: * ..
325: * .. Executable Statements ..
326: *
327: * Decode and test input parameters
328: *
329: INFO = 0
330: IERR = 0
331: NOTRAN = LSAME( TRANS, 'N' )
332: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
333: INFO = -1
334: ELSE IF( NOTRAN ) THEN
335: IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
336: INFO = -2
337: END IF
338: END IF
339: IF( INFO.EQ.0 ) THEN
340: IF( M.LE.0 ) THEN
341: INFO = -3
342: ELSE IF( N.LE.0 ) THEN
343: INFO = -4
344: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
345: INFO = -5
346: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
347: INFO = -8
348: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
349: INFO = -10
350: ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
351: INFO = -12
352: ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
353: INFO = -14
354: ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
355: INFO = -16
356: END IF
357: END IF
358: IF( INFO.NE.0 ) THEN
359: CALL XERBLA( 'DTGSY2', -INFO )
360: RETURN
361: END IF
362: *
363: * Determine block structure of A
364: *
365: PQ = 0
366: P = 0
367: I = 1
368: 10 CONTINUE
369: IF( I.GT.M )
370: $ GO TO 20
371: P = P + 1
372: IWORK( P ) = I
373: IF( I.EQ.M )
374: $ GO TO 20
375: IF( A( I+1, I ).NE.ZERO ) THEN
376: I = I + 2
377: ELSE
378: I = I + 1
379: END IF
380: GO TO 10
381: 20 CONTINUE
382: IWORK( P+1 ) = M + 1
383: *
384: * Determine block structure of B
385: *
386: Q = P + 1
387: J = 1
388: 30 CONTINUE
389: IF( J.GT.N )
390: $ GO TO 40
391: Q = Q + 1
392: IWORK( Q ) = J
393: IF( J.EQ.N )
394: $ GO TO 40
395: IF( B( J+1, J ).NE.ZERO ) THEN
396: J = J + 2
397: ELSE
398: J = J + 1
399: END IF
400: GO TO 30
401: 40 CONTINUE
402: IWORK( Q+1 ) = N + 1
403: PQ = P*( Q-P-1 )
404: *
405: IF( NOTRAN ) THEN
406: *
407: * Solve (I, J) - subsystem
408: * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
409: * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
410: * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
411: *
412: SCALE = ONE
413: SCALOC = ONE
414: DO 120 J = P + 2, Q
415: JS = IWORK( J )
416: JSP1 = JS + 1
417: JE = IWORK( J+1 ) - 1
418: NB = JE - JS + 1
419: DO 110 I = P, 1, -1
420: *
421: IS = IWORK( I )
422: ISP1 = IS + 1
423: IE = IWORK( I+1 ) - 1
424: MB = IE - IS + 1
425: ZDIM = MB*NB*2
426: *
427: IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
428: *
429: * Build a 2-by-2 system Z * x = RHS
430: *
431: Z( 1, 1 ) = A( IS, IS )
432: Z( 2, 1 ) = D( IS, IS )
433: Z( 1, 2 ) = -B( JS, JS )
434: Z( 2, 2 ) = -E( JS, JS )
435: *
436: * Set up right hand side(s)
437: *
438: RHS( 1 ) = C( IS, JS )
439: RHS( 2 ) = F( IS, JS )
440: *
441: * Solve Z * x = RHS
442: *
443: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
444: IF( IERR.GT.0 )
445: $ INFO = IERR
446: *
447: IF( IJOB.EQ.0 ) THEN
448: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
449: $ SCALOC )
450: IF( SCALOC.NE.ONE ) THEN
451: DO 50 K = 1, N
452: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
453: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
454: 50 CONTINUE
455: SCALE = SCALE*SCALOC
456: END IF
457: ELSE
458: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
459: $ RDSCAL, IPIV, JPIV )
460: END IF
461: *
462: * Unpack solution vector(s)
463: *
464: C( IS, JS ) = RHS( 1 )
465: F( IS, JS ) = RHS( 2 )
466: *
467: * Substitute R(I, J) and L(I, J) into remaining
468: * equation.
469: *
470: IF( I.GT.1 ) THEN
471: ALPHA = -RHS( 1 )
472: CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
473: $ 1 )
474: CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
475: $ 1 )
476: END IF
477: IF( J.LT.Q ) THEN
478: CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
479: $ C( IS, JE+1 ), LDC )
480: CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE,
481: $ F( IS, JE+1 ), LDF )
482: END IF
483: *
484: ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
485: *
486: * Build a 4-by-4 system Z * x = RHS
487: *
488: Z( 1, 1 ) = A( IS, IS )
489: Z( 2, 1 ) = ZERO
490: Z( 3, 1 ) = D( IS, IS )
491: Z( 4, 1 ) = ZERO
492: *
493: Z( 1, 2 ) = ZERO
494: Z( 2, 2 ) = A( IS, IS )
495: Z( 3, 2 ) = ZERO
496: Z( 4, 2 ) = D( IS, IS )
497: *
498: Z( 1, 3 ) = -B( JS, JS )
499: Z( 2, 3 ) = -B( JS, JSP1 )
500: Z( 3, 3 ) = -E( JS, JS )
501: Z( 4, 3 ) = -E( JS, JSP1 )
502: *
503: Z( 1, 4 ) = -B( JSP1, JS )
504: Z( 2, 4 ) = -B( JSP1, JSP1 )
505: Z( 3, 4 ) = ZERO
506: Z( 4, 4 ) = -E( JSP1, JSP1 )
507: *
508: * Set up right hand side(s)
509: *
510: RHS( 1 ) = C( IS, JS )
511: RHS( 2 ) = C( IS, JSP1 )
512: RHS( 3 ) = F( IS, JS )
513: RHS( 4 ) = F( IS, JSP1 )
514: *
515: * Solve Z * x = RHS
516: *
517: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
518: IF( IERR.GT.0 )
519: $ INFO = IERR
520: *
521: IF( IJOB.EQ.0 ) THEN
522: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
523: $ SCALOC )
524: IF( SCALOC.NE.ONE ) THEN
525: DO 60 K = 1, N
526: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
527: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
528: 60 CONTINUE
529: SCALE = SCALE*SCALOC
530: END IF
531: ELSE
532: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
533: $ RDSCAL, IPIV, JPIV )
534: END IF
535: *
536: * Unpack solution vector(s)
537: *
538: C( IS, JS ) = RHS( 1 )
539: C( IS, JSP1 ) = RHS( 2 )
540: F( IS, JS ) = RHS( 3 )
541: F( IS, JSP1 ) = RHS( 4 )
542: *
543: * Substitute R(I, J) and L(I, J) into remaining
544: * equation.
545: *
546: IF( I.GT.1 ) THEN
547: CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
548: $ 1, C( 1, JS ), LDC )
549: CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
550: $ 1, F( 1, JS ), LDF )
551: END IF
552: IF( J.LT.Q ) THEN
553: CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
554: $ C( IS, JE+1 ), LDC )
555: CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
556: $ F( IS, JE+1 ), LDF )
557: CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
558: $ C( IS, JE+1 ), LDC )
559: CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
560: $ F( IS, JE+1 ), LDF )
561: END IF
562: *
563: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
564: *
565: * Build a 4-by-4 system Z * x = RHS
566: *
567: Z( 1, 1 ) = A( IS, IS )
568: Z( 2, 1 ) = A( ISP1, IS )
569: Z( 3, 1 ) = D( IS, IS )
570: Z( 4, 1 ) = ZERO
571: *
572: Z( 1, 2 ) = A( IS, ISP1 )
573: Z( 2, 2 ) = A( ISP1, ISP1 )
574: Z( 3, 2 ) = D( IS, ISP1 )
575: Z( 4, 2 ) = D( ISP1, ISP1 )
576: *
577: Z( 1, 3 ) = -B( JS, JS )
578: Z( 2, 3 ) = ZERO
579: Z( 3, 3 ) = -E( JS, JS )
580: Z( 4, 3 ) = ZERO
581: *
582: Z( 1, 4 ) = ZERO
583: Z( 2, 4 ) = -B( JS, JS )
584: Z( 3, 4 ) = ZERO
585: Z( 4, 4 ) = -E( JS, JS )
586: *
587: * Set up right hand side(s)
588: *
589: RHS( 1 ) = C( IS, JS )
590: RHS( 2 ) = C( ISP1, JS )
591: RHS( 3 ) = F( IS, JS )
592: RHS( 4 ) = F( ISP1, JS )
593: *
594: * Solve Z * x = RHS
595: *
596: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
597: IF( IERR.GT.0 )
598: $ INFO = IERR
599: IF( IJOB.EQ.0 ) THEN
600: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
601: $ SCALOC )
602: IF( SCALOC.NE.ONE ) THEN
603: DO 70 K = 1, N
604: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
605: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
606: 70 CONTINUE
607: SCALE = SCALE*SCALOC
608: END IF
609: ELSE
610: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
611: $ RDSCAL, IPIV, JPIV )
612: END IF
613: *
614: * Unpack solution vector(s)
615: *
616: C( IS, JS ) = RHS( 1 )
617: C( ISP1, JS ) = RHS( 2 )
618: F( IS, JS ) = RHS( 3 )
619: F( ISP1, JS ) = RHS( 4 )
620: *
621: * Substitute R(I, J) and L(I, J) into remaining
622: * equation.
623: *
624: IF( I.GT.1 ) THEN
625: CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
626: $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
627: CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD,
628: $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
629: END IF
630: IF( J.LT.Q ) THEN
631: CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
632: $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
633: CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
634: $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
635: END IF
636: *
637: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
638: *
639: * Build an 8-by-8 system Z * x = RHS
640: *
641: CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
642: *
643: Z( 1, 1 ) = A( IS, IS )
644: Z( 2, 1 ) = A( ISP1, IS )
645: Z( 5, 1 ) = D( IS, IS )
646: *
647: Z( 1, 2 ) = A( IS, ISP1 )
648: Z( 2, 2 ) = A( ISP1, ISP1 )
649: Z( 5, 2 ) = D( IS, ISP1 )
650: Z( 6, 2 ) = D( ISP1, ISP1 )
651: *
652: Z( 3, 3 ) = A( IS, IS )
653: Z( 4, 3 ) = A( ISP1, IS )
654: Z( 7, 3 ) = D( IS, IS )
655: *
656: Z( 3, 4 ) = A( IS, ISP1 )
657: Z( 4, 4 ) = A( ISP1, ISP1 )
658: Z( 7, 4 ) = D( IS, ISP1 )
659: Z( 8, 4 ) = D( ISP1, ISP1 )
660: *
661: Z( 1, 5 ) = -B( JS, JS )
662: Z( 3, 5 ) = -B( JS, JSP1 )
663: Z( 5, 5 ) = -E( JS, JS )
664: Z( 7, 5 ) = -E( JS, JSP1 )
665: *
666: Z( 2, 6 ) = -B( JS, JS )
667: Z( 4, 6 ) = -B( JS, JSP1 )
668: Z( 6, 6 ) = -E( JS, JS )
669: Z( 8, 6 ) = -E( JS, JSP1 )
670: *
671: Z( 1, 7 ) = -B( JSP1, JS )
672: Z( 3, 7 ) = -B( JSP1, JSP1 )
673: Z( 7, 7 ) = -E( JSP1, JSP1 )
674: *
675: Z( 2, 8 ) = -B( JSP1, JS )
676: Z( 4, 8 ) = -B( JSP1, JSP1 )
677: Z( 8, 8 ) = -E( JSP1, JSP1 )
678: *
679: * Set up right hand side(s)
680: *
681: K = 1
682: II = MB*NB + 1
683: DO 80 JJ = 0, NB - 1
684: CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
685: CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
686: K = K + MB
687: II = II + MB
688: 80 CONTINUE
689: *
690: * Solve Z * x = RHS
691: *
692: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
693: IF( IERR.GT.0 )
694: $ INFO = IERR
695: IF( IJOB.EQ.0 ) THEN
696: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
697: $ SCALOC )
698: IF( SCALOC.NE.ONE ) THEN
699: DO 90 K = 1, N
700: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
701: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
702: 90 CONTINUE
703: SCALE = SCALE*SCALOC
704: END IF
705: ELSE
706: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
707: $ RDSCAL, IPIV, JPIV )
708: END IF
709: *
710: * Unpack solution vector(s)
711: *
712: K = 1
713: II = MB*NB + 1
714: DO 100 JJ = 0, NB - 1
715: CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
716: CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
717: K = K + MB
718: II = II + MB
719: 100 CONTINUE
720: *
721: * Substitute R(I, J) and L(I, J) into remaining
722: * equation.
723: *
724: IF( I.GT.1 ) THEN
725: CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
726: $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
727: $ C( 1, JS ), LDC )
728: CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
729: $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE,
730: $ F( 1, JS ), LDF )
731: END IF
732: IF( J.LT.Q ) THEN
733: K = MB*NB + 1
734: CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
735: $ MB, B( JS, JE+1 ), LDB, ONE,
736: $ C( IS, JE+1 ), LDC )
737: CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
738: $ MB, E( JS, JE+1 ), LDE, ONE,
739: $ F( IS, JE+1 ), LDF )
740: END IF
741: *
742: END IF
743: *
744: 110 CONTINUE
745: 120 CONTINUE
746: ELSE
747: *
748: * Solve (I, J) - subsystem
1.9 bertrand 749: * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
1.1 bertrand 750: * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
751: * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
752: *
753: SCALE = ONE
754: SCALOC = ONE
755: DO 200 I = 1, P
756: *
757: IS = IWORK( I )
758: ISP1 = IS + 1
1.5 bertrand 759: IE = IWORK ( I+1 ) - 1
1.1 bertrand 760: MB = IE - IS + 1
761: DO 190 J = Q, P + 2, -1
762: *
763: JS = IWORK( J )
764: JSP1 = JS + 1
765: JE = IWORK( J+1 ) - 1
766: NB = JE - JS + 1
767: ZDIM = MB*NB*2
768: IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
769: *
1.9 bertrand 770: * Build a 2-by-2 system Z**T * x = RHS
1.1 bertrand 771: *
772: Z( 1, 1 ) = A( IS, IS )
773: Z( 2, 1 ) = -B( JS, JS )
774: Z( 1, 2 ) = D( IS, IS )
775: Z( 2, 2 ) = -E( JS, JS )
776: *
777: * Set up right hand side(s)
778: *
779: RHS( 1 ) = C( IS, JS )
780: RHS( 2 ) = F( IS, JS )
781: *
1.9 bertrand 782: * Solve Z**T * x = RHS
1.1 bertrand 783: *
784: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
785: IF( IERR.GT.0 )
786: $ INFO = IERR
787: *
788: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
789: IF( SCALOC.NE.ONE ) THEN
790: DO 130 K = 1, N
791: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
792: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
793: 130 CONTINUE
794: SCALE = SCALE*SCALOC
795: END IF
796: *
797: * Unpack solution vector(s)
798: *
799: C( IS, JS ) = RHS( 1 )
800: F( IS, JS ) = RHS( 2 )
801: *
802: * Substitute R(I, J) and L(I, J) into remaining
803: * equation.
804: *
805: IF( J.GT.P+2 ) THEN
806: ALPHA = RHS( 1 )
807: CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
808: $ LDF )
809: ALPHA = RHS( 2 )
810: CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
811: $ LDF )
812: END IF
813: IF( I.LT.P ) THEN
814: ALPHA = -RHS( 1 )
815: CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
816: $ C( IE+1, JS ), 1 )
817: ALPHA = -RHS( 2 )
818: CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
819: $ C( IE+1, JS ), 1 )
820: END IF
821: *
822: ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
823: *
1.9 bertrand 824: * Build a 4-by-4 system Z**T * x = RHS
1.1 bertrand 825: *
826: Z( 1, 1 ) = A( IS, IS )
827: Z( 2, 1 ) = ZERO
828: Z( 3, 1 ) = -B( JS, JS )
829: Z( 4, 1 ) = -B( JSP1, JS )
830: *
831: Z( 1, 2 ) = ZERO
832: Z( 2, 2 ) = A( IS, IS )
833: Z( 3, 2 ) = -B( JS, JSP1 )
834: Z( 4, 2 ) = -B( JSP1, JSP1 )
835: *
836: Z( 1, 3 ) = D( IS, IS )
837: Z( 2, 3 ) = ZERO
838: Z( 3, 3 ) = -E( JS, JS )
839: Z( 4, 3 ) = ZERO
840: *
841: Z( 1, 4 ) = ZERO
842: Z( 2, 4 ) = D( IS, IS )
843: Z( 3, 4 ) = -E( JS, JSP1 )
844: Z( 4, 4 ) = -E( JSP1, JSP1 )
845: *
846: * Set up right hand side(s)
847: *
848: RHS( 1 ) = C( IS, JS )
849: RHS( 2 ) = C( IS, JSP1 )
850: RHS( 3 ) = F( IS, JS )
851: RHS( 4 ) = F( IS, JSP1 )
852: *
1.9 bertrand 853: * Solve Z**T * x = RHS
1.1 bertrand 854: *
855: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
856: IF( IERR.GT.0 )
857: $ INFO = IERR
858: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
859: IF( SCALOC.NE.ONE ) THEN
860: DO 140 K = 1, N
861: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
862: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
863: 140 CONTINUE
864: SCALE = SCALE*SCALOC
865: END IF
866: *
867: * Unpack solution vector(s)
868: *
869: C( IS, JS ) = RHS( 1 )
870: C( IS, JSP1 ) = RHS( 2 )
871: F( IS, JS ) = RHS( 3 )
872: F( IS, JSP1 ) = RHS( 4 )
873: *
874: * Substitute R(I, J) and L(I, J) into remaining
875: * equation.
876: *
877: IF( J.GT.P+2 ) THEN
878: CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
879: $ F( IS, 1 ), LDF )
880: CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
881: $ F( IS, 1 ), LDF )
882: CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
883: $ F( IS, 1 ), LDF )
884: CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
885: $ F( IS, 1 ), LDF )
886: END IF
887: IF( I.LT.P ) THEN
888: CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
889: $ RHS( 1 ), 1, C( IE+1, JS ), LDC )
890: CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
891: $ RHS( 3 ), 1, C( IE+1, JS ), LDC )
892: END IF
893: *
894: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
895: *
1.9 bertrand 896: * Build a 4-by-4 system Z**T * x = RHS
1.1 bertrand 897: *
898: Z( 1, 1 ) = A( IS, IS )
899: Z( 2, 1 ) = A( IS, ISP1 )
900: Z( 3, 1 ) = -B( JS, JS )
901: Z( 4, 1 ) = ZERO
902: *
903: Z( 1, 2 ) = A( ISP1, IS )
904: Z( 2, 2 ) = A( ISP1, ISP1 )
905: Z( 3, 2 ) = ZERO
906: Z( 4, 2 ) = -B( JS, JS )
907: *
908: Z( 1, 3 ) = D( IS, IS )
909: Z( 2, 3 ) = D( IS, ISP1 )
910: Z( 3, 3 ) = -E( JS, JS )
911: Z( 4, 3 ) = ZERO
912: *
913: Z( 1, 4 ) = ZERO
914: Z( 2, 4 ) = D( ISP1, ISP1 )
915: Z( 3, 4 ) = ZERO
916: Z( 4, 4 ) = -E( JS, JS )
917: *
918: * Set up right hand side(s)
919: *
920: RHS( 1 ) = C( IS, JS )
921: RHS( 2 ) = C( ISP1, JS )
922: RHS( 3 ) = F( IS, JS )
923: RHS( 4 ) = F( ISP1, JS )
924: *
1.9 bertrand 925: * Solve Z**T * x = RHS
1.1 bertrand 926: *
927: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
928: IF( IERR.GT.0 )
929: $ INFO = IERR
930: *
931: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
932: IF( SCALOC.NE.ONE ) THEN
933: DO 150 K = 1, N
934: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
935: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
936: 150 CONTINUE
937: SCALE = SCALE*SCALOC
938: END IF
939: *
940: * Unpack solution vector(s)
941: *
942: C( IS, JS ) = RHS( 1 )
943: C( ISP1, JS ) = RHS( 2 )
944: F( IS, JS ) = RHS( 3 )
945: F( ISP1, JS ) = RHS( 4 )
946: *
947: * Substitute R(I, J) and L(I, J) into remaining
948: * equation.
949: *
950: IF( J.GT.P+2 ) THEN
951: CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
952: $ 1, F( IS, 1 ), LDF )
953: CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
954: $ 1, F( IS, 1 ), LDF )
955: END IF
956: IF( I.LT.P ) THEN
957: CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
958: $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
959: $ 1 )
960: CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
961: $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
962: $ 1 )
963: END IF
964: *
965: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
966: *
1.9 bertrand 967: * Build an 8-by-8 system Z**T * x = RHS
1.1 bertrand 968: *
969: CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
970: *
971: Z( 1, 1 ) = A( IS, IS )
972: Z( 2, 1 ) = A( IS, ISP1 )
973: Z( 5, 1 ) = -B( JS, JS )
974: Z( 7, 1 ) = -B( JSP1, JS )
975: *
976: Z( 1, 2 ) = A( ISP1, IS )
977: Z( 2, 2 ) = A( ISP1, ISP1 )
978: Z( 6, 2 ) = -B( JS, JS )
979: Z( 8, 2 ) = -B( JSP1, JS )
980: *
981: Z( 3, 3 ) = A( IS, IS )
982: Z( 4, 3 ) = A( IS, ISP1 )
983: Z( 5, 3 ) = -B( JS, JSP1 )
984: Z( 7, 3 ) = -B( JSP1, JSP1 )
985: *
986: Z( 3, 4 ) = A( ISP1, IS )
987: Z( 4, 4 ) = A( ISP1, ISP1 )
988: Z( 6, 4 ) = -B( JS, JSP1 )
989: Z( 8, 4 ) = -B( JSP1, JSP1 )
990: *
991: Z( 1, 5 ) = D( IS, IS )
992: Z( 2, 5 ) = D( IS, ISP1 )
993: Z( 5, 5 ) = -E( JS, JS )
994: *
995: Z( 2, 6 ) = D( ISP1, ISP1 )
996: Z( 6, 6 ) = -E( JS, JS )
997: *
998: Z( 3, 7 ) = D( IS, IS )
999: Z( 4, 7 ) = D( IS, ISP1 )
1000: Z( 5, 7 ) = -E( JS, JSP1 )
1001: Z( 7, 7 ) = -E( JSP1, JSP1 )
1002: *
1003: Z( 4, 8 ) = D( ISP1, ISP1 )
1004: Z( 6, 8 ) = -E( JS, JSP1 )
1005: Z( 8, 8 ) = -E( JSP1, JSP1 )
1006: *
1007: * Set up right hand side(s)
1008: *
1009: K = 1
1010: II = MB*NB + 1
1011: DO 160 JJ = 0, NB - 1
1012: CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
1013: CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
1014: K = K + MB
1015: II = II + MB
1016: 160 CONTINUE
1017: *
1018: *
1.9 bertrand 1019: * Solve Z**T * x = RHS
1.1 bertrand 1020: *
1021: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
1022: IF( IERR.GT.0 )
1023: $ INFO = IERR
1024: *
1025: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
1026: IF( SCALOC.NE.ONE ) THEN
1027: DO 170 K = 1, N
1028: CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
1029: CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
1030: 170 CONTINUE
1031: SCALE = SCALE*SCALOC
1032: END IF
1033: *
1034: * Unpack solution vector(s)
1035: *
1036: K = 1
1037: II = MB*NB + 1
1038: DO 180 JJ = 0, NB - 1
1039: CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
1040: CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
1041: K = K + MB
1042: II = II + MB
1043: 180 CONTINUE
1044: *
1045: * Substitute R(I, J) and L(I, J) into remaining
1046: * equation.
1047: *
1048: IF( J.GT.P+2 ) THEN
1049: CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
1050: $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
1051: $ F( IS, 1 ), LDF )
1052: CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
1053: $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
1054: $ F( IS, 1 ), LDF )
1055: END IF
1056: IF( I.LT.P ) THEN
1057: CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
1058: $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
1059: $ ONE, C( IE+1, JS ), LDC )
1060: CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
1061: $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
1062: $ ONE, C( IE+1, JS ), LDC )
1063: END IF
1064: *
1065: END IF
1066: *
1067: 190 CONTINUE
1068: 200 CONTINUE
1069: *
1070: END IF
1071: RETURN
1072: *
1073: * End of DTGSY2
1074: *
1075: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgsy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack