LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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.```
Date
November 2011

Definition at line 136 of file dpst01.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: