LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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' or U'*U
69 *> factorization of A.
70 *> Overwritten with the reconstructed matrix, and then with the
71 *> difference L*L' - A (or U'*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' - A) / ( N * norm(A) * EPS )
89 *> If UPLO = 'U', norm(U'*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 *> \date November 2011
101 *
102 *> \ingroup double_lin
103 *
104 * =====================================================================
105  SUBROUTINE dpot01( UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID )
106 *
107 * -- LAPACK test routine (version 3.4.0) --
108 * -- LAPACK is a software package provided by Univ. of Tennessee, --
109 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
110 * November 2011
111 *
112 * .. Scalar Arguments ..
113  CHARACTER uplo
114  INTEGER lda, ldafac, n
115  DOUBLE PRECISION resid
116 * ..
117 * .. Array Arguments ..
118  DOUBLE PRECISION a( lda, * ), afac( ldafac, * ), rwork( * )
119 * ..
120 *
121 * =====================================================================
122 *
123 * .. Parameters ..
124  DOUBLE PRECISION zero, one
125  parameter( zero = 0.0d+0, one = 1.0d+0 )
126 * ..
127 * .. Local Scalars ..
128  INTEGER i, j, k
129  DOUBLE PRECISION anorm, eps, t
130 * ..
131 * .. External Functions ..
132  LOGICAL lsame
133  DOUBLE PRECISION ddot, dlamch, dlansy
134  EXTERNAL lsame, ddot, dlamch, dlansy
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL dscal, dsyr, dtrmv
138 * ..
139 * .. Intrinsic Functions ..
140  INTRINSIC dble
141 * ..
142 * .. Executable Statements ..
143 *
144 * Quick exit if N = 0.
145 *
146  IF( n.LE.0 ) THEN
147  resid = zero
148  return
149  END IF
150 *
151 * Exit with RESID = 1/EPS if ANORM = 0.
152 *
153  eps = dlamch( 'Epsilon' )
154  anorm = dlansy( '1', uplo, n, a, lda, rwork )
155  IF( anorm.LE.zero ) THEN
156  resid = one / eps
157  return
158  END IF
159 *
160 * Compute the product U'*U, overwriting U.
161 *
162  IF( lsame( uplo, 'U' ) ) THEN
163  DO 10 k = n, 1, -1
164 *
165 * Compute the (K,K) element of the result.
166 *
167  t = ddot( k, afac( 1, k ), 1, afac( 1, k ), 1 )
168  afac( k, k ) = t
169 *
170 * Compute the rest of column K.
171 *
172  CALL dtrmv( 'Upper', 'Transpose', 'Non-unit', k-1, afac,
173  $ ldafac, afac( 1, k ), 1 )
174 *
175  10 continue
176 *
177 * Compute the product L*L', overwriting L.
178 *
179  ELSE
180  DO 20 k = n, 1, -1
181 *
182 * Add a multiple of column K of the factor L to each of
183 * columns K+1 through N.
184 *
185  IF( k+1.LE.n )
186  $ CALL dsyr( 'Lower', n-k, one, afac( k+1, k ), 1,
187  $ afac( k+1, k+1 ), ldafac )
188 *
189 * Scale column K by the diagonal element.
190 *
191  t = afac( k, k )
192  CALL dscal( n-k+1, t, afac( k, k ), 1 )
193 *
194  20 continue
195  END IF
196 *
197 * Compute the difference L*L' - A (or U'*U - A).
198 *
199  IF( lsame( uplo, 'U' ) ) THEN
200  DO 40 j = 1, n
201  DO 30 i = 1, j
202  afac( i, j ) = afac( i, j ) - a( i, j )
203  30 continue
204  40 continue
205  ELSE
206  DO 60 j = 1, n
207  DO 50 i = j, n
208  afac( i, j ) = afac( i, j ) - a( i, j )
209  50 continue
210  60 continue
211  END IF
212 *
213 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
214 *
215  resid = dlansy( '1', uplo, n, afac, ldafac, rwork )
216 *
217  resid = ( ( resid / dble( n ) ) / anorm ) / eps
218 *
219  return
220 *
221 * End of DPOT01
222 *
223  END