1*> \brief <b> ZSYSV_RK computes the solution to system of linear equations A * X = B for SY matrices</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZSYSV_RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsysv_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsysv_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsysv_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZSYSV_RK( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
22*                            WORK, LWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          UPLO
26*       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IPIV( * )
30*       COMPLEX*16         A( LDA, * ), B( LDB, * ), E( * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*> ZSYSV_RK computes the solution to a complex system of linear
39*> equations A * X = B, where A is an N-by-N symmetric matrix
40*> and X and B are N-by-NRHS matrices.
41*>
42*> The bounded Bunch-Kaufman (rook) diagonal pivoting method is used
43*> to factor A as
44*>    A = P*U*D*(U**T)*(P**T),  if UPLO = 'U', or
45*>    A = P*L*D*(L**T)*(P**T),  if UPLO = 'L',
46*> where U (or L) is unit upper (or lower) triangular matrix,
47*> U**T (or L**T) is the transpose of U (or L), P is a permutation
48*> matrix, P**T is the transpose of P, and D is symmetric and block
49*> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
50*>
51*> ZSYTRF_RK is called to compute the factorization of a complex
52*> symmetric matrix.  The factored form of A is then used to solve
53*> the system of equations A * X = B by calling BLAS3 routine ZSYTRS_3.
54*> \endverbatim
55*
56*  Arguments:
57*  ==========
58*
59*> \param[in] UPLO
60*> \verbatim
61*>          UPLO is CHARACTER*1
62*>          Specifies whether the upper or lower triangular part of the
63*>          symmetric matrix A is stored:
64*>          = 'U':  Upper triangle of A is stored;
65*>          = 'L':  Lower triangle of A is stored.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*>          N is INTEGER
71*>          The number of linear equations, i.e., the order of the
72*>          matrix A.  N >= 0.
73*> \endverbatim
74*>
75*> \param[in] NRHS
76*> \verbatim
77*>          NRHS is INTEGER
78*>          The number of right hand sides, i.e., the number of columns
79*>          of the matrix B.  NRHS >= 0.
80*> \endverbatim
81*>
82*> \param[in,out] A
83*> \verbatim
84*>          A is COMPLEX*16 array, dimension (LDA,N)
85*>          On entry, the symmetric matrix A.
86*>            If UPLO = 'U': the leading N-by-N upper triangular part
87*>            of A contains the upper triangular part of the matrix A,
88*>            and the strictly lower triangular part of A is not
89*>            referenced.
90*>
91*>            If UPLO = 'L': the leading N-by-N lower triangular part
92*>            of A contains the lower triangular part of the matrix A,
93*>            and the strictly upper triangular part of A is not
94*>            referenced.
95*>
96*>          On exit, if INFO = 0, diagonal of the block diagonal
97*>          matrix D and factors U or L  as computed by ZSYTRF_RK:
98*>            a) ONLY diagonal elements of the symmetric block diagonal
99*>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
100*>               (superdiagonal (or subdiagonal) elements of D
101*>                are stored on exit in array E), and
102*>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
103*>               If UPLO = 'L': factor L in the subdiagonal part of A.
104*>
105*>          For more info see the description of ZSYTRF_RK routine.
106*> \endverbatim
107*>
108*> \param[in] LDA
109*> \verbatim
110*>          LDA is INTEGER
111*>          The leading dimension of the array A.  LDA >= max(1,N).
112*> \endverbatim
113*>
114*> \param[out] E
115*> \verbatim
116*>          E is COMPLEX*16 array, dimension (N)
117*>          On exit, contains the output computed by the factorization
118*>          routine ZSYTRF_RK, i.e. the superdiagonal (or subdiagonal)
119*>          elements of the symmetric block diagonal matrix D
120*>          with 1-by-1 or 2-by-2 diagonal blocks, where
121*>          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
122*>          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
123*>
124*>          NOTE: For 1-by-1 diagonal block D(k), where
125*>          1 <= k <= N, the element E(k) is set to 0 in both
126*>          UPLO = 'U' or UPLO = 'L' cases.
127*>
128*>          For more info see the description of ZSYTRF_RK routine.
129*> \endverbatim
130*>
131*> \param[out] IPIV
132*> \verbatim
133*>          IPIV is INTEGER array, dimension (N)
134*>          Details of the interchanges and the block structure of D,
135*>          as determined by ZSYTRF_RK.
136*>
137*>          For more info see the description of ZSYTRF_RK routine.
138*> \endverbatim
139*>
140*> \param[in,out] B
141*> \verbatim
142*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
143*>          On entry, the N-by-NRHS right hand side matrix B.
144*>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
145*> \endverbatim
146*>
147*> \param[in] LDB
148*> \verbatim
149*>          LDB is INTEGER
150*>          The leading dimension of the array B.  LDB >= max(1,N).
151*> \endverbatim
152*>
153*> \param[out] WORK
154*> \verbatim
155*>          WORK is COMPLEX*16 array, dimension ( MAX(1,LWORK) ).
156*>          Work array used in the factorization stage.
157*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
158*> \endverbatim
159*>
160*> \param[in] LWORK
161*> \verbatim
162*>          LWORK is INTEGER
163*>          The length of WORK.  LWORK >= 1. For best performance
164*>          of factorization stage LWORK >= max(1,N*NB), where NB is
165*>          the optimal blocksize for ZSYTRF_RK.
166*>
167*>          If LWORK = -1, then a workspace query is assumed;
168*>          the routine only calculates the optimal size of the WORK
169*>          array for factorization stage, returns this value as
170*>          the first entry of the WORK array, and no error message
171*>          related to LWORK is issued by XERBLA.
172*> \endverbatim
173*>
174*> \param[out] INFO
175*> \verbatim
176*>          INFO is INTEGER
177*>          = 0: successful exit
178*>
179*>          < 0: If INFO = -k, the k-th argument had an illegal value
180*>
181*>          > 0: If INFO = k, the matrix A is singular, because:
182*>                 If UPLO = 'U': column k in the upper
183*>                 triangular part of A contains all zeros.
184*>                 If UPLO = 'L': column k in the lower
185*>                 triangular part of A contains all zeros.
186*>
187*>               Therefore D(k,k) is exactly zero, and superdiagonal
188*>               elements of column k of U (or subdiagonal elements of
189*>               column k of L ) are all zeros. The factorization has
190*>               been completed, but the block diagonal matrix D is
191*>               exactly singular, and division by zero will occur if
192*>               it is used to solve a system of equations.
193*>
194*>               NOTE: INFO only stores the first occurrence of
195*>               a singularity, any subsequent occurrence of singularity
196*>               is not stored in INFO even though the factorization
197*>               always completes.
198*> \endverbatim
199*
200*  Authors:
201*  ========
202*
203*> \author Univ. of Tennessee
204*> \author Univ. of California Berkeley
205*> \author Univ. of Colorado Denver
206*> \author NAG Ltd.
207*
208*> \ingroup complex16SYsolve
209*
210*> \par Contributors:
211*  ==================
212*>
213*> \verbatim
214*>
215*>  December 2016,  Igor Kozachenko,
216*>                  Computer Science Division,
217*>                  University of California, Berkeley
218*>
219*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
220*>                  School of Mathematics,
221*>                  University of Manchester
222*>
223*> \endverbatim
224*
225*  =====================================================================
226      SUBROUTINE ZSYSV_RK( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, WORK,
227     $                     LWORK, INFO )
228*
229*  -- LAPACK driver routine --
230*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
231*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232*
233*     .. Scalar Arguments ..
234      CHARACTER          UPLO
235      INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
236*     ..
237*     .. Array Arguments ..
238      INTEGER            IPIV( * )
239      COMPLEX*16         A( LDA, * ), B( LDB, * ), E( * ), WORK( * )
240*     ..
241*
242*  =====================================================================
243*
244*     .. Local Scalars ..
245      LOGICAL            LQUERY
246      INTEGER            LWKOPT
247*     ..
248*     .. External Functions ..
249      LOGICAL            LSAME
250      EXTERNAL           LSAME
251*     ..
252*     .. External Subroutines ..
253      EXTERNAL           XERBLA, ZSYTRF_RK, ZSYTRS_3
254*     ..
255*     .. Intrinsic Functions ..
256      INTRINSIC          MAX
257*     ..
258*     .. Executable Statements ..
259*
260*     Test the input parameters.
261*
262      INFO = 0
263      LQUERY = ( LWORK.EQ.-1 )
264      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
265         INFO = -1
266      ELSE IF( N.LT.0 ) THEN
267         INFO = -2
268      ELSE IF( NRHS.LT.0 ) THEN
269         INFO = -3
270      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
271         INFO = -5
272      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
273         INFO = -9
274      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
275         INFO = -11
276      END IF
277*
278      IF( INFO.EQ.0 ) THEN
279         IF( N.EQ.0 ) THEN
280            LWKOPT = 1
281         ELSE
282            CALL ZSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, -1, INFO )
283            LWKOPT = DBLE( WORK(1) )
284         END IF
285         WORK( 1 ) = LWKOPT
286      END IF
287*
288      IF( INFO.NE.0 ) THEN
289         CALL XERBLA( 'ZSYSV_RK ', -INFO )
290         RETURN
291      ELSE IF( LQUERY ) THEN
292         RETURN
293      END IF
294*
295*     Compute the factorization A = P*U*D*(U**T)*(P**T) or
296*     A = P*U*D*(U**T)*(P**T).
297*
298      CALL ZSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK, INFO )
299*
300      IF( INFO.EQ.0 ) THEN
301*
302*        Solve the system A*X = B with BLAS3 solver, overwriting B with X.
303*
304         CALL ZSYTRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, INFO )
305*
306      END IF
307*
308      WORK( 1 ) = LWKOPT
309*
310      RETURN
311*
312*     End of ZSYSV_RK
313*
314      END
315