53
54
55
56
57
58
59 INTEGER NIN, NOUT
60
61
62
63
64
65 INTEGER LDA, LDB, LWORK
66 parameter( lda = 20, ldb = 20, lwork = 6*lda )
67 DOUBLE PRECISION ZERO
68 parameter( zero = 0.0d+0 )
69
70
71 INTEGER I, IHI, IHIIN, ILO, ILOIN, INFO, J, KNT, N,
72 $ NINFO
73 DOUBLE PRECISION ANORM, BNORM, EPS, RMAX, VMAX
74
75
76 INTEGER LMAX( 5 )
77 DOUBLE PRECISION A( LDA, LDA ), AIN( LDA, LDA ), B( LDB, LDB ),
78 $ BIN( LDB, LDB ), LSCALE( LDA ), LSCLIN( LDA ),
79 $ RSCALE( LDA ), RSCLIN( LDA ), WORK( LWORK )
80
81
82 DOUBLE PRECISION DLAMCH, DLANGE
84
85
87
88
89 INTRINSIC abs, max
90
91
92
93 lmax( 1 ) = 0
94 lmax( 2 ) = 0
95 lmax( 3 ) = 0
96 ninfo = 0
97 knt = 0
98 rmax = zero
99
100 eps =
dlamch(
'Precision' )
101
102 10 CONTINUE
103
104 READ( nin, fmt = * )n
105 IF( n.EQ.0 )
106 $ GO TO 90
107 DO 20 i = 1, n
108 READ( nin, fmt = * )( a( i, j ), j = 1, n )
109 20 CONTINUE
110
111 DO 30 i = 1, n
112 READ( nin, fmt = * )( b( i, j ), j = 1, n )
113 30 CONTINUE
114
115 READ( nin, fmt = * )iloin, ihiin
116 DO 40 i = 1, n
117 READ( nin, fmt = * )( ain( i, j ), j = 1, n )
118 40 CONTINUE
119 DO 50 i = 1, n
120 READ( nin, fmt = * )( bin( i, j ), j = 1, n )
121 50 CONTINUE
122
123 READ( nin, fmt = * )( lsclin( i ), i = 1, n )
124 READ( nin, fmt = * )( rsclin( i ), i = 1, n )
125
126 anorm =
dlange(
'M', n, n, a, lda, work )
127 bnorm =
dlange(
'M', n, n, b, ldb, work )
128
129 knt = knt + 1
130
131 CALL dggbal(
'B', n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
132 $ work, info )
133
134 IF( info.NE.0 ) THEN
135 ninfo = ninfo + 1
136 lmax( 1 ) = knt
137 END IF
138
139 IF( ilo.NE.iloin .OR. ihi.NE.ihiin ) THEN
140 ninfo = ninfo + 1
141 lmax( 2 ) = knt
142 END IF
143
144 vmax = zero
145 DO 70 i = 1, n
146 DO 60 j = 1, n
147 vmax = max( vmax, abs( a( i, j )-ain( i, j ) ) )
148 vmax = max( vmax, abs( b( i, j )-bin( i, j ) ) )
149 60 CONTINUE
150 70 CONTINUE
151
152 DO 80 i = 1, n
153 vmax = max( vmax, abs( lscale( i )-lsclin( i ) ) )
154 vmax = max( vmax, abs( rscale( i )-rsclin( i ) ) )
155 80 CONTINUE
156
157 vmax = vmax / ( eps*max( anorm, bnorm ) )
158
159 IF( vmax.GT.rmax ) THEN
160 lmax( 3 ) = knt
161 rmax = vmax
162 END IF
163
164 GO TO 10
165
166 90 CONTINUE
167
168 WRITE( nout, fmt = 9999 )
169 9999 FORMAT( 1x, '.. test output of DGGBAL .. ' )
170
171 WRITE( nout, fmt = 9998 )rmax
172 9998 FORMAT( 1x, 'value of largest test error = ', d12.3 )
173 WRITE( nout, fmt = 9997 )lmax( 1 )
174 9997 FORMAT( 1x, 'example number where info is not zero = ', i4 )
175 WRITE( nout, fmt = 9996 )lmax( 2 )
176 9996 FORMAT( 1x, 'example number where ILO or IHI wrong = ', i4 )
177 WRITE( nout, fmt = 9995 )lmax( 3 )
178 9995 FORMAT( 1x, 'example number having largest error = ', i4 )
179 WRITE( nout, fmt = 9994 )ninfo
180 9994 FORMAT( 1x, 'number of examples where info is not 0 = ', i4 )
181 WRITE( nout, fmt = 9993 )knt
182 9993 FORMAT( 1x, 'total number of examples tested = ', i4 )
183
184 RETURN
185
186
187
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...