2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480 DOUBLE PRECISION ZERO, ONE
2481 parameter( zero = 0.0d0, one = 1.0d0 )
2482 DOUBLE PRECISION ROGUE
2483 parameter( rogue = -1.0d10 )
2484
2485 DOUBLE PRECISION TRANSL
2486 INTEGER KL, KU, LDA, M, N, NMAX
2487 LOGICAL RESET
2488 CHARACTER*1 DIAG, UPLO
2489 CHARACTER*2 TYPE
2490
2491 DOUBLE PRECISION A( NMAX, * ), AA( * )
2492
2493 INTEGER I, I1, I2, I3, IBEG, IEND, IOFF, J, KK
2494 LOGICAL GEN, LOWER, SYM, TRI, UNIT, UPPER
2495
2496 DOUBLE PRECISION DBEG
2498
2499 INTRINSIC max, min
2500
2501 gen = TYPE( 1: 1 ).EQ.'g'
2502 sym = TYPE( 1: 1 ).EQ.'s'
2503 tri = TYPE( 1: 1 ).EQ.'t'
2504 upper = ( sym.OR.tri ).AND.uplo.EQ.'U'
2505 lower = ( sym.OR.tri ).AND.uplo.EQ.'L'
2506 unit = tri.AND.diag.EQ.'U'
2507
2508
2509
2510 DO 20 j = 1, n
2511 DO 10 i = 1, m
2512 IF( gen.OR.( upper.AND.i.LE.j ).OR.( lower.AND.i.GE.j ) )
2513 $ THEN
2514 IF( ( i.LE.j.AND.j - i.LE.ku ).OR.
2515 $ ( i.GE.j.AND.i - j.LE.kl ) )THEN
2516 a( i, j ) =
dbeg( reset ) + transl
2517 ELSE
2518 a( i, j ) = zero
2519 END IF
2520 IF( i.NE.j )THEN
2521 IF( sym )THEN
2522 a( j, i ) = a( i, j )
2523 ELSE IF( tri )THEN
2524 a( j, i ) = zero
2525 END IF
2526 END IF
2527 END IF
2528 10 CONTINUE
2529 IF( tri )
2530 $ a( j, j ) = a( j, j ) + one
2531 IF( unit )
2532 $ a( j, j ) = one
2533 20 CONTINUE
2534
2535
2536
2537 IF( type.EQ.'ge' )THEN
2538 DO 50 j = 1, n
2539 DO 30 i = 1, m
2540 aa( i + ( j - 1 )*lda ) = a( i, j )
2541 30 CONTINUE
2542 DO 40 i = m + 1, lda
2543 aa( i + ( j - 1 )*lda ) = rogue
2544 40 CONTINUE
2545 50 CONTINUE
2546 ELSE IF( type.EQ.'gb' )THEN
2547 DO 90 j = 1, n
2548 DO 60 i1 = 1, ku + 1 - j
2549 aa( i1 + ( j - 1 )*lda ) = rogue
2550 60 CONTINUE
2551 DO 70 i2 = i1, min( kl + ku + 1, ku + 1 + m - j )
2552 aa( i2 + ( j - 1 )*lda ) = a( i2 + j - ku - 1, j )
2553 70 CONTINUE
2554 DO 80 i3 = i2, lda
2555 aa( i3 + ( j - 1 )*lda ) = rogue
2556 80 CONTINUE
2557 90 CONTINUE
2558 ELSE IF( type.EQ.'sy'.OR.type.EQ.'tr' )THEN
2559 DO 130 j = 1, n
2560 IF( upper )THEN
2561 ibeg = 1
2562 IF( unit )THEN
2563 iend = j - 1
2564 ELSE
2565 iend = j
2566 END IF
2567 ELSE
2568 IF( unit )THEN
2569 ibeg = j + 1
2570 ELSE
2571 ibeg = j
2572 END IF
2573 iend = n
2574 END IF
2575 DO 100 i = 1, ibeg - 1
2576 aa( i + ( j - 1 )*lda ) = rogue
2577 100 CONTINUE
2578 DO 110 i = ibeg, iend
2579 aa( i + ( j - 1 )*lda ) = a( i, j )
2580 110 CONTINUE
2581 DO 120 i = iend + 1, lda
2582 aa( i + ( j - 1 )*lda ) = rogue
2583 120 CONTINUE
2584 130 CONTINUE
2585 ELSE IF( type.EQ.'sb'.OR.type.EQ.'tb' )THEN
2586 DO 170 j = 1, n
2587 IF( upper )THEN
2588 kk = kl + 1
2589 ibeg = max( 1, kl + 2 - j )
2590 IF( unit )THEN
2591 iend = kl
2592 ELSE
2593 iend = kl + 1
2594 END IF
2595 ELSE
2596 kk = 1
2597 IF( unit )THEN
2598 ibeg = 2
2599 ELSE
2600 ibeg = 1
2601 END IF
2602 iend = min( kl + 1, 1 + m - j )
2603 END IF
2604 DO 140 i = 1, ibeg - 1
2605 aa( i + ( j - 1 )*lda ) = rogue
2606 140 CONTINUE
2607 DO 150 i = ibeg, iend
2608 aa( i + ( j - 1 )*lda ) = a( i + j - kk, j )
2609 150 CONTINUE
2610 DO 160 i = iend + 1, lda
2611 aa( i + ( j - 1 )*lda ) = rogue
2612 160 CONTINUE
2613 170 CONTINUE
2614 ELSE IF( type.EQ.'sp'.OR.type.EQ.'tp' )THEN
2615 ioff = 0
2616 DO 190 j = 1, n
2617 IF( upper )THEN
2618 ibeg = 1
2619 iend = j
2620 ELSE
2621 ibeg = j
2622 iend = n
2623 END IF
2624 DO 180 i = ibeg, iend
2625 ioff = ioff + 1
2626 aa( ioff ) = a( i, j )
2627 IF( i.EQ.j )THEN
2628 IF( unit )
2629 $ aa( ioff ) = rogue
2630 END IF
2631 180 CONTINUE
2632 190 CONTINUE
2633 END IF
2634 RETURN
2635
2636
2637
double precision function dbeg(reset)