LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dspt01.f
Go to the documentation of this file.
1*> \brief \b DSPT01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
12*
13* .. Scalar Arguments ..
14* CHARACTER UPLO
15* INTEGER LDC, N
16* DOUBLE PRECISION RESID
17* ..
18* .. Array Arguments ..
19* INTEGER IPIV( * )
20* DOUBLE PRECISION A( * ), AFAC( * ), C( LDC, * ), RWORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DSPT01 reconstructs a symmetric indefinite packed matrix A from its
30*> block L*D*L' or U*D*U' factorization and computes the residual
31*> norm( C - A ) / ( N * norm(A) * EPS ),
32*> where C is the reconstructed matrix and EPS is the machine epsilon.
33*> \endverbatim
34*
35* Arguments:
36* ==========
37*
38*> \param[in] UPLO
39*> \verbatim
40*> UPLO is CHARACTER*1
41*> Specifies whether the upper or lower triangular part of the
42*> symmetric matrix A is stored:
43*> = 'U': Upper triangular
44*> = 'L': Lower triangular
45*> \endverbatim
46*>
47*> \param[in] N
48*> \verbatim
49*> N is INTEGER
50*> The number of rows and columns of the matrix A. N >= 0.
51*> \endverbatim
52*>
53*> \param[in] A
54*> \verbatim
55*> A is DOUBLE PRECISION array, dimension (N*(N+1)/2)
56*> The original symmetric matrix A, stored as a packed
57*> triangular matrix.
58*> \endverbatim
59*>
60*> \param[in] AFAC
61*> \verbatim
62*> AFAC is DOUBLE PRECISION array, dimension (N*(N+1)/2)
63*> The factored form of the matrix A, stored as a packed
64*> triangular matrix. AFAC contains the block diagonal matrix D
65*> and the multipliers used to obtain the factor L or U from the
66*> block L*D*L' or U*D*U' factorization as computed by DSPTRF.
67*> \endverbatim
68*>
69*> \param[in] IPIV
70*> \verbatim
71*> IPIV is INTEGER array, dimension (N)
72*> The pivot indices from DSPTRF.
73*> \endverbatim
74*>
75*> \param[out] C
76*> \verbatim
77*> C is DOUBLE PRECISION array, dimension (LDC,N)
78*> \endverbatim
79*>
80*> \param[in] LDC
81*> \verbatim
82*> LDC is INTEGER
83*> The leading dimension of the array C. LDC >= max(1,N).
84*> \endverbatim
85*>
86*> \param[out] RWORK
87*> \verbatim
88*> RWORK is DOUBLE PRECISION array, dimension (N)
89*> \endverbatim
90*>
91*> \param[out] RESID
92*> \verbatim
93*> RESID is DOUBLE PRECISION
94*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
95*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
96*> \endverbatim
97*
98* Authors:
99* ========
100*
101*> \author Univ. of Tennessee
102*> \author Univ. of California Berkeley
103*> \author Univ. of Colorado Denver
104*> \author NAG Ltd.
105*
106*> \ingroup double_lin
107*
108* =====================================================================
109 SUBROUTINE dspt01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID )
110*
111* -- LAPACK test routine --
112* -- LAPACK is a software package provided by Univ. of Tennessee, --
113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114*
115* .. Scalar Arguments ..
116 CHARACTER UPLO
117 INTEGER LDC, N
118 DOUBLE PRECISION RESID
119* ..
120* .. Array Arguments ..
121 INTEGER IPIV( * )
122 DOUBLE PRECISION A( * ), AFAC( * ), C( LDC, * ), RWORK( * )
123* ..
124*
125* =====================================================================
126*
127* .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE
129 parameter( zero = 0.0d+0, one = 1.0d+0 )
130* ..
131* .. Local Scalars ..
132 INTEGER I, INFO, J, JC
133 DOUBLE PRECISION ANORM, EPS
134* ..
135* .. External Functions ..
136 LOGICAL LSAME
137 DOUBLE PRECISION DLAMCH, DLANSP, DLANSY
138 EXTERNAL lsame, dlamch, dlansp, dlansy
139* ..
140* .. External Subroutines ..
141 EXTERNAL dlaset, dlavsp
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC dble
145* ..
146* .. Executable Statements ..
147*
148* Quick exit if N = 0.
149*
150 IF( n.LE.0 ) THEN
151 resid = zero
152 RETURN
153 END IF
154*
155* Determine EPS and the norm of A.
156*
157 eps = dlamch( 'Epsilon' )
158 anorm = dlansp( '1', uplo, n, a, rwork )
159*
160* Initialize C to the identity matrix.
161*
162 CALL dlaset( 'Full', n, n, zero, one, c, ldc )
163*
164* Call DLAVSP to form the product D * U' (or D * L' ).
165*
166 CALL dlavsp( uplo, 'Transpose', 'Non-unit', n, n, afac, ipiv, c,
167 $ ldc, info )
168*
169* Call DLAVSP again to multiply by U ( or L ).
170*
171 CALL dlavsp( uplo, 'No transpose', 'Unit', n, n, afac, ipiv, c,
172 $ ldc, info )
173*
174* Compute the difference C - A .
175*
176 IF( lsame( uplo, 'U' ) ) THEN
177 jc = 0
178 DO 20 j = 1, n
179 DO 10 i = 1, j
180 c( i, j ) = c( i, j ) - a( jc+i )
181 10 CONTINUE
182 jc = jc + j
183 20 CONTINUE
184 ELSE
185 jc = 1
186 DO 40 j = 1, n
187 DO 30 i = j, n
188 c( i, j ) = c( i, j ) - a( jc+i-j )
189 30 CONTINUE
190 jc = jc + n - j + 1
191 40 CONTINUE
192 END IF
193*
194* Compute norm( C - A ) / ( N * norm(A) * EPS )
195*
196 resid = dlansy( '1', uplo, n, c, ldc, rwork )
197*
198 IF( anorm.LE.zero ) THEN
199 IF( resid.NE.zero )
200 $ resid = one / eps
201 ELSE
202 resid = ( ( resid / dble( n ) ) / anorm ) / eps
203 END IF
204*
205 RETURN
206*
207* End of DSPT01
208*
209 END
subroutine dlavsp(uplo, trans, diag, n, nrhs, a, ipiv, b, ldb, info)
DLAVSP
Definition dlavsp.f:130
subroutine dspt01(uplo, n, a, afac, ipiv, c, ldc, rwork, resid)
DSPT01
Definition dspt01.f:110
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110