LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dpot01.f
Go to the documentation of this file.
1 *> \brief \b DPOT01
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 DPOT01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID )
12 *
13 * .. Scalar Arguments ..
14 * CHARACTER UPLO
15 * INTEGER LDA, LDAFAC, N
16 * DOUBLE PRECISION RESID
17 * ..
18 * .. Array Arguments ..
19 * DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> DPOT01 reconstructs a symmetric positive definite matrix A from
29 *> its L*L' or U'*U factorization and computes the residual
30 *> norm( L*L' - A ) / ( N * norm(A) * EPS ) or
31 *> norm( U'*U - A ) / ( N * norm(A) * EPS ),
32 *> where 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 (LDA,N)
56 *> The original symmetric matrix A.
57 *> \endverbatim
58 *>
59 *> \param[in] LDA
60 *> \verbatim
61 *> LDA is INTEGER
62 *> The leading dimension of the array A. LDA >= max(1,N)
63 *> \endverbatim
64 *>
65 *> \param[in,out] AFAC
66 *> \verbatim
67 *> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
68 *> On entry, the factor L or U from the L * L**T or U**T * U
69 *> factorization of A.
70 *> Overwritten with the reconstructed matrix, and then with
71 *> the difference L * L**T - A (or U**T * U - A).
72 *> \endverbatim
73 *>
74 *> \param[in] LDAFAC
75 *> \verbatim
76 *> LDAFAC is INTEGER
77 *> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
78 *> \endverbatim
79 *>
80 *> \param[out] RWORK
81 *> \verbatim
82 *> RWORK is DOUBLE PRECISION array, dimension (N)
83 *> \endverbatim
84 *>
85 *> \param[out] RESID
86 *> \verbatim
87 *> RESID is DOUBLE PRECISION
88 *> If UPLO = 'L', norm(L * L**T - A) / ( N * norm(A) * EPS )
89 *> If UPLO = 'U', norm(U**T * U - A) / ( N * norm(A) * EPS )
90 *> \endverbatim
91 *
92 * Authors:
93 * ========
94 *
95 *> \author Univ. of Tennessee
96 *> \author Univ. of California Berkeley
97 *> \author Univ. of Colorado Denver
98 *> \author NAG Ltd.
99 *
100 *> \ingroup double_lin
101 *
102 * =====================================================================
103  SUBROUTINE dpot01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID )
104 *
105 * -- LAPACK test routine --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 *
109 * .. Scalar Arguments ..
110  CHARACTER UPLO
111  INTEGER LDA, LDAFAC, N
112  DOUBLE PRECISION RESID
113 * ..
114 * .. Array Arguments ..
115  DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
116 * ..
117 *
118 * =====================================================================
119 *
120 * .. Parameters ..
121  DOUBLE PRECISION ZERO, ONE
122  parameter( zero = 0.0d+0, one = 1.0d+0 )
123 * ..
124 * .. Local Scalars ..
125  INTEGER I, J, K
126  DOUBLE PRECISION ANORM, EPS, T
127 * ..
128 * .. External Functions ..
129  LOGICAL LSAME
130  DOUBLE PRECISION DDOT, DLAMCH, DLANSY
131  EXTERNAL lsame, ddot, dlamch, dlansy
132 * ..
133 * .. External Subroutines ..
134  EXTERNAL dscal, dsyr, dtrmv
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC dble
138 * ..
139 * .. Executable Statements ..
140 *
141 * Quick exit if N = 0.
142 *
143  IF( n.LE.0 ) THEN
144  resid = zero
145  RETURN
146  END IF
147 *
148 * Exit with RESID = 1/EPS if ANORM = 0.
149 *
150  eps = dlamch( 'Epsilon' )
151  anorm = dlansy( '1', uplo, n, a, lda, rwork )
152  IF( anorm.LE.zero ) THEN
153  resid = one / eps
154  RETURN
155  END IF
156 *
157 * Compute the product U**T * U, overwriting U.
158 *
159  IF( lsame( uplo, 'U' ) ) THEN
160  DO 10 k = n, 1, -1
161 *
162 * Compute the (K,K) element of the result.
163 *
164  t = ddot( k, afac( 1, k ), 1, afac( 1, k ), 1 )
165  afac( k, k ) = t
166 *
167 * Compute the rest of column K.
168 *
169  CALL dtrmv( 'Upper', 'Transpose', 'Non-unit', k-1, afac,
170  $ ldafac, afac( 1, k ), 1 )
171 *
172  10 CONTINUE
173 *
174 * Compute the product L * L**T, overwriting L.
175 *
176  ELSE
177  DO 20 k = n, 1, -1
178 *
179 * Add a multiple of column K of the factor L to each of
180 * columns K+1 through N.
181 *
182  IF( k+1.LE.n )
183  $ CALL dsyr( 'Lower', n-k, one, afac( k+1, k ), 1,
184  $ afac( k+1, k+1 ), ldafac )
185 *
186 * Scale column K by the diagonal element.
187 *
188  t = afac( k, k )
189  CALL dscal( n-k+1, t, afac( k, k ), 1 )
190 *
191  20 CONTINUE
192  END IF
193 *
194 * Compute the difference L * L**T - A (or U**T * U - A).
195 *
196  IF( lsame( uplo, 'U' ) ) THEN
197  DO 40 j = 1, n
198  DO 30 i = 1, j
199  afac( i, j ) = afac( i, j ) - a( i, j )
200  30 CONTINUE
201  40 CONTINUE
202  ELSE
203  DO 60 j = 1, n
204  DO 50 i = j, n
205  afac( i, j ) = afac( i, j ) - a( i, j )
206  50 CONTINUE
207  60 CONTINUE
208  END IF
209 *
210 * Compute norm(L*U - A) / ( N * norm(A) * EPS )
211 *
212  resid = dlansy( '1', uplo, n, afac, ldafac, rwork )
213 *
214  resid = ( ( resid / dble( n ) ) / anorm ) / eps
215 *
216  RETURN
217 *
218 * End of DPOT01
219 *
220  END
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
subroutine dpot01(UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID)
DPOT01
Definition: dpot01.f:104