Annotation of rpl/lapack/lapack/dggbal.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DGGBAL
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGGBAL + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggbal.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggbal.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggbal.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGGBAL( 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 A( LDA, * ), B( LDB, * ), LSCALE( * ),
! 30: * $ RSCALE( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DGGBAL balances a pair of general real 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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)
! 119: *> is the scaling factor 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)
! 133: *> is the scaling factor applied to column j, then
! 134: *> LSCALE(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 November 2011
! 164: *
! 165: *> \ingroup doubleGBcomputational
! 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: * =====================================================================
1.1 bertrand 177: SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
178: $ RSCALE, WORK, INFO )
179: *
1.9 ! bertrand 180: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 181: * -- LAPACK is a software package provided by Univ. of Tennessee, --
182: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 183: * November 2011
1.1 bertrand 184: *
185: * .. Scalar Arguments ..
186: CHARACTER JOB
187: INTEGER IHI, ILO, INFO, LDA, LDB, N
188: * ..
189: * .. Array Arguments ..
190: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
191: $ RSCALE( * ), WORK( * )
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: * ..
202: * .. Local Scalars ..
203: INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
204: $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
205: $ M, NR, NRP2
206: DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
207: $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
208: $ SFMIN, SUM, T, TA, TB, TC
209: * ..
210: * .. External Functions ..
211: LOGICAL LSAME
212: INTEGER IDAMAX
213: DOUBLE PRECISION DDOT, DLAMCH
214: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
215: * ..
216: * .. External Subroutines ..
217: EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA
218: * ..
219: * .. Intrinsic Functions ..
220: INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
221: * ..
222: * .. Executable Statements ..
223: *
224: * Test the input parameters
225: *
226: INFO = 0
227: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
228: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
229: INFO = -1
230: ELSE IF( N.LT.0 ) THEN
231: INFO = -2
232: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
233: INFO = -4
234: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
235: INFO = -6
236: END IF
237: IF( INFO.NE.0 ) THEN
238: CALL XERBLA( 'DGGBAL', -INFO )
239: RETURN
240: END IF
241: *
242: * Quick return if possible
243: *
244: IF( N.EQ.0 ) THEN
245: ILO = 1
246: IHI = N
247: RETURN
248: END IF
249: *
250: IF( N.EQ.1 ) THEN
251: ILO = 1
252: IHI = N
253: LSCALE( 1 ) = ONE
254: RSCALE( 1 ) = ONE
255: RETURN
256: END IF
257: *
258: IF( LSAME( JOB, 'N' ) ) THEN
259: ILO = 1
260: IHI = N
261: DO 10 I = 1, N
262: LSCALE( I ) = ONE
263: RSCALE( I ) = ONE
264: 10 CONTINUE
265: RETURN
266: END IF
267: *
268: K = 1
269: L = N
270: IF( LSAME( JOB, 'S' ) )
271: $ GO TO 190
272: *
273: GO TO 30
274: *
275: * Permute the matrices A and B to isolate the eigenvalues.
276: *
277: * Find row with one nonzero in columns 1 through L
278: *
279: 20 CONTINUE
280: L = LM1
281: IF( L.NE.1 )
282: $ GO TO 30
283: *
284: RSCALE( 1 ) = ONE
285: LSCALE( 1 ) = ONE
286: GO TO 190
287: *
288: 30 CONTINUE
289: LM1 = L - 1
290: DO 80 I = L, 1, -1
291: DO 40 J = 1, LM1
292: JP1 = J + 1
293: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
294: $ GO TO 50
295: 40 CONTINUE
296: J = L
297: GO TO 70
298: *
299: 50 CONTINUE
300: DO 60 J = JP1, L
301: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
302: $ GO TO 80
303: 60 CONTINUE
304: J = JP1 - 1
305: *
306: 70 CONTINUE
307: M = L
308: IFLOW = 1
309: GO TO 160
310: 80 CONTINUE
311: GO TO 100
312: *
313: * Find column with one nonzero in rows K through N
314: *
315: 90 CONTINUE
316: K = K + 1
317: *
318: 100 CONTINUE
319: DO 150 J = K, L
320: DO 110 I = K, LM1
321: IP1 = I + 1
322: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
323: $ GO TO 120
324: 110 CONTINUE
325: I = L
326: GO TO 140
327: 120 CONTINUE
328: DO 130 I = IP1, L
329: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
330: $ GO TO 150
331: 130 CONTINUE
332: I = IP1 - 1
333: 140 CONTINUE
334: M = K
335: IFLOW = 2
336: GO TO 160
337: 150 CONTINUE
338: GO TO 190
339: *
340: * Permute rows M and I
341: *
342: 160 CONTINUE
343: LSCALE( M ) = I
344: IF( I.EQ.M )
345: $ GO TO 170
346: CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
347: CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
348: *
349: * Permute columns M and J
350: *
351: 170 CONTINUE
352: RSCALE( M ) = J
353: IF( J.EQ.M )
354: $ GO TO 180
355: CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
356: CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
357: *
358: 180 CONTINUE
359: GO TO ( 20, 90 )IFLOW
360: *
361: 190 CONTINUE
362: ILO = K
363: IHI = L
364: *
365: IF( LSAME( JOB, 'P' ) ) THEN
366: DO 195 I = ILO, IHI
367: LSCALE( I ) = ONE
368: RSCALE( I ) = ONE
369: 195 CONTINUE
370: RETURN
371: END IF
372: *
373: IF( ILO.EQ.IHI )
374: $ RETURN
375: *
376: * Balance the submatrix in rows ILO to IHI.
377: *
378: NR = IHI - ILO + 1
379: DO 200 I = ILO, IHI
380: RSCALE( I ) = ZERO
381: LSCALE( I ) = ZERO
382: *
383: WORK( I ) = ZERO
384: WORK( I+N ) = ZERO
385: WORK( I+2*N ) = ZERO
386: WORK( I+3*N ) = ZERO
387: WORK( I+4*N ) = ZERO
388: WORK( I+5*N ) = ZERO
389: 200 CONTINUE
390: *
391: * Compute right side vector in resulting linear equations
392: *
393: BASL = LOG10( SCLFAC )
394: DO 240 I = ILO, IHI
395: DO 230 J = ILO, IHI
396: TB = B( I, J )
397: TA = A( I, J )
398: IF( TA.EQ.ZERO )
399: $ GO TO 210
400: TA = LOG10( ABS( TA ) ) / BASL
401: 210 CONTINUE
402: IF( TB.EQ.ZERO )
403: $ GO TO 220
404: TB = LOG10( ABS( TB ) ) / BASL
405: 220 CONTINUE
406: WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
407: WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
408: 230 CONTINUE
409: 240 CONTINUE
410: *
411: COEF = ONE / DBLE( 2*NR )
412: COEF2 = COEF*COEF
413: COEF5 = HALF*COEF2
414: NRP2 = NR + 2
415: BETA = ZERO
416: IT = 1
417: *
418: * Start generalized conjugate gradient iteration
419: *
420: 250 CONTINUE
421: *
422: GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
423: $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
424: *
425: EW = ZERO
426: EWC = ZERO
427: DO 260 I = ILO, IHI
428: EW = EW + WORK( I+4*N )
429: EWC = EWC + WORK( I+5*N )
430: 260 CONTINUE
431: *
432: GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
433: IF( GAMMA.EQ.ZERO )
434: $ GO TO 350
435: IF( IT.NE.1 )
436: $ BETA = GAMMA / PGAMMA
437: T = COEF5*( EWC-THREE*EW )
438: TC = COEF5*( EW-THREE*EWC )
439: *
440: CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
441: CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
442: *
443: CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
444: CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
445: *
446: DO 270 I = ILO, IHI
447: WORK( I ) = WORK( I ) + TC
448: WORK( I+N ) = WORK( I+N ) + T
449: 270 CONTINUE
450: *
451: * Apply matrix to vector
452: *
453: DO 300 I = ILO, IHI
454: KOUNT = 0
455: SUM = ZERO
456: DO 290 J = ILO, IHI
457: IF( A( I, J ).EQ.ZERO )
458: $ GO TO 280
459: KOUNT = KOUNT + 1
460: SUM = SUM + WORK( J )
461: 280 CONTINUE
462: IF( B( I, J ).EQ.ZERO )
463: $ GO TO 290
464: KOUNT = KOUNT + 1
465: SUM = SUM + WORK( J )
466: 290 CONTINUE
467: WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
468: 300 CONTINUE
469: *
470: DO 330 J = ILO, IHI
471: KOUNT = 0
472: SUM = ZERO
473: DO 320 I = ILO, IHI
474: IF( A( I, J ).EQ.ZERO )
475: $ GO TO 310
476: KOUNT = KOUNT + 1
477: SUM = SUM + WORK( I+N )
478: 310 CONTINUE
479: IF( B( I, J ).EQ.ZERO )
480: $ GO TO 320
481: KOUNT = KOUNT + 1
482: SUM = SUM + WORK( I+N )
483: 320 CONTINUE
484: WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
485: 330 CONTINUE
486: *
487: SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
488: $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
489: ALPHA = GAMMA / SUM
490: *
491: * Determine correction to current iteration
492: *
493: CMAX = ZERO
494: DO 340 I = ILO, IHI
495: COR = ALPHA*WORK( I+N )
496: IF( ABS( COR ).GT.CMAX )
497: $ CMAX = ABS( COR )
498: LSCALE( I ) = LSCALE( I ) + COR
499: COR = ALPHA*WORK( I )
500: IF( ABS( COR ).GT.CMAX )
501: $ CMAX = ABS( COR )
502: RSCALE( I ) = RSCALE( I ) + COR
503: 340 CONTINUE
504: IF( CMAX.LT.HALF )
505: $ GO TO 350
506: *
507: CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
508: CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
509: *
510: PGAMMA = GAMMA
511: IT = IT + 1
512: IF( IT.LE.NRP2 )
513: $ GO TO 250
514: *
515: * End generalized conjugate gradient iteration
516: *
517: 350 CONTINUE
518: SFMIN = DLAMCH( 'S' )
519: SFMAX = ONE / SFMIN
520: LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
521: LSFMAX = INT( LOG10( SFMAX ) / BASL )
522: DO 360 I = ILO, IHI
523: IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
524: RAB = ABS( A( I, IRAB+ILO-1 ) )
525: IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
526: RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
527: LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
528: IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
529: IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
530: LSCALE( I ) = SCLFAC**IR
531: ICAB = IDAMAX( IHI, A( 1, I ), 1 )
532: CAB = ABS( A( ICAB, I ) )
533: ICAB = IDAMAX( IHI, B( 1, I ), 1 )
534: CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
535: LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
536: JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
537: JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
538: RSCALE( I ) = SCLFAC**JC
539: 360 CONTINUE
540: *
541: * Row scaling of matrices A and B
542: *
543: DO 370 I = ILO, IHI
544: CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
545: CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
546: 370 CONTINUE
547: *
548: * Column scaling of matrices A and B
549: *
550: DO 380 J = ILO, IHI
551: CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
552: CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
553: 380 CONTINUE
554: *
555: RETURN
556: *
557: * End of DGGBAL
558: *
559: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggbal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack