177
178
179
180
181
182
183 CHARACTER JOB
184 INTEGER IHI, ILO, INFO, LDA, LDB, N
185
186
187 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
188 $ RSCALE( * ), WORK( * )
189
190
191
192
193
194 DOUBLE PRECISION ZERO, HALF, ONE
195 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
196 DOUBLE PRECISION THREE, SCLFAC
197 parameter( three = 3.0d+0, sclfac = 1.0d+1 )
198
199
200 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
201 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
202 $ M, NR, NRP2
203 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
204 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
205 $ SFMIN, SUM, T, TA, TB, TC
206
207
208 LOGICAL LSAME
209 INTEGER IDAMAX
210 DOUBLE PRECISION DDOT, DLAMCH
212
213
215
216
217 INTRINSIC abs, dble, int, log10, max, min, sign
218
219
220
221
222
223 info = 0
224 IF( .NOT.
lsame( job,
'N' ) .AND. .NOT.
lsame( job,
'P' ) .AND.
225 $ .NOT.
lsame( job,
'S' ) .AND. .NOT.
lsame( job,
'B' ) )
THEN
226 info = -1
227 ELSE IF( n.LT.0 ) THEN
228 info = -2
229 ELSE IF( lda.LT.max( 1, n ) ) THEN
230 info = -4
231 ELSE IF( ldb.LT.max( 1, n ) ) THEN
232 info = -6
233 END IF
234 IF( info.NE.0 ) THEN
235 CALL xerbla(
'DGGBAL', -info )
236 RETURN
237 END IF
238
239
240
241 IF( n.EQ.0 ) THEN
242 ilo = 1
243 ihi = n
244 RETURN
245 END IF
246
247 IF( n.EQ.1 ) THEN
248 ilo = 1
249 ihi = n
250 lscale( 1 ) = one
251 rscale( 1 ) = one
252 RETURN
253 END IF
254
255 IF(
lsame( job,
'N' ) )
THEN
256 ilo = 1
257 ihi = n
258 DO 10 i = 1, n
259 lscale( i ) = one
260 rscale( i ) = one
261 10 CONTINUE
262 RETURN
263 END IF
264
265 k = 1
266 l = n
267 IF(
lsame( job,
'S' ) )
268 $ GO TO 190
269
270 GO TO 30
271
272
273
274
275
276 20 CONTINUE
277 l = lm1
278 IF( l.NE.1 )
279 $ GO TO 30
280
281 rscale( 1 ) = one
282 lscale( 1 ) = one
283 GO TO 190
284
285 30 CONTINUE
286 lm1 = l - 1
287 DO 80 i = l, 1, -1
288 DO 40 j = 1, lm1
289 jp1 = j + 1
290 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
291 $ GO TO 50
292 40 CONTINUE
293 j = l
294 GO TO 70
295
296 50 CONTINUE
297 DO 60 j = jp1, l
298 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
299 $ GO TO 80
300 60 CONTINUE
301 j = jp1 - 1
302
303 70 CONTINUE
304 m = l
305 iflow = 1
306 GO TO 160
307 80 CONTINUE
308 GO TO 100
309
310
311
312 90 CONTINUE
313 k = k + 1
314
315 100 CONTINUE
316 DO 150 j = k, l
317 DO 110 i = k, lm1
318 ip1 = i + 1
319 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
320 $ GO TO 120
321 110 CONTINUE
322 i = l
323 GO TO 140
324 120 CONTINUE
325 DO 130 i = ip1, l
326 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
327 $ GO TO 150
328 130 CONTINUE
329 i = ip1 - 1
330 140 CONTINUE
331 m = k
332 iflow = 2
333 GO TO 160
334 150 CONTINUE
335 GO TO 190
336
337
338
339 160 CONTINUE
340 lscale( m ) = i
341 IF( i.EQ.m )
342 $ GO TO 170
343 CALL dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344 CALL dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
345
346
347
348 170 CONTINUE
349 rscale( m ) = j
350 IF( j.EQ.m )
351 $ GO TO 180
352 CALL dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353 CALL dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
354
355 180 CONTINUE
356 GO TO ( 20, 90 )iflow
357
358 190 CONTINUE
359 ilo = k
360 ihi = l
361
362 IF(
lsame( job,
'P' ) )
THEN
363 DO 195 i = ilo, ihi
364 lscale( i ) = one
365 rscale( i ) = one
366 195 CONTINUE
367 RETURN
368 END IF
369
370 IF( ilo.EQ.ihi )
371 $ RETURN
372
373
374
375 nr = ihi - ilo + 1
376 DO 200 i = ilo, ihi
377 rscale( i ) = zero
378 lscale( i ) = zero
379
380 work( i ) = zero
381 work( i+n ) = zero
382 work( i+2*n ) = zero
383 work( i+3*n ) = zero
384 work( i+4*n ) = zero
385 work( i+5*n ) = zero
386 200 CONTINUE
387
388
389
390 basl = log10( sclfac )
391 DO 240 i = ilo, ihi
392 DO 230 j = ilo, ihi
393 tb = b( i, j )
394 ta = a( i, j )
395 IF( ta.EQ.zero )
396 $ GO TO 210
397 ta = log10( abs( ta ) ) / basl
398 210 CONTINUE
399 IF( tb.EQ.zero )
400 $ GO TO 220
401 tb = log10( abs( tb ) ) / basl
402 220 CONTINUE
403 work( i+4*n ) = work( i+4*n ) - ta - tb
404 work( j+5*n ) = work( j+5*n ) - ta - tb
405 230 CONTINUE
406 240 CONTINUE
407
408 coef = one / dble( 2*nr )
409 coef2 = coef*coef
410 coef5 = half*coef2
411 nrp2 = nr + 2
412 beta = zero
413 it = 1
414
415
416
417 250 CONTINUE
418
419 gamma =
ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420 $
ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
421
422 ew = zero
423 ewc = zero
424 DO 260 i = ilo, ihi
425 ew = ew + work( i+4*n )
426 ewc = ewc + work( i+5*n )
427 260 CONTINUE
428
429 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
430 IF( gamma.EQ.zero )
431 $ GO TO 350
432 IF( it.NE.1 )
433 $ beta = gamma / pgamma
434 t = coef5*( ewc-three*ew )
435 tc = coef5*( ew-three*ewc )
436
437 CALL dscal( nr, beta, work( ilo ), 1 )
438 CALL dscal( nr, beta, work( ilo+n ), 1 )
439
440 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
442
443 DO 270 i = ilo, ihi
444 work( i ) = work( i ) + tc
445 work( i+n ) = work( i+n ) + t
446 270 CONTINUE
447
448
449
450 DO 300 i = ilo, ihi
451 kount = 0
452 sum = zero
453 DO 290 j = ilo, ihi
454 IF( a( i, j ).EQ.zero )
455 $ GO TO 280
456 kount = kount + 1
457 sum = sum + work( j )
458 280 CONTINUE
459 IF( b( i, j ).EQ.zero )
460 $ GO TO 290
461 kount = kount + 1
462 sum = sum + work( j )
463 290 CONTINUE
464 work( i+2*n ) = dble( kount )*work( i+n ) + sum
465 300 CONTINUE
466
467 DO 330 j = ilo, ihi
468 kount = 0
469 sum = zero
470 DO 320 i = ilo, ihi
471 IF( a( i, j ).EQ.zero )
472 $ GO TO 310
473 kount = kount + 1
474 sum = sum + work( i+n )
475 310 CONTINUE
476 IF( b( i, j ).EQ.zero )
477 $ GO TO 320
478 kount = kount + 1
479 sum = sum + work( i+n )
480 320 CONTINUE
481 work( j+3*n ) = dble( kount )*work( j ) + sum
482 330 CONTINUE
483
484 sum =
ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485 $
ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
486 alpha = gamma / sum
487
488
489
490 cmax = zero
491 DO 340 i = ilo, ihi
492 cor = alpha*work( i+n )
493 IF( abs( cor ).GT.cmax )
494 $ cmax = abs( cor )
495 lscale( i ) = lscale( i ) + cor
496 cor = alpha*work( i )
497 IF( abs( cor ).GT.cmax )
498 $ cmax = abs( cor )
499 rscale( i ) = rscale( i ) + cor
500 340 CONTINUE
501 IF( cmax.LT.half )
502 $ GO TO 350
503
504 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
505 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
506
507 pgamma = gamma
508 it = it + 1
509 IF( it.LE.nrp2 )
510 $ GO TO 250
511
512
513
514 350 CONTINUE
516 sfmax = one / sfmin
517 lsfmin = int( log10( sfmin ) / basl+one )
518 lsfmax = int( log10( sfmax ) / basl )
519 DO 360 i = ilo, ihi
520 irab =
idamax( n-ilo+1, a( i, ilo ), lda )
521 rab = abs( a( i, irab+ilo-1 ) )
522 irab =
idamax( n-ilo+1, b( i, ilo ), ldb )
523 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
524 lrab = int( log10( rab+sfmin ) / basl+one )
525 ir = int(lscale( i ) + sign( half, lscale( i ) ))
526 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
527 lscale( i ) = sclfac**ir
528 icab =
idamax( ihi, a( 1, i ), 1 )
529 cab = abs( a( icab, i ) )
530 icab =
idamax( ihi, b( 1, i ), 1 )
531 cab = max( cab, abs( b( icab, i ) ) )
532 lcab = int( log10( cab+sfmin ) / basl+one )
533 jc = int(rscale( i ) + sign( half, rscale( i ) ))
534 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
535 rscale( i ) = sclfac**jc
536 360 CONTINUE
537
538
539
540 DO 370 i = ilo, ihi
541 CALL dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
542 CALL dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
543 370 CONTINUE
544
545
546
547 DO 380 j = ilo, ihi
548 CALL dscal( ihi, rscale( j ), a( 1, j ), 1 )
549 CALL dscal( ihi, rscale( j ), b( 1, j ), 1 )
550 380 CONTINUE
551
552 RETURN
553
554
555
subroutine xerbla(srname, info)
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
double precision function ddot(n, dx, incx, dy, incy)
DDOT
integer function idamax(n, dx, incx)
IDAMAX
double precision function dlamch(cmach)
DLAMCH
logical function lsame(ca, cb)
LSAME
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine dswap(n, dx, incx, dy, incy)
DSWAP