2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766 COMPLEX*16 ZERO, ONE
2767 parameter( zero = ( 0.0d0, 0.0d0 ),
2768 $ one = ( 1.0d0, 0.0d0 ) )
2769 COMPLEX*16 ROGUE
2770 parameter( rogue = ( -1.0d10, 1.0d10 ) )
2771 DOUBLE PRECISION RZERO
2772 parameter( rzero = 0.0d0 )
2773 DOUBLE PRECISION RROGUE
2774 parameter( rrogue = -1.0d10 )
2775
2776 COMPLEX*16 TRANSL
2777 INTEGER KL, KU, LDA, M, N, NMAX
2778 LOGICAL RESET
2779 CHARACTER*1 DIAG, UPLO
2780 CHARACTER*2 TYPE
2781
2782 COMPLEX*16 A( NMAX, * ), AA( * )
2783
2784 INTEGER I, I1, I2, I3, IBEG, IEND, IOFF, J, JJ, KK
2785 LOGICAL GEN, LOWER, SYM, TRI, UNIT, UPPER
2786
2787 COMPLEX*16 ZBEG
2789
2790 INTRINSIC dble, dcmplx, dconjg, max, min
2791
2792 gen = TYPE( 1: 1 ).EQ.'G'
2793 sym = TYPE( 1: 1 ).EQ.'H'
2794 tri = TYPE( 1: 1 ).EQ.'T'
2795 upper = ( sym.OR.tri ).AND.uplo.EQ.'U'
2796 lower = ( sym.OR.tri ).AND.uplo.EQ.'L'
2797 unit = tri.AND.diag.EQ.'U'
2798
2799
2800
2801 DO 20 j = 1, n
2802 DO 10 i = 1, m
2803 IF( gen.OR.( upper.AND.i.LE.j ).OR.( lower.AND.i.GE.j ) )
2804 $ THEN
2805 IF( ( i.LE.j.AND.j - i.LE.ku ).OR.
2806 $ ( i.GE.j.AND.i - j.LE.kl ) )THEN
2807 a( i, j ) =
zbeg( reset ) + transl
2808 ELSE
2809 a( i, j ) = zero
2810 END IF
2811 IF( i.NE.j )THEN
2812 IF( sym )THEN
2813 a( j, i ) = dconjg( a( i, j ) )
2814 ELSE IF( tri )THEN
2815 a( j, i ) = zero
2816 END IF
2817 END IF
2818 END IF
2819 10 CONTINUE
2820 IF( sym )
2821 $ a( j, j ) = dcmplx( dble( a( j, j ) ), rzero )
2822 IF( tri )
2823 $ a( j, j ) = a( j, j ) + one
2824 IF( unit )
2825 $ a( j, j ) = one
2826 20 CONTINUE
2827
2828
2829
2830 IF( type.EQ.'GE' )THEN
2831 DO 50 j = 1, n
2832 DO 30 i = 1, m
2833 aa( i + ( j - 1 )*lda ) = a( i, j )
2834 30 CONTINUE
2835 DO 40 i = m + 1, lda
2836 aa( i + ( j - 1 )*lda ) = rogue
2837 40 CONTINUE
2838 50 CONTINUE
2839 ELSE IF( type.EQ.'GB' )THEN
2840 DO 90 j = 1, n
2841 DO 60 i1 = 1, ku + 1 - j
2842 aa( i1 + ( j - 1 )*lda ) = rogue
2843 60 CONTINUE
2844 DO 70 i2 = i1, min( kl + ku + 1, ku + 1 + m - j )
2845 aa( i2 + ( j - 1 )*lda ) = a( i2 + j - ku - 1, j )
2846 70 CONTINUE
2847 DO 80 i3 = i2, lda
2848 aa( i3 + ( j - 1 )*lda ) = rogue
2849 80 CONTINUE
2850 90 CONTINUE
2851 ELSE IF( type.EQ.'HE'.OR.type.EQ.'TR' )THEN
2852 DO 130 j = 1, n
2853 IF( upper )THEN
2854 ibeg = 1
2855 IF( unit )THEN
2856 iend = j - 1
2857 ELSE
2858 iend = j
2859 END IF
2860 ELSE
2861 IF( unit )THEN
2862 ibeg = j + 1
2863 ELSE
2864 ibeg = j
2865 END IF
2866 iend = n
2867 END IF
2868 DO 100 i = 1, ibeg - 1
2869 aa( i + ( j - 1 )*lda ) = rogue
2870 100 CONTINUE
2871 DO 110 i = ibeg, iend
2872 aa( i + ( j - 1 )*lda ) = a( i, j )
2873 110 CONTINUE
2874 DO 120 i = iend + 1, lda
2875 aa( i + ( j - 1 )*lda ) = rogue
2876 120 CONTINUE
2877 IF( sym )THEN
2878 jj = j + ( j - 1 )*lda
2879 aa( jj ) = dcmplx( dble( aa( jj ) ), rrogue )
2880 END IF
2881 130 CONTINUE
2882 ELSE IF( type.EQ.'HB'.OR.type.EQ.'TB' )THEN
2883 DO 170 j = 1, n
2884 IF( upper )THEN
2885 kk = kl + 1
2886 ibeg = max( 1, kl + 2 - j )
2887 IF( unit )THEN
2888 iend = kl
2889 ELSE
2890 iend = kl + 1
2891 END IF
2892 ELSE
2893 kk = 1
2894 IF( unit )THEN
2895 ibeg = 2
2896 ELSE
2897 ibeg = 1
2898 END IF
2899 iend = min( kl + 1, 1 + m - j )
2900 END IF
2901 DO 140 i = 1, ibeg - 1
2902 aa( i + ( j - 1 )*lda ) = rogue
2903 140 CONTINUE
2904 DO 150 i = ibeg, iend
2905 aa( i + ( j - 1 )*lda ) = a( i + j - kk, j )
2906 150 CONTINUE
2907 DO 160 i = iend + 1, lda
2908 aa( i + ( j - 1 )*lda ) = rogue
2909 160 CONTINUE
2910 IF( sym )THEN
2911 jj = kk + ( j - 1 )*lda
2912 aa( jj ) = dcmplx( dble( aa( jj ) ), rrogue )
2913 END IF
2914 170 CONTINUE
2915 ELSE IF( type.EQ.'HP'.OR.type.EQ.'TP' )THEN
2916 ioff = 0
2917 DO 190 j = 1, n
2918 IF( upper )THEN
2919 ibeg = 1
2920 iend = j
2921 ELSE
2922 ibeg = j
2923 iend = n
2924 END IF
2925 DO 180 i = ibeg, iend
2926 ioff = ioff + 1
2927 aa( ioff ) = a( i, j )
2928 IF( i.EQ.j )THEN
2929 IF( unit )
2930 $ aa( ioff ) = rogue
2931 IF( sym )
2932 $ aa( ioff ) = dcmplx( dble( aa( ioff ) ), rrogue )
2933 END IF
2934 180 CONTINUE
2935 190 CONTINUE
2936 END IF
2937 RETURN
2938
2939
2940
complex *16 function zbeg(reset)