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

◆ dpst01()

subroutine dpst01 ( character  uplo,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldafac, * )  afac,
integer  ldafac,
double precision, dimension( ldperm, * )  perm,
integer  ldperm,
integer, dimension( * )  piv,
double precision, dimension( * )  rwork,
double precision  resid,
integer  rank 
)

DPST01

Purpose:
 DPST01 reconstructs a symmetric positive semidefinite matrix A
 from its L or U factors and the permutation matrix P and computes
 the residual
    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
 where EPS is the machine epsilon.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          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 DOUBLE PRECISION array, dimension (LDA,N)
          The original symmetric matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N)
[in]AFAC
          AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
          The factor L or U from the L*L' or U'*U
          factorization of A.
[in]LDAFAC
          LDAFAC is INTEGER
          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
[out]PERM
          PERM is DOUBLE PRECISION array, dimension (LDPERM,N)
          Overwritten with the reconstructed matrix, and then with the
          difference P*L*L'*P' - A (or P*U'*U*P' - A)
[in]LDPERM
          LDPERM is INTEGER
          The leading dimension of the array PERM.
          LDAPERM >= max(1,N).
[in]PIV
          PIV is INTEGER array, dimension (N)
          PIV is such that the nonzero entries are
          P( PIV( K ), K ) = 1.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]RESID
          RESID is DOUBLE PRECISION
          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
[in]RANK
          RANK is INTEGER
          number of nonzero singular values of A.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 132 of file dpst01.f.

134*
135* -- LAPACK test routine --
136* -- LAPACK is a software package provided by Univ. of Tennessee, --
137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139* .. Scalar Arguments ..
140 DOUBLE PRECISION RESID
141 INTEGER LDA, LDAFAC, LDPERM, N, RANK
142 CHARACTER UPLO
143* ..
144* .. Array Arguments ..
145 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ),
146 $ PERM( LDPERM, * ), RWORK( * )
147 INTEGER PIV( * )
148* ..
149*
150* =====================================================================
151*
152* .. Parameters ..
153 DOUBLE PRECISION ZERO, ONE
154 parameter( zero = 0.0d+0, one = 1.0d+0 )
155* ..
156* .. Local Scalars ..
157 DOUBLE PRECISION ANORM, EPS, T
158 INTEGER I, J, K
159* ..
160* .. External Functions ..
161 DOUBLE PRECISION DDOT, DLAMCH, DLANSY
162 LOGICAL LSAME
163 EXTERNAL ddot, dlamch, dlansy, lsame
164* ..
165* .. External Subroutines ..
166 EXTERNAL dscal, dsyr, dtrmv
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC dble
170* ..
171* .. Executable Statements ..
172*
173* Quick exit if N = 0.
174*
175 IF( n.LE.0 ) THEN
176 resid = zero
177 RETURN
178 END IF
179*
180* Exit with RESID = 1/EPS if ANORM = 0.
181*
182 eps = dlamch( 'Epsilon' )
183 anorm = dlansy( '1', uplo, n, a, lda, rwork )
184 IF( anorm.LE.zero ) THEN
185 resid = one / eps
186 RETURN
187 END IF
188*
189* Compute the product U'*U, overwriting U.
190*
191 IF( lsame( uplo, 'U' ) ) THEN
192*
193 IF( rank.LT.n ) THEN
194 DO 110 j = rank + 1, n
195 DO 100 i = rank + 1, j
196 afac( i, j ) = zero
197 100 CONTINUE
198 110 CONTINUE
199 END IF
200*
201 DO 120 k = n, 1, -1
202*
203* Compute the (K,K) element of the result.
204*
205 t = ddot( k, afac( 1, k ), 1, afac( 1, k ), 1 )
206 afac( k, k ) = t
207*
208* Compute the rest of column K.
209*
210 CALL dtrmv( 'Upper', 'Transpose', 'Non-unit', k-1, afac,
211 $ ldafac, afac( 1, k ), 1 )
212*
213 120 CONTINUE
214*
215* Compute the product L*L', overwriting L.
216*
217 ELSE
218*
219 IF( rank.LT.n ) THEN
220 DO 140 j = rank + 1, n
221 DO 130 i = j, n
222 afac( i, j ) = zero
223 130 CONTINUE
224 140 CONTINUE
225 END IF
226*
227 DO 150 k = n, 1, -1
228* Add a multiple of column K of the factor L to each of
229* columns K+1 through N.
230*
231 IF( k+1.LE.n )
232 $ CALL dsyr( 'Lower', n-k, one, afac( k+1, k ), 1,
233 $ afac( k+1, k+1 ), ldafac )
234*
235* Scale column K by the diagonal element.
236*
237 t = afac( k, k )
238 CALL dscal( n-k+1, t, afac( k, k ), 1 )
239 150 CONTINUE
240*
241 END IF
242*
243* Form P*L*L'*P' or P*U'*U*P'
244*
245 IF( lsame( uplo, 'U' ) ) THEN
246*
247 DO 170 j = 1, n
248 DO 160 i = 1, n
249 IF( piv( i ).LE.piv( j ) ) THEN
250 IF( i.LE.j ) THEN
251 perm( piv( i ), piv( j ) ) = afac( i, j )
252 ELSE
253 perm( piv( i ), piv( j ) ) = afac( j, i )
254 END IF
255 END IF
256 160 CONTINUE
257 170 CONTINUE
258*
259*
260 ELSE
261*
262 DO 190 j = 1, n
263 DO 180 i = 1, n
264 IF( piv( i ).GE.piv( j ) ) THEN
265 IF( i.GE.j ) THEN
266 perm( piv( i ), piv( j ) ) = afac( i, j )
267 ELSE
268 perm( piv( i ), piv( j ) ) = afac( j, i )
269 END IF
270 END IF
271 180 CONTINUE
272 190 CONTINUE
273*
274 END IF
275*
276* Compute the difference P*L*L'*P' - A (or P*U'*U*P' - A).
277*
278 IF( lsame( uplo, 'U' ) ) THEN
279 DO 210 j = 1, n
280 DO 200 i = 1, j
281 perm( i, j ) = perm( i, j ) - a( i, j )
282 200 CONTINUE
283 210 CONTINUE
284 ELSE
285 DO 230 j = 1, n
286 DO 220 i = j, n
287 perm( i, j ) = perm( i, j ) - a( i, j )
288 220 CONTINUE
289 230 CONTINUE
290 END IF
291*
292* Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
293* ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
294*
295 resid = dlansy( '1', uplo, n, perm, ldafac, rwork )
296*
297 resid = ( ( resid / dble( n ) ) / anorm ) / eps
298*
299 RETURN
300*
301* End of DPST01
302*
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
dsyr
subroutine dsyr(uplo, n, alpha, x, incx, a, lda)
DSYR
Definition dsyr.f:132
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlansy
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansy.f:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtrmv
subroutine dtrmv(uplo, trans, diag, n, a, lda, x, incx)
DTRMV
Definition dtrmv.f:147
Here is the call graph for this function:
Here is the caller graph for this function: