summaryrefslogtreecommitdiff
path: root/2.3-1/src/fortran/lapack/dggbal.f
diff options
context:
space:
mode:
Diffstat (limited to '2.3-1/src/fortran/lapack/dggbal.f')
-rw-r--r--2.3-1/src/fortran/lapack/dggbal.f469
1 files changed, 469 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dggbal.f b/2.3-1/src/fortran/lapack/dggbal.f
new file mode 100644
index 00000000..2034880a
--- /dev/null
+++ b/2.3-1/src/fortran/lapack/dggbal.f
@@ -0,0 +1,469 @@
+ SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
+ $ RSCALE, WORK, INFO )
+*
+* -- LAPACK routine (version 3.1) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ CHARACTER JOB
+ INTEGER IHI, ILO, INFO, LDA, LDB, N
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
+ $ RSCALE( * ), WORK( * )
+* ..
+*
+* Purpose
+* =======
+*
+* DGGBAL balances a pair of general real matrices (A,B). This
+* involves, first, permuting A and B by similarity transformations to
+* isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
+* elements on the diagonal; and second, applying a diagonal similarity
+* transformation to rows and columns ILO to IHI to make the rows
+* and columns as close in norm as possible. Both steps are optional.
+*
+* Balancing may reduce the 1-norm of the matrices, and improve the
+* accuracy of the computed eigenvalues and/or eigenvectors in the
+* generalized eigenvalue problem A*x = lambda*B*x.
+*
+* Arguments
+* =========
+*
+* JOB (input) CHARACTER*1
+* Specifies the operations to be performed on A and B:
+* = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
+* and RSCALE(I) = 1.0 for i = 1,...,N.
+* = 'P': permute only;
+* = 'S': scale only;
+* = 'B': both permute and scale.
+*
+* N (input) INTEGER
+* The order of the matrices A and B. N >= 0.
+*
+* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+* On entry, the input matrix A.
+* On exit, A is overwritten by the balanced matrix.
+* If JOB = 'N', A is not referenced.
+*
+* LDA (input) INTEGER
+* The leading dimension of the array A. LDA >= max(1,N).
+*
+* B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
+* On entry, the input matrix B.
+* On exit, B is overwritten by the balanced matrix.
+* If JOB = 'N', B is not referenced.
+*
+* LDB (input) INTEGER
+* The leading dimension of the array B. LDB >= max(1,N).
+*
+* ILO (output) INTEGER
+* IHI (output) INTEGER
+* ILO and IHI are set to integers such that on exit
+* A(i,j) = 0 and B(i,j) = 0 if i > j and
+* j = 1,...,ILO-1 or i = IHI+1,...,N.
+* If JOB = 'N' or 'S', ILO = 1 and IHI = N.
+*
+* LSCALE (output) DOUBLE PRECISION array, dimension (N)
+* Details of the permutations and scaling factors applied
+* to the left side of A and B. If P(j) is the index of the
+* row interchanged with row j, and D(j)
+* is the scaling factor applied to row j, then
+* LSCALE(j) = P(j) for J = 1,...,ILO-1
+* = D(j) for J = ILO,...,IHI
+* = P(j) for J = IHI+1,...,N.
+* The order in which the interchanges are made is N to IHI+1,
+* then 1 to ILO-1.
+*
+* RSCALE (output) DOUBLE PRECISION array, dimension (N)
+* Details of the permutations and scaling factors applied
+* to the right side of A and B. If P(j) is the index of the
+* column interchanged with column j, and D(j)
+* is the scaling factor applied to column j, then
+* LSCALE(j) = P(j) for J = 1,...,ILO-1
+* = D(j) for J = ILO,...,IHI
+* = P(j) for J = IHI+1,...,N.
+* The order in which the interchanges are made is N to IHI+1,
+* then 1 to ILO-1.
+*
+* WORK (workspace) REAL array, dimension (lwork)
+* lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
+* at least 1 when JOB = 'N' or 'P'.
+*
+* INFO (output) INTEGER
+* = 0: successful exit
+* < 0: if INFO = -i, the i-th argument had an illegal value.
+*
+* Further Details
+* ===============
+*
+* See R.C. WARD, Balancing the generalized eigenvalue problem,
+* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, HALF, ONE
+ PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
+ DOUBLE PRECISION THREE, SCLFAC
+ PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
+* ..
+* .. Local Scalars ..
+ INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
+ $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
+ $ M, NR, NRP2
+ DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
+ $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
+ $ SFMIN, SUM, T, TA, TB, TC
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ DOUBLE PRECISION DDOT, DLAMCH
+ EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters
+*
+ INFO = 0
+ IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
+ $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -4
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DGGBAL', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 ) THEN
+ ILO = 1
+ IHI = N
+ RETURN
+ END IF
+*
+ IF( N.EQ.1 ) THEN
+ ILO = 1
+ IHI = N
+ LSCALE( 1 ) = ONE
+ RSCALE( 1 ) = ONE
+ RETURN
+ END IF
+*
+ IF( LSAME( JOB, 'N' ) ) THEN
+ ILO = 1
+ IHI = N
+ DO 10 I = 1, N
+ LSCALE( I ) = ONE
+ RSCALE( I ) = ONE
+ 10 CONTINUE
+ RETURN
+ END IF
+*
+ K = 1
+ L = N
+ IF( LSAME( JOB, 'S' ) )
+ $ GO TO 190
+*
+ GO TO 30
+*
+* Permute the matrices A and B to isolate the eigenvalues.
+*
+* Find row with one nonzero in columns 1 through L
+*
+ 20 CONTINUE
+ L = LM1
+ IF( L.NE.1 )
+ $ GO TO 30
+*
+ RSCALE( 1 ) = ONE
+ LSCALE( 1 ) = ONE
+ GO TO 190
+*
+ 30 CONTINUE
+ LM1 = L - 1
+ DO 80 I = L, 1, -1
+ DO 40 J = 1, LM1
+ JP1 = J + 1
+ IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
+ $ GO TO 50
+ 40 CONTINUE
+ J = L
+ GO TO 70
+*
+ 50 CONTINUE
+ DO 60 J = JP1, L
+ IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
+ $ GO TO 80
+ 60 CONTINUE
+ J = JP1 - 1
+*
+ 70 CONTINUE
+ M = L
+ IFLOW = 1
+ GO TO 160
+ 80 CONTINUE
+ GO TO 100
+*
+* Find column with one nonzero in rows K through N
+*
+ 90 CONTINUE
+ K = K + 1
+*
+ 100 CONTINUE
+ DO 150 J = K, L
+ DO 110 I = K, LM1
+ IP1 = I + 1
+ IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
+ $ GO TO 120
+ 110 CONTINUE
+ I = L
+ GO TO 140
+ 120 CONTINUE
+ DO 130 I = IP1, L
+ IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
+ $ GO TO 150
+ 130 CONTINUE
+ I = IP1 - 1
+ 140 CONTINUE
+ M = K
+ IFLOW = 2
+ GO TO 160
+ 150 CONTINUE
+ GO TO 190
+*
+* Permute rows M and I
+*
+ 160 CONTINUE
+ LSCALE( M ) = I
+ IF( I.EQ.M )
+ $ GO TO 170
+ CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
+ CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
+*
+* Permute columns M and J
+*
+ 170 CONTINUE
+ RSCALE( M ) = J
+ IF( J.EQ.M )
+ $ GO TO 180
+ CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
+ CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
+*
+ 180 CONTINUE
+ GO TO ( 20, 90 )IFLOW
+*
+ 190 CONTINUE
+ ILO = K
+ IHI = L
+*
+ IF( LSAME( JOB, 'P' ) ) THEN
+ DO 195 I = ILO, IHI
+ LSCALE( I ) = ONE
+ RSCALE( I ) = ONE
+ 195 CONTINUE
+ RETURN
+ END IF
+*
+ IF( ILO.EQ.IHI )
+ $ RETURN
+*
+* Balance the submatrix in rows ILO to IHI.
+*
+ NR = IHI - ILO + 1
+ DO 200 I = ILO, IHI
+ RSCALE( I ) = ZERO
+ LSCALE( I ) = ZERO
+*
+ WORK( I ) = ZERO
+ WORK( I+N ) = ZERO
+ WORK( I+2*N ) = ZERO
+ WORK( I+3*N ) = ZERO
+ WORK( I+4*N ) = ZERO
+ WORK( I+5*N ) = ZERO
+ 200 CONTINUE
+*
+* Compute right side vector in resulting linear equations
+*
+ BASL = LOG10( SCLFAC )
+ DO 240 I = ILO, IHI
+ DO 230 J = ILO, IHI
+ TB = B( I, J )
+ TA = A( I, J )
+ IF( TA.EQ.ZERO )
+ $ GO TO 210
+ TA = LOG10( ABS( TA ) ) / BASL
+ 210 CONTINUE
+ IF( TB.EQ.ZERO )
+ $ GO TO 220
+ TB = LOG10( ABS( TB ) ) / BASL
+ 220 CONTINUE
+ WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
+ WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
+ 230 CONTINUE
+ 240 CONTINUE
+*
+ COEF = ONE / DBLE( 2*NR )
+ COEF2 = COEF*COEF
+ COEF5 = HALF*COEF2
+ NRP2 = NR + 2
+ BETA = ZERO
+ IT = 1
+*
+* Start generalized conjugate gradient iteration
+*
+ 250 CONTINUE
+*
+ GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
+ $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
+*
+ EW = ZERO
+ EWC = ZERO
+ DO 260 I = ILO, IHI
+ EW = EW + WORK( I+4*N )
+ EWC = EWC + WORK( I+5*N )
+ 260 CONTINUE
+*
+ GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
+ IF( GAMMA.EQ.ZERO )
+ $ GO TO 350
+ IF( IT.NE.1 )
+ $ BETA = GAMMA / PGAMMA
+ T = COEF5*( EWC-THREE*EW )
+ TC = COEF5*( EW-THREE*EWC )
+*
+ CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
+ CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
+*
+ CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
+ CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
+*
+ DO 270 I = ILO, IHI
+ WORK( I ) = WORK( I ) + TC
+ WORK( I+N ) = WORK( I+N ) + T
+ 270 CONTINUE
+*
+* Apply matrix to vector
+*
+ DO 300 I = ILO, IHI
+ KOUNT = 0
+ SUM = ZERO
+ DO 290 J = ILO, IHI
+ IF( A( I, J ).EQ.ZERO )
+ $ GO TO 280
+ KOUNT = KOUNT + 1
+ SUM = SUM + WORK( J )
+ 280 CONTINUE
+ IF( B( I, J ).EQ.ZERO )
+ $ GO TO 290
+ KOUNT = KOUNT + 1
+ SUM = SUM + WORK( J )
+ 290 CONTINUE
+ WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
+ 300 CONTINUE
+*
+ DO 330 J = ILO, IHI
+ KOUNT = 0
+ SUM = ZERO
+ DO 320 I = ILO, IHI
+ IF( A( I, J ).EQ.ZERO )
+ $ GO TO 310
+ KOUNT = KOUNT + 1
+ SUM = SUM + WORK( I+N )
+ 310 CONTINUE
+ IF( B( I, J ).EQ.ZERO )
+ $ GO TO 320
+ KOUNT = KOUNT + 1
+ SUM = SUM + WORK( I+N )
+ 320 CONTINUE
+ WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
+ 330 CONTINUE
+*
+ SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
+ $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
+ ALPHA = GAMMA / SUM
+*
+* Determine correction to current iteration
+*
+ CMAX = ZERO
+ DO 340 I = ILO, IHI
+ COR = ALPHA*WORK( I+N )
+ IF( ABS( COR ).GT.CMAX )
+ $ CMAX = ABS( COR )
+ LSCALE( I ) = LSCALE( I ) + COR
+ COR = ALPHA*WORK( I )
+ IF( ABS( COR ).GT.CMAX )
+ $ CMAX = ABS( COR )
+ RSCALE( I ) = RSCALE( I ) + COR
+ 340 CONTINUE
+ IF( CMAX.LT.HALF )
+ $ GO TO 350
+*
+ CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
+ CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
+*
+ PGAMMA = GAMMA
+ IT = IT + 1
+ IF( IT.LE.NRP2 )
+ $ GO TO 250
+*
+* End generalized conjugate gradient iteration
+*
+ 350 CONTINUE
+ SFMIN = DLAMCH( 'S' )
+ SFMAX = ONE / SFMIN
+ LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
+ LSFMAX = INT( LOG10( SFMAX ) / BASL )
+ DO 360 I = ILO, IHI
+ IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
+ RAB = ABS( A( I, IRAB+ILO-1 ) )
+ IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
+ RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
+ LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
+ IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
+ IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
+ LSCALE( I ) = SCLFAC**IR
+ ICAB = IDAMAX( IHI, A( 1, I ), 1 )
+ CAB = ABS( A( ICAB, I ) )
+ ICAB = IDAMAX( IHI, B( 1, I ), 1 )
+ CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
+ LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
+ JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
+ JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
+ RSCALE( I ) = SCLFAC**JC
+ 360 CONTINUE
+*
+* Row scaling of matrices A and B
+*
+ DO 370 I = ILO, IHI
+ CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
+ CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
+ 370 CONTINUE
+*
+* Column scaling of matrices A and B
+*
+ DO 380 J = ILO, IHI
+ CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
+ CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
+ 380 CONTINUE
+*
+ RETURN
+*
+* End of DGGBAL
+*
+ END
generated by cgit (git 2.34.1) at 2024-12-20 11:35:51 +0000