![[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|
Return to dtgsy2.f CVS log | Up to [local] / rpl / lapack / lapack |
1.13 bertrand 1: *> \brief \b DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
1.10 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.18 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.18 bertrand 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">
1.10 bertrand 15: *> [TXT]</a>
1.18 bertrand 16: *> \endhtmlonly
1.10 bertrand 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 )
1.18 bertrand 24: *
1.10 bertrand 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: * ..
1.18 bertrand 36: *
1.10 bertrand 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)] =
1.21 bertrand 74: *> sigma_min(Z) using reverse communication with DLACON.
1.10 bertrand 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
1.21 bertrand 88: *> = 'N': solve the generalized Sylvester equation (1).
1.10 bertrand 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: *
1.18 bertrand 257: *> \author Univ. of Tennessee
258: *> \author Univ. of California Berkeley
259: *> \author Univ. of Colorado Denver
260: *> \author NAG Ltd.
1.10 bertrand 261: *
262: *> \ingroup doubleSYauxiliary
263: *
264: *> \par Contributors:
265: * ==================
266: *>
267: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
268: *> Umea University, S-901 87 Umea, Sweden.
269: *
270: * =====================================================================
1.1 bertrand 271: SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
272: $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
273: $ IWORK, PQ, INFO )
274: *
1.22 ! bertrand 275: * -- LAPACK auxiliary routine --
1.1 bertrand 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
1.16 bertrand 342: INFO = -6
1.1 bertrand 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
1.9 bertrand 746: * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
1.1 bertrand 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
1.5 bertrand 756: IE = IWORK ( I+1 ) - 1
1.1 bertrand 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: *
1.9 bertrand 767: * Build a 2-by-2 system Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 779: * Solve Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 821: * Build a 4-by-4 system Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 850: * Solve Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 893: * Build a 4-by-4 system Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 922: * Solve Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 964: * Build an 8-by-8 system Z**T * x = RHS
1.1 bertrand 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: *
1.9 bertrand 1016: * Solve Z**T * x = RHS
1.1 bertrand 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: *
1072: END