1*> \brief \b CPST01
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 CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
12*                          PIV, RWORK, RESID, RANK )
13*
14*       .. Scalar Arguments ..
15*       REAL               RESID
16*       INTEGER            LDA, LDAFAC, LDPERM, N, RANK
17*       CHARACTER          UPLO
18*       ..
19*       .. Array Arguments ..
20*       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ),
21*      $                   PERM( LDPERM, * )
22*       REAL               RWORK( * )
23*       INTEGER            PIV( * )
24*       ..
25*
26*
27*> \par Purpose:
28*  =============
29*>
30*> \verbatim
31*>
32*> CPST01 reconstructs an Hermitian positive semidefinite matrix A
33*> from its L or U factors and the permutation matrix P and computes
34*> the residual
35*>    norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or
36*>    norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ),
37*> where EPS is the machine epsilon, L' is the conjugate transpose of L,
38*> and U' is the conjugate transpose of U.
39*> \endverbatim
40*
41*  Arguments:
42*  ==========
43*
44*> \param[in] UPLO
45*> \verbatim
46*>          UPLO is CHARACTER*1
47*>          Specifies whether the upper or lower triangular part of the
48*>          Hermitian matrix A is stored:
49*>          = 'U':  Upper triangular
50*>          = 'L':  Lower triangular
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*>          N is INTEGER
56*>          The number of rows and columns of the matrix A.  N >= 0.
57*> \endverbatim
58*>
59*> \param[in] A
60*> \verbatim
61*>          A is COMPLEX array, dimension (LDA,N)
62*>          The original Hermitian matrix A.
63*> \endverbatim
64*>
65*> \param[in] LDA
66*> \verbatim
67*>          LDA is INTEGER
68*>          The leading dimension of the array A.  LDA >= max(1,N)
69*> \endverbatim
70*>
71*> \param[in] AFAC
72*> \verbatim
73*>          AFAC is COMPLEX array, dimension (LDAFAC,N)
74*>          The factor L or U from the L*L' or U'*U
75*>          factorization of A.
76*> \endverbatim
77*>
78*> \param[in] LDAFAC
79*> \verbatim
80*>          LDAFAC is INTEGER
81*>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
82*> \endverbatim
83*>
84*> \param[out] PERM
85*> \verbatim
86*>          PERM is COMPLEX array, dimension (LDPERM,N)
87*>          Overwritten with the reconstructed matrix, and then with the
88*>          difference P*L*L'*P' - A (or P*U'*U*P' - A)
89*> \endverbatim
90*>
91*> \param[in] LDPERM
92*> \verbatim
93*>          LDPERM is INTEGER
94*>          The leading dimension of the array PERM.
95*>          LDAPERM >= max(1,N).
96*> \endverbatim
97*>
98*> \param[in] PIV
99*> \verbatim
100*>          PIV is INTEGER array, dimension (N)
101*>          PIV is such that the nonzero entries are
102*>          P( PIV( K ), K ) = 1.
103*> \endverbatim
104*>
105*> \param[out] RWORK
106*> \verbatim
107*>          RWORK is REAL array, dimension (N)
108*> \endverbatim
109*>
110*> \param[out] RESID
111*> \verbatim
112*>          RESID is REAL
113*>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
114*>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
115*> \endverbatim
116*>
117*> \param[in] RANK
118*> \verbatim
119*>          RANK is INTEGER
120*>          number of nonzero singular values of A.
121*> \endverbatim
122*
123*  Authors:
124*  ========
125*
126*> \author Univ. of Tennessee
127*> \author Univ. of California Berkeley
128*> \author Univ. of Colorado Denver
129*> \author NAG Ltd.
130*
131*> \ingroup complex_lin
132*
133*  =====================================================================
134      SUBROUTINE CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
135     $                   PIV, RWORK, RESID, RANK )
136*
137*  -- LAPACK test routine --
138*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
139*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*
141*     .. Scalar Arguments ..
142      REAL               RESID
143      INTEGER            LDA, LDAFAC, LDPERM, N, RANK
144      CHARACTER          UPLO
145*     ..
146*     .. Array Arguments ..
147      COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ),
148     $                   PERM( LDPERM, * )
149      REAL               RWORK( * )
150      INTEGER            PIV( * )
151*     ..
152*
153*  =====================================================================
154*
155*     .. Parameters ..
156      REAL               ZERO, ONE
157      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
158      COMPLEX            CZERO
159      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
160*     ..
161*     .. Local Scalars ..
162      COMPLEX            TC
163      REAL               ANORM, EPS, TR
164      INTEGER            I, J, K
165*     ..
166*     .. External Functions ..
167      COMPLEX            CDOTC
168      REAL               CLANHE, SLAMCH
169      LOGICAL            LSAME
170      EXTERNAL           CDOTC, CLANHE, SLAMCH, LSAME
171*     ..
172*     .. External Subroutines ..
173      EXTERNAL           CHER, CSCAL, CTRMV
174*     ..
175*     .. Intrinsic Functions ..
176      INTRINSIC          AIMAG, CONJG, REAL
177*     ..
178*     .. Executable Statements ..
179*
180*     Quick exit if N = 0.
181*
182      IF( N.LE.0 ) THEN
183         RESID = ZERO
184         RETURN
185      END IF
186*
187*     Exit with RESID = 1/EPS if ANORM = 0.
188*
189      EPS = SLAMCH( 'Epsilon' )
190      ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
191      IF( ANORM.LE.ZERO ) THEN
192         RESID = ONE / EPS
193         RETURN
194      END IF
195*
196*     Check the imaginary parts of the diagonal elements and return with
197*     an error code if any are nonzero.
198*
199      DO 100 J = 1, N
200         IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
201            RESID = ONE / EPS
202            RETURN
203         END IF
204  100 CONTINUE
205*
206*     Compute the product U'*U, overwriting U.
207*
208      IF( LSAME( UPLO, 'U' ) ) THEN
209*
210         IF( RANK.LT.N ) THEN
211            DO 120 J = RANK + 1, N
212               DO 110 I = RANK + 1, J
213                  AFAC( I, J ) = CZERO
214  110          CONTINUE
215  120       CONTINUE
216         END IF
217*
218         DO 130 K = N, 1, -1
219*
220*           Compute the (K,K) element of the result.
221*
222            TR = CDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
223            AFAC( K, K ) = TR
224*
225*           Compute the rest of column K.
226*
227            CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
228     $                  LDAFAC, AFAC( 1, K ), 1 )
229*
230  130    CONTINUE
231*
232*     Compute the product L*L', overwriting L.
233*
234      ELSE
235*
236         IF( RANK.LT.N ) THEN
237            DO 150 J = RANK + 1, N
238               DO 140 I = J, N
239                  AFAC( I, J ) = CZERO
240  140          CONTINUE
241  150       CONTINUE
242         END IF
243*
244         DO 160 K = N, 1, -1
245*           Add a multiple of column K of the factor L to each of
246*           columns K+1 through N.
247*
248            IF( K+1.LE.N )
249     $         CALL CHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
250     $                    AFAC( K+1, K+1 ), LDAFAC )
251*
252*           Scale column K by the diagonal element.
253*
254            TC = AFAC( K, K )
255            CALL CSCAL( N-K+1, TC, AFAC( K, K ), 1 )
256  160    CONTINUE
257*
258      END IF
259*
260*        Form P*L*L'*P' or P*U'*U*P'
261*
262      IF( LSAME( UPLO, 'U' ) ) THEN
263*
264         DO 180 J = 1, N
265            DO 170 I = 1, N
266               IF( PIV( I ).LE.PIV( J ) ) THEN
267                  IF( I.LE.J ) THEN
268                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
269                  ELSE
270                     PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
271                  END IF
272               END IF
273  170       CONTINUE
274  180    CONTINUE
275*
276*
277      ELSE
278*
279         DO 200 J = 1, N
280            DO 190 I = 1, N
281               IF( PIV( I ).GE.PIV( J ) ) THEN
282                  IF( I.GE.J ) THEN
283                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
284                  ELSE
285                     PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
286                  END IF
287               END IF
288  190       CONTINUE
289  200    CONTINUE
290*
291      END IF
292*
293*     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
294*
295      IF( LSAME( UPLO, 'U' ) ) THEN
296         DO 220 J = 1, N
297            DO 210 I = 1, J - 1
298               PERM( I, J ) = PERM( I, J ) - A( I, J )
299  210       CONTINUE
300            PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
301  220    CONTINUE
302      ELSE
303         DO 240 J = 1, N
304            PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
305            DO 230 I = J + 1, N
306               PERM( I, J ) = PERM( I, J ) - A( I, J )
307  230       CONTINUE
308  240    CONTINUE
309      END IF
310*
311*     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
312*     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
313*
314      RESID = CLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK )
315*
316      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
317*
318      RETURN
319*
320*     End of CPST01
321*
322      END
323