LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cptrfs()

subroutine cptrfs ( character  uplo,
integer  n,
integer  nrhs,
real, dimension( * )  d,
complex, dimension( * )  e,
real, dimension( * )  df,
complex, dimension( * )  ef,
complex, dimension( ldb, * )  b,
integer  ldb,
complex, dimension( ldx, * )  x,
integer  ldx,
real, dimension( * )  ferr,
real, dimension( * )  berr,
complex, dimension( * )  work,
real, dimension( * )  rwork,
integer  info 
)

CPTRFS

Download CPTRFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CPTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian positive definite
 and tridiagonal, and provides error bounds and backward error
 estimates for the solution.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the superdiagonal or the subdiagonal of the
          tridiagonal matrix A is stored and the form of the
          factorization:
          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
          (The two forms are equivalent if A is real.)
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.
[in]D
          D is REAL array, dimension (N)
          The n real diagonal elements of the tridiagonal matrix A.
[in]E
          E is COMPLEX array, dimension (N-1)
          The (n-1) off-diagonal elements of the tridiagonal matrix A
          (see UPLO).
[in]DF
          DF is REAL array, dimension (N)
          The n diagonal elements of the diagonal matrix D from
          the factorization computed by CPTTRF.
[in]EF
          EF is COMPLEX array, dimension (N-1)
          The (n-1) off-diagonal elements of the unit bidiagonal
          factor U or L from the factorization computed by CPTTRF
          (see UPLO).
[in]B
          B is COMPLEX array, dimension (LDB,NRHS)
          The right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by CPTTRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is REAL array, dimension (NRHS)
          The forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).
[out]BERR
          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is COMPLEX array, dimension (N)
[out]RWORK
          RWORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Internal Parameters:
  ITMAX is the maximum number of steps of iterative refinement.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 181 of file cptrfs.f.

183*
184* -- LAPACK computational routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER UPLO
190 INTEGER INFO, LDB, LDX, N, NRHS
191* ..
192* .. Array Arguments ..
193 REAL BERR( * ), D( * ), DF( * ), FERR( * ),
194 $ RWORK( * )
195 COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
196 $ X( LDX, * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 INTEGER ITMAX
203 parameter( itmax = 5 )
204 REAL ZERO
205 parameter( zero = 0.0e+0 )
206 REAL ONE
207 parameter( one = 1.0e+0 )
208 REAL TWO
209 parameter( two = 2.0e+0 )
210 REAL THREE
211 parameter( three = 3.0e+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL UPPER
215 INTEGER COUNT, I, IX, J, NZ
216 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
217 COMPLEX BI, CX, DX, EX, ZDUM
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 INTEGER ISAMAX
222 REAL SLAMCH
223 EXTERNAL lsame, isamax, slamch
224* ..
225* .. External Subroutines ..
226 EXTERNAL caxpy, cpttrs, xerbla
227* ..
228* .. Intrinsic Functions ..
229 INTRINSIC abs, aimag, cmplx, conjg, max, real
230* ..
231* .. Statement Functions ..
232 REAL CABS1
233* ..
234* .. Statement Function definitions ..
235 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
236* ..
237* .. Executable Statements ..
238*
239* Test the input parameters.
240*
241 info = 0
242 upper = lsame( uplo, 'U' )
243 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
244 info = -1
245 ELSE IF( n.LT.0 ) THEN
246 info = -2
247 ELSE IF( nrhs.LT.0 ) THEN
248 info = -3
249 ELSE IF( ldb.LT.max( 1, n ) ) THEN
250 info = -9
251 ELSE IF( ldx.LT.max( 1, n ) ) THEN
252 info = -11
253 END IF
254 IF( info.NE.0 ) THEN
255 CALL xerbla( 'CPTRFS', -info )
256 RETURN
257 END IF
258*
259* Quick return if possible
260*
261 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
262 DO 10 j = 1, nrhs
263 ferr( j ) = zero
264 berr( j ) = zero
265 10 CONTINUE
266 RETURN
267 END IF
268*
269* NZ = maximum number of nonzero elements in each row of A, plus 1
270*
271 nz = 4
272 eps = slamch( 'Epsilon' )
273 safmin = slamch( 'Safe minimum' )
274 safe1 = nz*safmin
275 safe2 = safe1 / eps
276*
277* Do for each right hand side
278*
279 DO 100 j = 1, nrhs
280*
281 count = 1
282 lstres = three
283 20 CONTINUE
284*
285* Loop until stopping criterion is satisfied.
286*
287* Compute residual R = B - A * X. Also compute
288* abs(A)*abs(x) + abs(b) for use in the backward error bound.
289*
290 IF( upper ) THEN
291 IF( n.EQ.1 ) THEN
292 bi = b( 1, j )
293 dx = d( 1 )*x( 1, j )
294 work( 1 ) = bi - dx
295 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
296 ELSE
297 bi = b( 1, j )
298 dx = d( 1 )*x( 1, j )
299 ex = e( 1 )*x( 2, j )
300 work( 1 ) = bi - dx - ex
301 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
302 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
303 DO 30 i = 2, n - 1
304 bi = b( i, j )
305 cx = conjg( e( i-1 ) )*x( i-1, j )
306 dx = d( i )*x( i, j )
307 ex = e( i )*x( i+1, j )
308 work( i ) = bi - cx - dx - ex
309 rwork( i ) = cabs1( bi ) +
310 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
311 $ cabs1( dx ) + cabs1( e( i ) )*
312 $ cabs1( x( i+1, j ) )
313 30 CONTINUE
314 bi = b( n, j )
315 cx = conjg( e( n-1 ) )*x( n-1, j )
316 dx = d( n )*x( n, j )
317 work( n ) = bi - cx - dx
318 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
319 $ cabs1( x( n-1, j ) ) + cabs1( dx )
320 END IF
321 ELSE
322 IF( n.EQ.1 ) THEN
323 bi = b( 1, j )
324 dx = d( 1 )*x( 1, j )
325 work( 1 ) = bi - dx
326 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
327 ELSE
328 bi = b( 1, j )
329 dx = d( 1 )*x( 1, j )
330 ex = conjg( e( 1 ) )*x( 2, j )
331 work( 1 ) = bi - dx - ex
332 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
333 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
334 DO 40 i = 2, n - 1
335 bi = b( i, j )
336 cx = e( i-1 )*x( i-1, j )
337 dx = d( i )*x( i, j )
338 ex = conjg( e( i ) )*x( i+1, j )
339 work( i ) = bi - cx - dx - ex
340 rwork( i ) = cabs1( bi ) +
341 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
342 $ cabs1( dx ) + cabs1( e( i ) )*
343 $ cabs1( x( i+1, j ) )
344 40 CONTINUE
345 bi = b( n, j )
346 cx = e( n-1 )*x( n-1, j )
347 dx = d( n )*x( n, j )
348 work( n ) = bi - cx - dx
349 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
350 $ cabs1( x( n-1, j ) ) + cabs1( dx )
351 END IF
352 END IF
353*
354* Compute componentwise relative backward error from formula
355*
356* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
357*
358* where abs(Z) is the componentwise absolute value of the matrix
359* or vector Z. If the i-th component of the denominator is less
360* than SAFE2, then SAFE1 is added to the i-th components of the
361* numerator and denominator before dividing.
362*
363 s = zero
364 DO 50 i = 1, n
365 IF( rwork( i ).GT.safe2 ) THEN
366 s = max( s, cabs1( work( i ) ) / rwork( i ) )
367 ELSE
368 s = max( s, ( cabs1( work( i ) )+safe1 ) /
369 $ ( rwork( i )+safe1 ) )
370 END IF
371 50 CONTINUE
372 berr( j ) = s
373*
374* Test stopping criterion. Continue iterating if
375* 1) The residual BERR(J) is larger than machine epsilon, and
376* 2) BERR(J) decreased by at least a factor of 2 during the
377* last iteration, and
378* 3) At most ITMAX iterations tried.
379*
380 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
381 $ count.LE.itmax ) THEN
382*
383* Update solution and try again.
384*
385 CALL cpttrs( uplo, n, 1, df, ef, work, n, info )
386 CALL caxpy( n, cmplx( one ), work, 1, x( 1, j ), 1 )
387 lstres = berr( j )
388 count = count + 1
389 GO TO 20
390 END IF
391*
392* Bound error from formula
393*
394* norm(X - XTRUE) / norm(X) .le. FERR =
395* norm( abs(inv(A))*
396* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
397*
398* where
399* norm(Z) is the magnitude of the largest component of Z
400* inv(A) is the inverse of A
401* abs(Z) is the componentwise absolute value of the matrix or
402* vector Z
403* NZ is the maximum number of nonzeros in any row of A, plus 1
404* EPS is machine epsilon
405*
406* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
407* is incremented by SAFE1 if the i-th component of
408* abs(A)*abs(X) + abs(B) is less than SAFE2.
409*
410 DO 60 i = 1, n
411 IF( rwork( i ).GT.safe2 ) THEN
412 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
413 ELSE
414 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
415 $ safe1
416 END IF
417 60 CONTINUE
418 ix = isamax( n, rwork, 1 )
419 ferr( j ) = rwork( ix )
420*
421* Estimate the norm of inv(A).
422*
423* Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
424*
425* m(i,j) = abs(A(i,j)), i = j,
426* m(i,j) = -abs(A(i,j)), i .ne. j,
427*
428* and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
429*
430* Solve M(L) * x = e.
431*
432 rwork( 1 ) = one
433 DO 70 i = 2, n
434 rwork( i ) = one + rwork( i-1 )*abs( ef( i-1 ) )
435 70 CONTINUE
436*
437* Solve D * M(L)**H * x = b.
438*
439 rwork( n ) = rwork( n ) / df( n )
440 DO 80 i = n - 1, 1, -1
441 rwork( i ) = rwork( i ) / df( i ) +
442 $ rwork( i+1 )*abs( ef( i ) )
443 80 CONTINUE
444*
445* Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
446*
447 ix = isamax( n, rwork, 1 )
448 ferr( j ) = ferr( j )*abs( rwork( ix ) )
449*
450* Normalize error.
451*
452 lstres = zero
453 DO 90 i = 1, n
454 lstres = max( lstres, abs( x( i, j ) ) )
455 90 CONTINUE
456 IF( lstres.NE.zero )
457 $ ferr( j ) = ferr( j ) / lstres
458*
459 100 CONTINUE
460*
461 RETURN
462*
463* End of CPTRFS
464*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cpttrs
subroutine cpttrs(uplo, n, nrhs, d, e, b, ldb, info)
CPTTRS
Definition cpttrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: