2835
2836
2837
2838
2839
2840
2841
2842
2843
2844 DOUBLE PRECISION ZERO, ONE
2845 parameter( zero = 0.0d0, one = 1.0d0 )
2846
2847 DOUBLE PRECISION ALPHA, BETA, EPS, ERR
2848 INTEGER KK, LDA, LDB, LDC, LDCC, N, NOUT
2849 LOGICAL FATAL, MV
2850 CHARACTER*1 UPLO, TRANSA, TRANSB
2851
2852 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
2853 $ CC( LDCC, * ), CT( * ), G( * )
2854
2855 DOUBLE PRECISION ERRI
2856 INTEGER I, J, K, ISTART, ISTOP
2857 LOGICAL TRANA, TRANB, UPPER
2858
2859 INTRINSIC abs, max, sqrt
2860
2861 upper = uplo.EQ.'U'
2862 trana = transa.EQ.'T'.OR.transa.EQ.'C'
2863 tranb = transb.EQ.'T'.OR.transb.EQ.'C'
2864
2865
2866
2867
2868
2869 istart = 1
2870 istop = n
2871
2872 DO 120 j = 1, n
2873
2874 IF ( upper ) THEN
2875 istart = 1
2876 istop = j
2877 ELSE
2878 istart = j
2879 istop = n
2880 END IF
2881 DO 10 i = istart, istop
2882 ct( i ) = zero
2883 g( i ) = zero
2884 10 CONTINUE
2885 IF( .NOT.trana.AND..NOT.tranb )THEN
2886 DO 30 k = 1, kk
2887 DO 20 i = istart, istop
2888 ct( i ) = ct( i ) + a( i, k )*b( k, j )
2889 g( i ) = g( i ) + abs( a( i, k ) )*abs( b( k, j ) )
2890 20 CONTINUE
2891 30 CONTINUE
2892 ELSE IF( trana.AND..NOT.tranb )THEN
2893 DO 50 k = 1, kk
2894 DO 40 i = istart, istop
2895 ct( i ) = ct( i ) + a( k, i )*b( k, j )
2896 g( i ) = g( i ) + abs( a( k, i ) )*abs( b( k, j ) )
2897 40 CONTINUE
2898 50 CONTINUE
2899 ELSE IF( .NOT.trana.AND.tranb )THEN
2900 DO 70 k = 1, kk
2901 DO 60 i = istart, istop
2902 ct( i ) = ct( i ) + a( i, k )*b( j, k )
2903 g( i ) = g( i ) + abs( a( i, k ) )*abs( b( j, k ) )
2904 60 CONTINUE
2905 70 CONTINUE
2906 ELSE IF( trana.AND.tranb )THEN
2907 DO 90 k = 1, kk
2908 DO 80 i = istart, istop
2909 ct( i ) = ct( i ) + a( k, i )*b( j, k )
2910 g( i ) = g( i ) + abs( a( k, i ) )*abs( b( j, k ) )
2911 80 CONTINUE
2912 90 CONTINUE
2913 END IF
2914 DO 100 i = istart, istop
2915 ct( i ) = alpha*ct( i ) + beta*c( i, j )
2916 g( i ) = abs( alpha )*g( i ) + abs( beta )*abs( c( i, j ) )
2917 100 CONTINUE
2918
2919
2920
2921 err = zero
2922 DO 110 i = istart, istop
2923 erri = abs( ct( i ) - cc( i, j ) )/eps
2924 IF( g( i ).NE.zero )
2925 $ erri = erri/g( i )
2926 err = max( err, erri )
2927 IF( err*sqrt( eps ).GE.one )
2928 $ GO TO 130
2929 110 CONTINUE
2930
2931 120 CONTINUE
2932
2933
2934 GO TO 150
2935
2936
2937
2938 130 fatal = .true.
2939 WRITE( nout, fmt = 9999 )
2940 DO 140 i = istart, istop
2941 IF( mv )THEN
2942 WRITE( nout, fmt = 9998 )i, ct( i ), cc( i, j )
2943 ELSE
2944 WRITE( nout, fmt = 9998 )i, cc( i, j ), ct( i )
2945 END IF
2946 140 CONTINUE
2947 IF( n.GT.1 )
2948 $ WRITE( nout, fmt = 9997 )j
2949
2950 150 CONTINUE
2951 RETURN
2952
2953 9999 FORMAT( ' ******* FATAL ERROR - COMPUTED RESULT IS LESS THAN HAL',
2954 $ 'F ACCURATE *******', /' EXPECTED RESULT COMPU',
2955 $ 'TED RESULT' )
2956 9998 FORMAT( 1x, i7, 2g18.6 )
2957 9997 FORMAT( ' THESE ARE THE RESULTS FOR COLUMN ', i3 )
2958
2959
2960