 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ 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.```

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:132
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:147
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
Here is the call graph for this function:
Here is the caller graph for this function: