LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ stpt01()

subroutine stpt01 ( character  UPLO,
character  DIAG,
integer  N,
real, dimension( * )  AP,
real, dimension( * )  AINVP,
real  RCOND,
real, dimension( * )  WORK,
real  RESID 
)

STPT01

Purpose:
 STPT01 computes the residual for a triangular matrix A times its
 inverse when A is stored in packed format:
    RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
 where EPS is the machine epsilon.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the matrix A is upper or lower triangular.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]DIAG
          DIAG is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]AP
          AP is REAL array, dimension (N*(N+1)/2)
          The original upper or lower triangular matrix A, packed
          columnwise in a linear array.  The j-th column of A is stored
          in the array AP as follows:
          if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
          if UPLO = 'L',
             AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
[in,out]AINVP
          AINVP is REAL array, dimension (N*(N+1)/2)
          On entry, the (triangular) inverse of the matrix A, packed
          columnwise in a linear array as in AP.
          On exit, the contents of AINVP are destroyed.
[out]RCOND
          RCOND is REAL
          The reciprocal condition number of A, computed as
          1/(norm(A) * norm(AINV)).
[out]WORK
          WORK is REAL array, dimension (N)
[out]RESID
          RESID is REAL
          norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 107 of file stpt01.f.

108 *
109 * -- LAPACK test routine --
110 * -- LAPACK is a software package provided by Univ. of Tennessee, --
111 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
112 *
113 * .. Scalar Arguments ..
114  CHARACTER DIAG, UPLO
115  INTEGER N
116  REAL RCOND, RESID
117 * ..
118 * .. Array Arguments ..
119  REAL AINVP( * ), AP( * ), WORK( * )
120 * ..
121 *
122 * =====================================================================
123 *
124 * .. Parameters ..
125  REAL ZERO, ONE
126  parameter( zero = 0.0e+0, one = 1.0e+0 )
127 * ..
128 * .. Local Scalars ..
129  LOGICAL UNITD
130  INTEGER J, JC
131  REAL AINVNM, ANORM, EPS
132 * ..
133 * .. External Functions ..
134  LOGICAL LSAME
135  REAL SLAMCH, SLANTP
136  EXTERNAL lsame, slamch, slantp
137 * ..
138 * .. External Subroutines ..
139  EXTERNAL stpmv
140 * ..
141 * .. Intrinsic Functions ..
142  INTRINSIC real
143 * ..
144 * .. Executable Statements ..
145 *
146 * Quick exit if N = 0.
147 *
148  IF( n.LE.0 ) THEN
149  rcond = one
150  resid = zero
151  RETURN
152  END IF
153 *
154 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
155 *
156  eps = slamch( 'Epsilon' )
157  anorm = slantp( '1', uplo, diag, n, ap, work )
158  ainvnm = slantp( '1', uplo, diag, n, ainvp, work )
159  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
160  rcond = zero
161  resid = one / eps
162  RETURN
163  END IF
164  rcond = ( one / anorm ) / ainvnm
165 *
166 * Compute A * AINV, overwriting AINV.
167 *
168  unitd = lsame( diag, 'U' )
169  IF( lsame( uplo, 'U' ) ) THEN
170  jc = 1
171  DO 10 j = 1, n
172  IF( unitd )
173  $ ainvp( jc+j-1 ) = one
174 *
175 * Form the j-th column of A*AINV
176 *
177  CALL stpmv( 'Upper', 'No transpose', diag, j, ap,
178  $ ainvp( jc ), 1 )
179 *
180 * Subtract 1 from the diagonal
181 *
182  ainvp( jc+j-1 ) = ainvp( jc+j-1 ) - one
183  jc = jc + j
184  10 CONTINUE
185  ELSE
186  jc = 1
187  DO 20 j = 1, n
188  IF( unitd )
189  $ ainvp( jc ) = one
190 *
191 * Form the j-th column of A*AINV
192 *
193  CALL stpmv( 'Lower', 'No transpose', diag, n-j+1, ap( jc ),
194  $ ainvp( jc ), 1 )
195 *
196 * Subtract 1 from the diagonal
197 *
198  ainvp( jc ) = ainvp( jc ) - one
199  jc = jc + n - j + 1
200  20 CONTINUE
201  END IF
202 *
203 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
204 *
205  resid = slantp( '1', uplo, 'Non-unit', n, ainvp, work )
206 *
207  resid = ( ( resid*rcond ) / real( n ) ) / eps
208 *
209  RETURN
210 *
211 * End of STPT01
212 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slantp(NORM, UPLO, DIAG, N, AP, WORK)
SLANTP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slantp.f:124
subroutine stpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
STPMV
Definition: stpmv.f:142
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: