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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 *
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, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
dsyr
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:134
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
dtrmv
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:149
lsame
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: