1: *> \brief \b DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
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 communication 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 December 2016
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: * =====================================================================
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: *
277: * -- LAPACK auxiliary routine (version 3.7.0) --
278: * -- LAPACK is a software package provided by Univ. of Tennessee, --
279: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280: * December 2016
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 = -6
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
749: * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
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
759: IE = IWORK ( I+1 ) - 1
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: *
770: * Build a 2-by-2 system Z**T * x = RHS
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: *
782: * Solve Z**T * x = RHS
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: *
824: * Build a 4-by-4 system Z**T * x = RHS
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: *
853: * Solve Z**T * x = RHS
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: *
896: * Build a 4-by-4 system Z**T * x = RHS
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: *
925: * Solve Z**T * x = RHS
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: *
967: * Build an 8-by-8 system Z**T * x = RHS
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: *
1019: * Solve Z**T * x = RHS
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