001:       SUBROUTINE SPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR,
002:      $                   BERR, WORK, IWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH.
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          UPLO
012:       INTEGER            INFO, LDB, LDX, N, NRHS
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IWORK( * )
016:       REAL               AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
017:      $                   FERR( * ), WORK( * ), X( LDX, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  SPPRFS improves the computed solution to a system of linear
024: *  equations when the coefficient matrix is symmetric positive definite
025: *  and packed, and provides error bounds and backward error estimates
026: *  for the solution.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  UPLO    (input) CHARACTER*1
032: *          = 'U':  Upper triangle of A is stored;
033: *          = 'L':  Lower triangle of A is stored.
034: *
035: *  N       (input) INTEGER
036: *          The order of the matrix A.  N >= 0.
037: *
038: *  NRHS    (input) INTEGER
039: *          The number of right hand sides, i.e., the number of columns
040: *          of the matrices B and X.  NRHS >= 0.
041: *
042: *  AP      (input) REAL array, dimension (N*(N+1)/2)
043: *          The upper or lower triangle of the symmetric matrix A, packed
044: *          columnwise in a linear array.  The j-th column of A is stored
045: *          in the array AP as follows:
046: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
047: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
048: *
049: *  AFP     (input) REAL array, dimension (N*(N+1)/2)
050: *          The triangular factor U or L from the Cholesky factorization
051: *          A = U**T*U or A = L*L**T, as computed by SPPTRF/CPPTRF,
052: *          packed columnwise in a linear array in the same format as A
053: *          (see AP).
054: *
055: *  B       (input) REAL array, dimension (LDB,NRHS)
056: *          The right hand side matrix B.
057: *
058: *  LDB     (input) INTEGER
059: *          The leading dimension of the array B.  LDB >= max(1,N).
060: *
061: *  X       (input/output) REAL array, dimension (LDX,NRHS)
062: *          On entry, the solution matrix X, as computed by SPPTRS.
063: *          On exit, the improved solution matrix X.
064: *
065: *  LDX     (input) INTEGER
066: *          The leading dimension of the array X.  LDX >= max(1,N).
067: *
068: *  FERR    (output) REAL array, dimension (NRHS)
069: *          The estimated forward error bound for each solution vector
070: *          X(j) (the j-th column of the solution matrix X).
071: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
072: *          is an estimated upper bound for the magnitude of the largest
073: *          element in (X(j) - XTRUE) divided by the magnitude of the
074: *          largest element in X(j).  The estimate is as reliable as
075: *          the estimate for RCOND, and is almost always a slight
076: *          overestimate of the true error.
077: *
078: *  BERR    (output) REAL array, dimension (NRHS)
079: *          The componentwise relative backward error of each solution
080: *          vector X(j) (i.e., the smallest relative change in
081: *          any element of A or B that makes X(j) an exact solution).
082: *
083: *  WORK    (workspace) REAL array, dimension (3*N)
084: *
085: *  IWORK   (workspace) INTEGER array, dimension (N)
086: *
087: *  INFO    (output) INTEGER
088: *          = 0:  successful exit
089: *          < 0:  if INFO = -i, the i-th argument had an illegal value
090: *
091: *  Internal Parameters
092: *  ===================
093: *
094: *  ITMAX is the maximum number of steps of iterative refinement.
095: *
096: *  =====================================================================
097: *
098: *     .. Parameters ..
099:       INTEGER            ITMAX
100:       PARAMETER          ( ITMAX = 5 )
101:       REAL               ZERO
102:       PARAMETER          ( ZERO = 0.0E+0 )
103:       REAL               ONE
104:       PARAMETER          ( ONE = 1.0E+0 )
105:       REAL               TWO
106:       PARAMETER          ( TWO = 2.0E+0 )
107:       REAL               THREE
108:       PARAMETER          ( THREE = 3.0E+0 )
109: *     ..
110: *     .. Local Scalars ..
111:       LOGICAL            UPPER
112:       INTEGER            COUNT, I, IK, J, K, KASE, KK, NZ
113:       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
114: *     ..
115: *     .. Local Arrays ..
116:       INTEGER            ISAVE( 3 )
117: *     ..
118: *     .. External Subroutines ..
119:       EXTERNAL           SAXPY, SCOPY, SLACN2, SPPTRS, SSPMV, XERBLA
120: *     ..
121: *     .. Intrinsic Functions ..
122:       INTRINSIC          ABS, MAX
123: *     ..
124: *     .. External Functions ..
125:       LOGICAL            LSAME
126:       REAL               SLAMCH
127:       EXTERNAL           LSAME, SLAMCH
128: *     ..
129: *     .. Executable Statements ..
130: *
131: *     Test the input parameters.
132: *
133:       INFO = 0
134:       UPPER = LSAME( UPLO, 'U' )
135:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
136:          INFO = -1
137:       ELSE IF( N.LT.0 ) THEN
138:          INFO = -2
139:       ELSE IF( NRHS.LT.0 ) THEN
140:          INFO = -3
141:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
142:          INFO = -7
143:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
144:          INFO = -9
145:       END IF
146:       IF( INFO.NE.0 ) THEN
147:          CALL XERBLA( 'SPPRFS', -INFO )
148:          RETURN
149:       END IF
150: *
151: *     Quick return if possible
152: *
153:       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
154:          DO 10 J = 1, NRHS
155:             FERR( J ) = ZERO
156:             BERR( J ) = ZERO
157:    10    CONTINUE
158:          RETURN
159:       END IF
160: *
161: *     NZ = maximum number of nonzero elements in each row of A, plus 1
162: *
163:       NZ = N + 1
164:       EPS = SLAMCH( 'Epsilon' )
165:       SAFMIN = SLAMCH( 'Safe minimum' )
166:       SAFE1 = NZ*SAFMIN
167:       SAFE2 = SAFE1 / EPS
168: *
169: *     Do for each right hand side
170: *
171:       DO 140 J = 1, NRHS
172: *
173:          COUNT = 1
174:          LSTRES = THREE
175:    20    CONTINUE
176: *
177: *        Loop until stopping criterion is satisfied.
178: *
179: *        Compute residual R = B - A * X
180: *
181:          CALL SCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
182:          CALL SSPMV( UPLO, N, -ONE, AP, X( 1, J ), 1, ONE, WORK( N+1 ),
183:      $               1 )
184: *
185: *        Compute componentwise relative backward error from formula
186: *
187: *        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
188: *
189: *        where abs(Z) is the componentwise absolute value of the matrix
190: *        or vector Z.  If the i-th component of the denominator is less
191: *        than SAFE2, then SAFE1 is added to the i-th components of the
192: *        numerator and denominator before dividing.
193: *
194:          DO 30 I = 1, N
195:             WORK( I ) = ABS( B( I, J ) )
196:    30    CONTINUE
197: *
198: *        Compute abs(A)*abs(X) + abs(B).
199: *
200:          KK = 1
201:          IF( UPPER ) THEN
202:             DO 50 K = 1, N
203:                S = ZERO
204:                XK = ABS( X( K, J ) )
205:                IK = KK
206:                DO 40 I = 1, K - 1
207:                   WORK( I ) = WORK( I ) + ABS( AP( IK ) )*XK
208:                   S = S + ABS( AP( IK ) )*ABS( X( I, J ) )
209:                   IK = IK + 1
210:    40          CONTINUE
211:                WORK( K ) = WORK( K ) + ABS( AP( KK+K-1 ) )*XK + S
212:                KK = KK + K
213:    50       CONTINUE
214:          ELSE
215:             DO 70 K = 1, N
216:                S = ZERO
217:                XK = ABS( X( K, J ) )
218:                WORK( K ) = WORK( K ) + ABS( AP( KK ) )*XK
219:                IK = KK + 1
220:                DO 60 I = K + 1, N
221:                   WORK( I ) = WORK( I ) + ABS( AP( IK ) )*XK
222:                   S = S + ABS( AP( IK ) )*ABS( X( I, J ) )
223:                   IK = IK + 1
224:    60          CONTINUE
225:                WORK( K ) = WORK( K ) + S
226:                KK = KK + ( N-K+1 )
227:    70       CONTINUE
228:          END IF
229:          S = ZERO
230:          DO 80 I = 1, N
231:             IF( WORK( I ).GT.SAFE2 ) THEN
232:                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
233:             ELSE
234:                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
235:      $             ( WORK( I )+SAFE1 ) )
236:             END IF
237:    80    CONTINUE
238:          BERR( J ) = S
239: *
240: *        Test stopping criterion. Continue iterating if
241: *           1) The residual BERR(J) is larger than machine epsilon, and
242: *           2) BERR(J) decreased by at least a factor of 2 during the
243: *              last iteration, and
244: *           3) At most ITMAX iterations tried.
245: *
246:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
247:      $       COUNT.LE.ITMAX ) THEN
248: *
249: *           Update solution and try again.
250: *
251:             CALL SPPTRS( UPLO, N, 1, AFP, WORK( N+1 ), N, INFO )
252:             CALL SAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
253:             LSTRES = BERR( J )
254:             COUNT = COUNT + 1
255:             GO TO 20
256:          END IF
257: *
258: *        Bound error from formula
259: *
260: *        norm(X - XTRUE) / norm(X) .le. FERR =
261: *        norm( abs(inv(A))*
262: *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
263: *
264: *        where
265: *          norm(Z) is the magnitude of the largest component of Z
266: *          inv(A) is the inverse of A
267: *          abs(Z) is the componentwise absolute value of the matrix or
268: *             vector Z
269: *          NZ is the maximum number of nonzeros in any row of A, plus 1
270: *          EPS is machine epsilon
271: *
272: *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
273: *        is incremented by SAFE1 if the i-th component of
274: *        abs(A)*abs(X) + abs(B) is less than SAFE2.
275: *
276: *        Use SLACN2 to estimate the infinity-norm of the matrix
277: *           inv(A) * diag(W),
278: *        where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
279: *
280:          DO 90 I = 1, N
281:             IF( WORK( I ).GT.SAFE2 ) THEN
282:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
283:             ELSE
284:                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
285:             END IF
286:    90    CONTINUE
287: *
288:          KASE = 0
289:   100    CONTINUE
290:          CALL SLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
291:      $                KASE, ISAVE )
292:          IF( KASE.NE.0 ) THEN
293:             IF( KASE.EQ.1 ) THEN
294: *
295: *              Multiply by diag(W)*inv(A').
296: *
297:                CALL SPPTRS( UPLO, N, 1, AFP, WORK( N+1 ), N, INFO )
298:                DO 110 I = 1, N
299:                   WORK( N+I ) = WORK( I )*WORK( N+I )
300:   110          CONTINUE
301:             ELSE IF( KASE.EQ.2 ) THEN
302: *
303: *              Multiply by inv(A)*diag(W).
304: *
305:                DO 120 I = 1, N
306:                   WORK( N+I ) = WORK( I )*WORK( N+I )
307:   120          CONTINUE
308:                CALL SPPTRS( UPLO, N, 1, AFP, WORK( N+1 ), N, INFO )
309:             END IF
310:             GO TO 100
311:          END IF
312: *
313: *        Normalize error.
314: *
315:          LSTRES = ZERO
316:          DO 130 I = 1, N
317:             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
318:   130    CONTINUE
319:          IF( LSTRES.NE.ZERO )
320:      $      FERR( J ) = FERR( J ) / LSTRES
321: *
322:   140 CONTINUE
323: *
324:       RETURN
325: *
326: *     End of SPPRFS
327: *
328:       END
329: