1: *> \brief \b ZGGBAL
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGGBAL + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggbal.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggbal.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggbal.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
22: * RSCALE, WORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOB
26: * INTEGER IHI, ILO, INFO, LDA, LDB, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * )
30: * COMPLEX*16 A( LDA, * ), B( LDB, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZGGBAL balances a pair of general complex matrices (A,B). This
40: *> involves, first, permuting A and B by similarity transformations to
41: *> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
42: *> elements on the diagonal; and second, applying a diagonal similarity
43: *> transformation to rows and columns ILO to IHI to make the rows
44: *> and columns as close in norm as possible. Both steps are optional.
45: *>
46: *> Balancing may reduce the 1-norm of the matrices, and improve the
47: *> accuracy of the computed eigenvalues and/or eigenvectors in the
48: *> generalized eigenvalue problem A*x = lambda*B*x.
49: *> \endverbatim
50: *
51: * Arguments:
52: * ==========
53: *
54: *> \param[in] JOB
55: *> \verbatim
56: *> JOB is CHARACTER*1
57: *> Specifies the operations to be performed on A and B:
58: *> = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
59: *> and RSCALE(I) = 1.0 for i=1,...,N;
60: *> = 'P': permute only;
61: *> = 'S': scale only;
62: *> = 'B': both permute and scale.
63: *> \endverbatim
64: *>
65: *> \param[in] N
66: *> \verbatim
67: *> N is INTEGER
68: *> The order of the matrices A and B. N >= 0.
69: *> \endverbatim
70: *>
71: *> \param[in,out] A
72: *> \verbatim
73: *> A is COMPLEX*16 array, dimension (LDA,N)
74: *> On entry, the input matrix A.
75: *> On exit, A is overwritten by the balanced matrix.
76: *> If JOB = 'N', A is not referenced.
77: *> \endverbatim
78: *>
79: *> \param[in] LDA
80: *> \verbatim
81: *> LDA is INTEGER
82: *> The leading dimension of the array A. LDA >= max(1,N).
83: *> \endverbatim
84: *>
85: *> \param[in,out] B
86: *> \verbatim
87: *> B is COMPLEX*16 array, dimension (LDB,N)
88: *> On entry, the input matrix B.
89: *> On exit, B is overwritten by the balanced matrix.
90: *> If JOB = 'N', B is not referenced.
91: *> \endverbatim
92: *>
93: *> \param[in] LDB
94: *> \verbatim
95: *> LDB is INTEGER
96: *> The leading dimension of the array B. LDB >= max(1,N).
97: *> \endverbatim
98: *>
99: *> \param[out] ILO
100: *> \verbatim
101: *> ILO is INTEGER
102: *> \endverbatim
103: *>
104: *> \param[out] IHI
105: *> \verbatim
106: *> IHI is INTEGER
107: *> ILO and IHI are set to integers such that on exit
108: *> A(i,j) = 0 and B(i,j) = 0 if i > j and
109: *> j = 1,...,ILO-1 or i = IHI+1,...,N.
110: *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
111: *> \endverbatim
112: *>
113: *> \param[out] LSCALE
114: *> \verbatim
115: *> LSCALE is DOUBLE PRECISION array, dimension (N)
116: *> Details of the permutations and scaling factors applied
117: *> to the left side of A and B. If P(j) is the index of the
118: *> row interchanged with row j, and D(j) is the scaling factor
119: *> applied to row j, then
120: *> LSCALE(j) = P(j) for J = 1,...,ILO-1
121: *> = D(j) for J = ILO,...,IHI
122: *> = P(j) for J = IHI+1,...,N.
123: *> The order in which the interchanges are made is N to IHI+1,
124: *> then 1 to ILO-1.
125: *> \endverbatim
126: *>
127: *> \param[out] RSCALE
128: *> \verbatim
129: *> RSCALE is DOUBLE PRECISION array, dimension (N)
130: *> Details of the permutations and scaling factors applied
131: *> to the right side of A and B. If P(j) is the index of the
132: *> column interchanged with column j, and D(j) is the scaling
133: *> factor applied to column j, then
134: *> RSCALE(j) = P(j) for J = 1,...,ILO-1
135: *> = D(j) for J = ILO,...,IHI
136: *> = P(j) for J = IHI+1,...,N.
137: *> The order in which the interchanges are made is N to IHI+1,
138: *> then 1 to ILO-1.
139: *> \endverbatim
140: *>
141: *> \param[out] WORK
142: *> \verbatim
143: *> WORK is DOUBLE PRECISION array, dimension (lwork)
144: *> lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
145: *> at least 1 when JOB = 'N' or 'P'.
146: *> \endverbatim
147: *>
148: *> \param[out] INFO
149: *> \verbatim
150: *> INFO is INTEGER
151: *> = 0: successful exit
152: *> < 0: if INFO = -i, the i-th argument had an illegal value.
153: *> \endverbatim
154: *
155: * Authors:
156: * ========
157: *
158: *> \author Univ. of Tennessee
159: *> \author Univ. of California Berkeley
160: *> \author Univ. of Colorado Denver
161: *> \author NAG Ltd.
162: *
163: *> \date June 2016
164: *
165: *> \ingroup complex16GBcomputational
166: *
167: *> \par Further Details:
168: * =====================
169: *>
170: *> \verbatim
171: *>
172: *> See R.C. WARD, Balancing the generalized eigenvalue problem,
173: *> SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
174: *> \endverbatim
175: *>
176: * =====================================================================
177: SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
178: $ RSCALE, WORK, INFO )
179: *
180: * -- LAPACK computational routine (version 3.6.1) --
181: * -- LAPACK is a software package provided by Univ. of Tennessee, --
182: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183: * June 2016
184: *
185: * .. Scalar Arguments ..
186: CHARACTER JOB
187: INTEGER IHI, ILO, INFO, LDA, LDB, N
188: * ..
189: * .. Array Arguments ..
190: DOUBLE PRECISION LSCALE( * ), RSCALE( * ), WORK( * )
191: COMPLEX*16 A( LDA, * ), B( LDB, * )
192: * ..
193: *
194: * =====================================================================
195: *
196: * .. Parameters ..
197: DOUBLE PRECISION ZERO, HALF, ONE
198: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
199: DOUBLE PRECISION THREE, SCLFAC
200: PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
201: COMPLEX*16 CZERO
202: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
203: * ..
204: * .. Local Scalars ..
205: INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
206: $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
207: $ M, NR, NRP2
208: DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
209: $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
210: $ SFMIN, SUM, T, TA, TB, TC
211: COMPLEX*16 CDUM
212: * ..
213: * .. External Functions ..
214: LOGICAL LSAME
215: INTEGER IZAMAX
216: DOUBLE PRECISION DDOT, DLAMCH
217: EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH
218: * ..
219: * .. External Subroutines ..
220: EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
221: * ..
222: * .. Intrinsic Functions ..
223: INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
224: * ..
225: * .. Statement Functions ..
226: DOUBLE PRECISION CABS1
227: * ..
228: * .. Statement Function definitions ..
229: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
230: * ..
231: * .. Executable Statements ..
232: *
233: * Test the input parameters
234: *
235: INFO = 0
236: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
237: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
238: INFO = -1
239: ELSE IF( N.LT.0 ) THEN
240: INFO = -2
241: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
242: INFO = -4
243: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
244: INFO = -6
245: END IF
246: IF( INFO.NE.0 ) THEN
247: CALL XERBLA( 'ZGGBAL', -INFO )
248: RETURN
249: END IF
250: *
251: * Quick return if possible
252: *
253: IF( N.EQ.0 ) THEN
254: ILO = 1
255: IHI = N
256: RETURN
257: END IF
258: *
259: IF( N.EQ.1 ) THEN
260: ILO = 1
261: IHI = N
262: LSCALE( 1 ) = ONE
263: RSCALE( 1 ) = ONE
264: RETURN
265: END IF
266: *
267: IF( LSAME( JOB, 'N' ) ) THEN
268: ILO = 1
269: IHI = N
270: DO 10 I = 1, N
271: LSCALE( I ) = ONE
272: RSCALE( I ) = ONE
273: 10 CONTINUE
274: RETURN
275: END IF
276: *
277: K = 1
278: L = N
279: IF( LSAME( JOB, 'S' ) )
280: $ GO TO 190
281: *
282: GO TO 30
283: *
284: * Permute the matrices A and B to isolate the eigenvalues.
285: *
286: * Find row with one nonzero in columns 1 through L
287: *
288: 20 CONTINUE
289: L = LM1
290: IF( L.NE.1 )
291: $ GO TO 30
292: *
293: RSCALE( 1 ) = 1
294: LSCALE( 1 ) = 1
295: GO TO 190
296: *
297: 30 CONTINUE
298: LM1 = L - 1
299: DO 80 I = L, 1, -1
300: DO 40 J = 1, LM1
301: JP1 = J + 1
302: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
303: $ GO TO 50
304: 40 CONTINUE
305: J = L
306: GO TO 70
307: *
308: 50 CONTINUE
309: DO 60 J = JP1, L
310: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
311: $ GO TO 80
312: 60 CONTINUE
313: J = JP1 - 1
314: *
315: 70 CONTINUE
316: M = L
317: IFLOW = 1
318: GO TO 160
319: 80 CONTINUE
320: GO TO 100
321: *
322: * Find column with one nonzero in rows K through N
323: *
324: 90 CONTINUE
325: K = K + 1
326: *
327: 100 CONTINUE
328: DO 150 J = K, L
329: DO 110 I = K, LM1
330: IP1 = I + 1
331: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
332: $ GO TO 120
333: 110 CONTINUE
334: I = L
335: GO TO 140
336: 120 CONTINUE
337: DO 130 I = IP1, L
338: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
339: $ GO TO 150
340: 130 CONTINUE
341: I = IP1 - 1
342: 140 CONTINUE
343: M = K
344: IFLOW = 2
345: GO TO 160
346: 150 CONTINUE
347: GO TO 190
348: *
349: * Permute rows M and I
350: *
351: 160 CONTINUE
352: LSCALE( M ) = I
353: IF( I.EQ.M )
354: $ GO TO 170
355: CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
356: CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
357: *
358: * Permute columns M and J
359: *
360: 170 CONTINUE
361: RSCALE( M ) = J
362: IF( J.EQ.M )
363: $ GO TO 180
364: CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
365: CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
366: *
367: 180 CONTINUE
368: GO TO ( 20, 90 )IFLOW
369: *
370: 190 CONTINUE
371: ILO = K
372: IHI = L
373: *
374: IF( LSAME( JOB, 'P' ) ) THEN
375: DO 195 I = ILO, IHI
376: LSCALE( I ) = ONE
377: RSCALE( I ) = ONE
378: 195 CONTINUE
379: RETURN
380: END IF
381: *
382: IF( ILO.EQ.IHI )
383: $ RETURN
384: *
385: * Balance the submatrix in rows ILO to IHI.
386: *
387: NR = IHI - ILO + 1
388: DO 200 I = ILO, IHI
389: RSCALE( I ) = ZERO
390: LSCALE( I ) = ZERO
391: *
392: WORK( I ) = ZERO
393: WORK( I+N ) = ZERO
394: WORK( I+2*N ) = ZERO
395: WORK( I+3*N ) = ZERO
396: WORK( I+4*N ) = ZERO
397: WORK( I+5*N ) = ZERO
398: 200 CONTINUE
399: *
400: * Compute right side vector in resulting linear equations
401: *
402: BASL = LOG10( SCLFAC )
403: DO 240 I = ILO, IHI
404: DO 230 J = ILO, IHI
405: IF( A( I, J ).EQ.CZERO ) THEN
406: TA = ZERO
407: GO TO 210
408: END IF
409: TA = LOG10( CABS1( A( I, J ) ) ) / BASL
410: *
411: 210 CONTINUE
412: IF( B( I, J ).EQ.CZERO ) THEN
413: TB = ZERO
414: GO TO 220
415: END IF
416: TB = LOG10( CABS1( B( I, J ) ) ) / BASL
417: *
418: 220 CONTINUE
419: WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
420: WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
421: 230 CONTINUE
422: 240 CONTINUE
423: *
424: COEF = ONE / DBLE( 2*NR )
425: COEF2 = COEF*COEF
426: COEF5 = HALF*COEF2
427: NRP2 = NR + 2
428: BETA = ZERO
429: IT = 1
430: *
431: * Start generalized conjugate gradient iteration
432: *
433: 250 CONTINUE
434: *
435: GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
436: $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
437: *
438: EW = ZERO
439: EWC = ZERO
440: DO 260 I = ILO, IHI
441: EW = EW + WORK( I+4*N )
442: EWC = EWC + WORK( I+5*N )
443: 260 CONTINUE
444: *
445: GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
446: IF( GAMMA.EQ.ZERO )
447: $ GO TO 350
448: IF( IT.NE.1 )
449: $ BETA = GAMMA / PGAMMA
450: T = COEF5*( EWC-THREE*EW )
451: TC = COEF5*( EW-THREE*EWC )
452: *
453: CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
454: CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
455: *
456: CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
457: CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
458: *
459: DO 270 I = ILO, IHI
460: WORK( I ) = WORK( I ) + TC
461: WORK( I+N ) = WORK( I+N ) + T
462: 270 CONTINUE
463: *
464: * Apply matrix to vector
465: *
466: DO 300 I = ILO, IHI
467: KOUNT = 0
468: SUM = ZERO
469: DO 290 J = ILO, IHI
470: IF( A( I, J ).EQ.CZERO )
471: $ GO TO 280
472: KOUNT = KOUNT + 1
473: SUM = SUM + WORK( J )
474: 280 CONTINUE
475: IF( B( I, J ).EQ.CZERO )
476: $ GO TO 290
477: KOUNT = KOUNT + 1
478: SUM = SUM + WORK( J )
479: 290 CONTINUE
480: WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
481: 300 CONTINUE
482: *
483: DO 330 J = ILO, IHI
484: KOUNT = 0
485: SUM = ZERO
486: DO 320 I = ILO, IHI
487: IF( A( I, J ).EQ.CZERO )
488: $ GO TO 310
489: KOUNT = KOUNT + 1
490: SUM = SUM + WORK( I+N )
491: 310 CONTINUE
492: IF( B( I, J ).EQ.CZERO )
493: $ GO TO 320
494: KOUNT = KOUNT + 1
495: SUM = SUM + WORK( I+N )
496: 320 CONTINUE
497: WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
498: 330 CONTINUE
499: *
500: SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
501: $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
502: ALPHA = GAMMA / SUM
503: *
504: * Determine correction to current iteration
505: *
506: CMAX = ZERO
507: DO 340 I = ILO, IHI
508: COR = ALPHA*WORK( I+N )
509: IF( ABS( COR ).GT.CMAX )
510: $ CMAX = ABS( COR )
511: LSCALE( I ) = LSCALE( I ) + COR
512: COR = ALPHA*WORK( I )
513: IF( ABS( COR ).GT.CMAX )
514: $ CMAX = ABS( COR )
515: RSCALE( I ) = RSCALE( I ) + COR
516: 340 CONTINUE
517: IF( CMAX.LT.HALF )
518: $ GO TO 350
519: *
520: CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
521: CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
522: *
523: PGAMMA = GAMMA
524: IT = IT + 1
525: IF( IT.LE.NRP2 )
526: $ GO TO 250
527: *
528: * End generalized conjugate gradient iteration
529: *
530: 350 CONTINUE
531: SFMIN = DLAMCH( 'S' )
532: SFMAX = ONE / SFMIN
533: LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
534: LSFMAX = INT( LOG10( SFMAX ) / BASL )
535: DO 360 I = ILO, IHI
536: IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
537: RAB = ABS( A( I, IRAB+ILO-1 ) )
538: IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB )
539: RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
540: LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
541: IR = INT(LSCALE( I ) + SIGN( HALF, LSCALE( I ) ))
542: IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
543: LSCALE( I ) = SCLFAC**IR
544: ICAB = IZAMAX( IHI, A( 1, I ), 1 )
545: CAB = ABS( A( ICAB, I ) )
546: ICAB = IZAMAX( IHI, B( 1, I ), 1 )
547: CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
548: LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
549: JC = INT(RSCALE( I ) + SIGN( HALF, RSCALE( I ) ))
550: JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
551: RSCALE( I ) = SCLFAC**JC
552: 360 CONTINUE
553: *
554: * Row scaling of matrices A and B
555: *
556: DO 370 I = ILO, IHI
557: CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
558: CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
559: 370 CONTINUE
560: *
561: * Column scaling of matrices A and B
562: *
563: DO 380 J = ILO, IHI
564: CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
565: CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
566: 380 CONTINUE
567: *
568: RETURN
569: *
570: * End of ZGGBAL
571: *
572: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggbal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack