1       SUBROUTINE SGET07( TRANS, N, NRHS, A, LDA, B, LDB, X, LDX, XACT,
  2      $                   LDXACT, FERR, CHKFERR, BERR, RESLTS )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          TRANS
 10       LOGICAL            CHKFERR
 11       INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
 12 *     ..
 13 *     .. Array Arguments ..
 14       REAL               A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
 15      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  SGET07 tests the error bounds from iterative refinement for the
 22 *  computed solution to a system of equations op(A)*X = B, where A is a
 23 *  general n by n matrix and op(A) = A or A**T, depending on TRANS.
 24 *
 25 *  RESLTS(1) = test of the error bound
 26 *            = norm(X - XACT) / ( norm(X) * FERR )
 27 *
 28 *  A large value is returned if this ratio is not less than one.
 29 *
 30 *  RESLTS(2) = residual from the iterative refinement routine
 31 *            = the maximum of BERR / ( (n+1)*EPS + (*) ), where
 32 *              (*) = (n+1)*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
 33 *
 34 *  Arguments
 35 *  =========
 36 *
 37 *  TRANS   (input) CHARACTER*1
 38 *          Specifies the form of the system of equations.
 39 *          = 'N':  A * X = B     (No transpose)
 40 *          = 'T':  A**T * X = B  (Transpose)
 41 *          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
 42 *
 43 *  N       (input) INTEGER
 44 *          The number of rows of the matrices X and XACT.  N >= 0.
 45 *
 46 *  NRHS    (input) INTEGER
 47 *          The number of columns of the matrices X and XACT.  NRHS >= 0.
 48 *
 49 *  A       (input) REAL array, dimension (LDA,N)
 50 *          The original n by n matrix A.
 51 *
 52 *  LDA     (input) INTEGER
 53 *          The leading dimension of the array A.  LDA >= max(1,N).
 54 *
 55 *  B       (input) REAL array, dimension (LDB,NRHS)
 56 *          The right hand side vectors for the system of linear
 57 *          equations.
 58 *
 59 *  LDB     (input) INTEGER
 60 *          The leading dimension of the array B.  LDB >= max(1,N).
 61 *
 62 *  X       (input) REAL array, dimension (LDX,NRHS)
 63 *          The computed solution vectors.  Each vector is stored as a
 64 *          column of the matrix X.
 65 *
 66 *  LDX     (input) INTEGER
 67 *          The leading dimension of the array X.  LDX >= max(1,N).
 68 *
 69 *  XACT    (input) REAL array, dimension (LDX,NRHS)
 70 *          The exact solution vectors.  Each vector is stored as a
 71 *          column of the matrix XACT.
 72 *
 73 *  LDXACT  (input) INTEGER
 74 *          The leading dimension of the array XACT.  LDXACT >= max(1,N).
 75 *
 76 *  FERR    (input) REAL array, dimension (NRHS)
 77 *          The estimated forward error bounds for each solution vector
 78 *          X.  If XTRUE is the true solution, FERR bounds the magnitude
 79 *          of the largest entry in (X - XTRUE) divided by the magnitude
 80 *          of the largest entry in X.
 81 *
 82 *  CHKFERR (input) LOGICAL
 83 *          Set to .TRUE. to check FERR, .FALSE. not to check FERR.
 84 *          When the test system is ill-conditioned, the "true"
 85 *          solution in XACT may be incorrect.
 86 *
 87 *  BERR    (input) REAL array, dimension (NRHS)
 88 *          The componentwise relative backward error of each solution
 89 *          vector (i.e., the smallest relative change in any entry of A
 90 *          or B that makes X an exact solution).
 91 *
 92 *  RESLTS  (output) REAL array, dimension (2)
 93 *          The maximum over the NRHS solution vectors of the ratios:
 94 *          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
 95 *          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
 96 *
 97 *  =====================================================================
 98 *
 99 *     .. Parameters ..
100       REAL               ZERO, ONE
101       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
102 *     ..
103 *     .. Local Scalars ..
104       LOGICAL            NOTRAN
105       INTEGER            I, IMAX, J, K
106       REAL               AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
107 *     ..
108 *     .. External Functions ..
109       LOGICAL            LSAME
110       INTEGER            ISAMAX
111       REAL               SLAMCH
112       EXTERNAL           LSAME, ISAMAX, SLAMCH
113 *     ..
114 *     .. Intrinsic Functions ..
115       INTRINSIC          ABSMAXMIN
116 *     ..
117 *     .. Executable Statements ..
118 *
119 *     Quick exit if N = 0 or NRHS = 0.
120 *
121       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
122          RESLTS( 1 ) = ZERO
123          RESLTS( 2 ) = ZERO
124          RETURN
125       END IF
126 *
127       EPS = SLAMCH( 'Epsilon' )
128       UNFL = SLAMCH( 'Safe minimum' )
129       OVFL = ONE / UNFL
130       NOTRAN = LSAME( TRANS, 'N' )
131 *
132 *     Test 1:  Compute the maximum of
133 *        norm(X - XACT) / ( norm(X) * FERR )
134 *     over all the vectors X and XACT using the infinity-norm.
135 *
136       ERRBND = ZERO
137       IF( CHKFERR ) THEN
138          DO 30 J = 1, NRHS
139             IMAX = ISAMAX( N, X( 1, J ), 1 )
140             XNORM = MAXABS( X( IMAX, J ) ), UNFL )
141             DIFF = ZERO
142             DO 10 I = 1, N
143                DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
144  10         CONTINUE
145 *
146             IF( XNORM.GT.ONE ) THEN
147                GO TO 20
148             ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
149                GO TO 20
150             ELSE
151                ERRBND = ONE / EPS
152                GO TO 30
153             END IF
154 *
155  20         CONTINUE
156             IF( DIFF / XNORM.LE.FERR( J ) ) THEN
157                ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
158             ELSE
159                ERRBND = ONE / EPS
160             END IF
161  30      CONTINUE
162       END IF
163       RESLTS( 1 ) = ERRBND
164 *
165 *     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
166 *     (*) = (n+1)*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
167 *
168       DO 70 K = 1, NRHS
169          DO 60 I = 1, N
170             TMP = ABS( B( I, K ) )
171             IF( NOTRAN ) THEN
172                DO 40 J = 1, N
173                   TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
174    40          CONTINUE
175             ELSE
176                DO 50 J = 1, N
177                   TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
178    50          CONTINUE
179             END IF
180             IF( I.EQ.1 ) THEN
181                AXBI = TMP
182             ELSE
183                AXBI = MIN( AXBI, TMP )
184             END IF
185    60    CONTINUE
186          TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
187      $         MAX( AXBI, ( N+1 )*UNFL ) )
188          IF( K.EQ.1 ) THEN
189             RESLTS( 2 ) = TMP
190          ELSE
191             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
192          END IF
193    70 CONTINUE
194 *
195       RETURN
196 *
197 *     End of SGET07
198 *
199       END