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*> \date December 2016
132*
133*> \ingroup complex_lin
134*
135*  =====================================================================
136      SUBROUTINE CPST01( UPLO, N, A, LDA, AFAC, LDAFAC, PERM, LDPERM,
137     $                   PIV, RWORK, RESID, RANK )
138*
139*  -- LAPACK test routine (version 3.7.0) --
140*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
141*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142*     December 2016
143*
144*     .. Scalar Arguments ..
145      REAL               RESID
146      INTEGER            LDA, LDAFAC, LDPERM, N, RANK
147      CHARACTER          UPLO
148*     ..
149*     .. Array Arguments ..
150      COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ),
151     $                   PERM( LDPERM, * )
152      REAL               RWORK( * )
153      INTEGER            PIV( * )
154*     ..
155*
156*  =====================================================================
157*
158*     .. Parameters ..
159      REAL               ZERO, ONE
160      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
161      COMPLEX            CZERO
162      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
163*     ..
164*     .. Local Scalars ..
165      COMPLEX            TC
166      REAL               ANORM, EPS, TR
167      INTEGER            I, J, K
168*     ..
169*     .. External Functions ..
170      COMPLEX            CDOTC
171      REAL               CLANHE, SLAMCH
172      LOGICAL            LSAME
173      EXTERNAL           CDOTC, CLANHE, SLAMCH, LSAME
174*     ..
175*     .. External Subroutines ..
176      EXTERNAL           CHER, CSCAL, CTRMV
177*     ..
178*     .. Intrinsic Functions ..
179      INTRINSIC          AIMAG, CONJG, REAL
180*     ..
181*     .. Executable Statements ..
182*
183*     Quick exit if N = 0.
184*
185      IF( N.LE.0 ) THEN
186         RESID = ZERO
187         RETURN
188      END IF
189*
190*     Exit with RESID = 1/EPS if ANORM = 0.
191*
192      EPS = SLAMCH( 'Epsilon' )
193      ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
194      IF( ANORM.LE.ZERO ) THEN
195         RESID = ONE / EPS
196         RETURN
197      END IF
198*
199*     Check the imaginary parts of the diagonal elements and return with
200*     an error code if any are nonzero.
201*
202      DO 100 J = 1, N
203         IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
204            RESID = ONE / EPS
205            RETURN
206         END IF
207  100 CONTINUE
208*
209*     Compute the product U'*U, overwriting U.
210*
211      IF( LSAME( UPLO, 'U' ) ) THEN
212*
213         IF( RANK.LT.N ) THEN
214            DO 120 J = RANK + 1, N
215               DO 110 I = RANK + 1, J
216                  AFAC( I, J ) = CZERO
217  110          CONTINUE
218  120       CONTINUE
219         END IF
220*
221         DO 130 K = N, 1, -1
222*
223*           Compute the (K,K) element of the result.
224*
225            TR = CDOTC( K, AFAC( 1, K ), 1, AFAC( 1, K ), 1 )
226            AFAC( K, K ) = TR
227*
228*           Compute the rest of column K.
229*
230            CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', K-1, AFAC,
231     $                  LDAFAC, AFAC( 1, K ), 1 )
232*
233  130    CONTINUE
234*
235*     Compute the product L*L', overwriting L.
236*
237      ELSE
238*
239         IF( RANK.LT.N ) THEN
240            DO 150 J = RANK + 1, N
241               DO 140 I = J, N
242                  AFAC( I, J ) = CZERO
243  140          CONTINUE
244  150       CONTINUE
245         END IF
246*
247         DO 160 K = N, 1, -1
248*           Add a multiple of column K of the factor L to each of
249*           columns K+1 through N.
250*
251            IF( K+1.LE.N )
252     $         CALL CHER( 'Lower', N-K, ONE, AFAC( K+1, K ), 1,
253     $                    AFAC( K+1, K+1 ), LDAFAC )
254*
255*           Scale column K by the diagonal element.
256*
257            TC = AFAC( K, K )
258            CALL CSCAL( N-K+1, TC, AFAC( K, K ), 1 )
259  160    CONTINUE
260*
261      END IF
262*
263*        Form P*L*L'*P' or P*U'*U*P'
264*
265      IF( LSAME( UPLO, 'U' ) ) THEN
266*
267         DO 180 J = 1, N
268            DO 170 I = 1, N
269               IF( PIV( I ).LE.PIV( J ) ) THEN
270                  IF( I.LE.J ) THEN
271                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
272                  ELSE
273                     PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
274                  END IF
275               END IF
276  170       CONTINUE
277  180    CONTINUE
278*
279*
280      ELSE
281*
282         DO 200 J = 1, N
283            DO 190 I = 1, N
284               IF( PIV( I ).GE.PIV( J ) ) THEN
285                  IF( I.GE.J ) THEN
286                     PERM( PIV( I ), PIV( J ) ) = AFAC( I, J )
287                  ELSE
288                     PERM( PIV( I ), PIV( J ) ) = CONJG( AFAC( J, I ) )
289                  END IF
290               END IF
291  190       CONTINUE
292  200    CONTINUE
293*
294      END IF
295*
296*     Compute the difference  P*L*L'*P' - A (or P*U'*U*P' - A).
297*
298      IF( LSAME( UPLO, 'U' ) ) THEN
299         DO 220 J = 1, N
300            DO 210 I = 1, J - 1
301               PERM( I, J ) = PERM( I, J ) - A( I, J )
302  210       CONTINUE
303            PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
304  220    CONTINUE
305      ELSE
306         DO 240 J = 1, N
307            PERM( J, J ) = PERM( J, J ) - REAL( A( J, J ) )
308            DO 230 I = J + 1, N
309               PERM( I, J ) = PERM( I, J ) - A( I, J )
310  230       CONTINUE
311  240    CONTINUE
312      END IF
313*
314*     Compute norm( P*L*L'P - A ) / ( N * norm(A) * EPS ), or
315*     ( P*U'*U*P' - A )/ ( N * norm(A) * EPS ).
316*
317      RESID = CLANHE( '1', UPLO, N, PERM, LDAFAC, RWORK )
318*
319      RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
320*
321      RETURN
322*
323*     End of CPST01
324*
325      END
326