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