cspt03.f

Go to the documentation of this file.

*> \brief \b CSPT03
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
*                          RESID )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            LDW, N
*       REAL               RCOND, RESID
*       ..
*       .. Array Arguments ..
*       REAL               RWORK( * )
*       COMPLEX            A( * ), AINV( * ), WORK( LDW, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          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
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of rows and columns of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is COMPLEX array, dimension (N*(N+1)/2)
*>          The original complex symmetric matrix A, stored as a packed
*>          triangular matrix.
*> \endverbatim
*>
*> \param[in] AINV
*> \verbatim
*>          AINV is COMPLEX array, dimension (N*(N+1)/2)
*>          The (symmetric) inverse of the matrix A, stored as a packed
*>          triangular matrix.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX array, dimension (LDW,N)
*> \endverbatim
*>
*> \param[in] LDW
*> \verbatim
*>          LDW is INTEGER
*>          The leading dimension of the array WORK.  LDW >= max(1,N).
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is REAL array, dimension (N)
*> \endverbatim
*>
*> \param[out] RCOND
*> \verbatim
*>          RCOND is REAL
*>          The reciprocal of the condition number of A, computed as
*>          ( 1/norm(A) ) / norm(AINV).
*> \endverbatim
*>
*> \param[out] RESID
*> \verbatim
*>          RESID is REAL
*>          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex_lin
*
*  =====================================================================
      SUBROUTINE cspt03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
     $                   RESID )
*
*  -- 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
*
      END