cspt03()

subroutine cspt03	(	character	uplo,
		integer	n,
		complex, dimension( * )	a,
		complex, dimension( * )	ainv,
		complex, dimension( ldw, * )	work,
		integer	ldw,
		real, dimension( * )	rwork,
		real	rcond,
		real	resid
	)

CSPT03

Purpose:

 CSPT03 computes the residual for a complex symmetric packed matrix
 times its inverse:
    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
 where EPS is the machine epsilon.

Parameters

[in]	UPLO	UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the complex symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular
[in]	N	N is INTEGER The number of rows and columns of the matrix A. N >= 0.
[in]	A	A is COMPLEX array, dimension (N*(N+1)/2) The original complex symmetric matrix A, stored as a packed triangular matrix.
[in]	AINV	AINV is COMPLEX array, dimension (N*(N+1)/2) The (symmetric) inverse of the matrix A, stored as a packed triangular matrix.
[out]	WORK	WORK is COMPLEX array, dimension (LDW,N)
[in]	LDW	LDW is INTEGER The leading dimension of the array WORK. LDW >= max(1,N).
[out]	RWORK	RWORK is REAL array, dimension (N)
[out]	RCOND	RCOND is REAL The reciprocal of the condition number of A, computed as ( 1/norm(A) ) / norm(AINV).
[out]	RESID	RESID is REAL norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 108 of file cspt03.f.

*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            LDW, N
      REAL               RCOND, RESID
*     ..
*     .. Array Arguments ..
      REAL               RWORK( * )
      COMPLEX            A( * ), AINV( * ), WORK( LDW, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ICOL, J, JCOL, K, KCOL, NALL
      REAL               AINVNM, ANORM, EPS
      COMPLEX            T
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               CLANGE, CLANSP, SLAMCH
      COMPLEX            CDOTU
      EXTERNAL           lsame, clange, clansp, slamch, cdotu
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          real
*     ..
*     .. Executable Statements ..
*
*     Quick exit if N = 0.
*
      IF( n.LE.0 ) THEN
         rcond = one
         resid = zero
         RETURN
      END IF
*
*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
*
      eps = slamch( 'Epsilon' )
      anorm = clansp( '1', uplo, n, a, rwork )
      ainvnm = clansp( '1', uplo, n, ainv, rwork )
      IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
         rcond = zero
         resid = one / eps
         RETURN
      END IF
      rcond = ( one/anorm ) / ainvnm
*
*     Case where both A and AINV are upper triangular:
*     Each element of - A * AINV is computed by taking the dot product
*     of a row of A with a column of AINV.
*
      IF( lsame( uplo, 'U' ) ) THEN
         DO 70 i = 1, n
            icol = ( ( i-1 )*i ) / 2 + 1
*
*           Code when J <= I
*
            DO 30 j = 1, i
               jcol = ( ( j-1 )*j ) / 2 + 1
               t = cdotu( j, a( icol ), 1, ainv( jcol ), 1 )
               jcol = jcol + 2*j - 1
               kcol = icol - 1
               DO 10 k = j + 1, i
                  t = t + a( kcol+k )*ainv( jcol )
                  jcol = jcol + k
   10          CONTINUE
               kcol = kcol + 2*i
               DO 20 k = i + 1, n
                  t = t + a( kcol )*ainv( jcol )
                  kcol = kcol + k
                  jcol = jcol + k
   20          CONTINUE
               work( i, j ) = -t
   30       CONTINUE
*
*           Code when J > I
*
            DO 60 j = i + 1, n
               jcol = ( ( j-1 )*j ) / 2 + 1
               t = cdotu( i, a( icol ), 1, ainv( jcol ), 1 )
               jcol = jcol - 1
               kcol = icol + 2*i - 1
               DO 40 k = i + 1, j
                  t = t + a( kcol )*ainv( jcol+k )
                  kcol = kcol + k
   40          CONTINUE
               jcol = jcol + 2*j
               DO 50 k = j + 1, n
                  t = t + a( kcol )*ainv( jcol )
                  kcol = kcol + k
                  jcol = jcol + k
   50          CONTINUE
               work( i, j ) = -t
   60       CONTINUE
   70    CONTINUE
      ELSE
*
*        Case where both A and AINV are lower triangular
*
         nall = ( n*( n+1 ) ) / 2
         DO 140 i = 1, n
*
*           Code when J <= I
*
            icol = nall - ( ( n-i+1 )*( n-i+2 ) ) / 2 + 1
            DO 100 j = 1, i
               jcol = nall - ( ( n-j )*( n-j+1 ) ) / 2 - ( n-i )
               t = cdotu( n-i+1, a( icol ), 1, ainv( jcol ), 1 )
               kcol = i
               jcol = j
               DO 80 k = 1, j - 1
                  t = t + a( kcol )*ainv( jcol )
                  jcol = jcol + n - k
                  kcol = kcol + n - k
   80          CONTINUE
               jcol = jcol - j
               DO 90 k = j, i - 1
                  t = t + a( kcol )*ainv( jcol+k )
                  kcol = kcol + n - k
   90          CONTINUE
               work( i, j ) = -t
  100       CONTINUE
*
*           Code when J > I
*
            icol = nall - ( ( n-i )*( n-i+1 ) ) / 2
            DO 130 j = i + 1, n
               jcol = nall - ( ( n-j+1 )*( n-j+2 ) ) / 2 + 1
               t = cdotu( n-j+1, a( icol-n+j ), 1, ainv( jcol ), 1 )
               kcol = i
               jcol = j
               DO 110 k = 1, i - 1
                  t = t + a( kcol )*ainv( jcol )
                  jcol = jcol + n - k
                  kcol = kcol + n - k
  110          CONTINUE
               kcol = kcol - i
               DO 120 k = i, j - 1
                  t = t + a( kcol+k )*ainv( jcol )
                  jcol = jcol + n - k
  120          CONTINUE
               work( i, j ) = -t
  130       CONTINUE
  140    CONTINUE
      END IF
*
*     Add the identity matrix to WORK .
*
      DO 150 i = 1, n
         work( i, i ) = work( i, i ) + one
  150 CONTINUE
*
*     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
*
      resid = clange( '1', n, n, work, ldw, rwork )
*
      resid = ( ( resid*rcond )/eps ) / real( n )
*
      RETURN
*
*     End of CSPT03
*

Here is the caller graph for this function: