1       SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
  2      $                   LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
  3      $                   RWORK, INFO )
  4 *
  5 *  -- LAPACK driver routine (version 3.3.1) --
  6 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  7 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  8 *  -- April 2011                                                      --
  9 *
 10 *     .. Scalar Arguments ..
 11       CHARACTER          FACT, UPLO
 12       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
 13       DOUBLE PRECISION   RCOND
 14 *     ..
 15 *     .. Array Arguments ..
 16       INTEGER            IPIV( * )
 17       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
 18       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
 19      $                   WORK( * ), X( LDX, * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  ZSYSVX uses the diagonal pivoting factorization to compute the
 26 *  solution to a complex system of linear equations A * X = B,
 27 *  where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
 28 *  matrices.
 29 *
 30 *  Error bounds on the solution and a condition estimate are also
 31 *  provided.
 32 *
 33 *  Description
 34 *  ===========
 35 *
 36 *  The following steps are performed:
 37 *
 38 *  1. If FACT = 'N', the diagonal pivoting method is used to factor A.
 39 *     The form of the factorization is
 40 *        A = U * D * U**T,  if UPLO = 'U', or
 41 *        A = L * D * L**T,  if UPLO = 'L',
 42 *     where U (or L) is a product of permutation and unit upper (lower)
 43 *     triangular matrices, and D is symmetric and block diagonal with
 44 *     1-by-1 and 2-by-2 diagonal blocks.
 45 *
 46 *  2. If some D(i,i)=0, so that D is exactly singular, then the routine
 47 *     returns with INFO = i. Otherwise, the factored form of A is used
 48 *     to estimate the condition number of the matrix A.  If the
 49 *     reciprocal of the condition number is less than machine precision,
 50 *     INFO = N+1 is returned as a warning, but the routine still goes on
 51 *     to solve for X and compute error bounds as described below.
 52 *
 53 *  3. The system of equations is solved for X using the factored form
 54 *     of A.
 55 *
 56 *  4. Iterative refinement is applied to improve the computed solution
 57 *     matrix and calculate error bounds and backward error estimates
 58 *     for it.
 59 *
 60 *  Arguments
 61 *  =========
 62 *
 63 *  FACT    (input) CHARACTER*1
 64 *          Specifies whether or not the factored form of A has been
 65 *          supplied on entry.
 66 *          = 'F':  On entry, AF and IPIV contain the factored form
 67 *                  of A.  A, AF and IPIV will not be modified.
 68 *          = 'N':  The matrix A will be copied to AF and factored.
 69 *
 70 *  UPLO    (input) CHARACTER*1
 71 *          = 'U':  Upper triangle of A is stored;
 72 *          = 'L':  Lower triangle of A is stored.
 73 *
 74 *  N       (input) INTEGER
 75 *          The number of linear equations, i.e., the order of the
 76 *          matrix A.  N >= 0.
 77 *
 78 *  NRHS    (input) INTEGER
 79 *          The number of right hand sides, i.e., the number of columns
 80 *          of the matrices B and X.  NRHS >= 0.
 81 *
 82 *  A       (input) COMPLEX*16 array, dimension (LDA,N)
 83 *          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
 84 *          upper triangular part of A contains the upper triangular part
 85 *          of the matrix A, and the strictly lower triangular part of A
 86 *          is not referenced.  If UPLO = 'L', the leading N-by-N lower
 87 *          triangular part of A contains the lower triangular part of
 88 *          the matrix A, and the strictly upper triangular part of A is
 89 *          not referenced.
 90 *
 91 *  LDA     (input) INTEGER
 92 *          The leading dimension of the array A.  LDA >= max(1,N).
 93 *
 94 *  AF      (input or output) COMPLEX*16 array, dimension (LDAF,N)
 95 *          If FACT = 'F', then AF is an input argument and on entry
 96 *          contains the block diagonal matrix D and the multipliers used
 97 *          to obtain the factor U or L from the factorization
 98 *          A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF.
 99 *
100 *          If FACT = 'N', then AF is an output argument and on exit
101 *          returns the block diagonal matrix D and the multipliers used
102 *          to obtain the factor U or L from the factorization
103 *          A = U*D*U**T or A = L*D*L**T.
104 *
105 *  LDAF    (input) INTEGER
106 *          The leading dimension of the array AF.  LDAF >= max(1,N).
107 *
108 *  IPIV    (input or output) INTEGER array, dimension (N)
109 *          If FACT = 'F', then IPIV is an input argument and on entry
110 *          contains details of the interchanges and the block structure
111 *          of D, as determined by ZSYTRF.
112 *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
113 *          interchanged and D(k,k) is a 1-by-1 diagonal block.
114 *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
115 *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
116 *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
117 *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
118 *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
119 *
120 *          If FACT = 'N', then IPIV is an output argument and on exit
121 *          contains details of the interchanges and the block structure
122 *          of D, as determined by ZSYTRF.
123 *
124 *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
125 *          The N-by-NRHS right hand side matrix B.
126 *
127 *  LDB     (input) INTEGER
128 *          The leading dimension of the array B.  LDB >= max(1,N).
129 *
130 *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
131 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
132 *
133 *  LDX     (input) INTEGER
134 *          The leading dimension of the array X.  LDX >= max(1,N).
135 *
136 *  RCOND   (output) DOUBLE PRECISION
137 *          The estimate of the reciprocal condition number of the matrix
138 *          A.  If RCOND is less than the machine precision (in
139 *          particular, if RCOND = 0), the matrix is singular to working
140 *          precision.  This condition is indicated by a return code of
141 *          INFO > 0.
142 *
143 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
144 *          The estimated forward error bound for each solution vector
145 *          X(j) (the j-th column of the solution matrix X).
146 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
147 *          is an estimated upper bound for the magnitude of the largest
148 *          element in (X(j) - XTRUE) divided by the magnitude of the
149 *          largest element in X(j).  The estimate is as reliable as
150 *          the estimate for RCOND, and is almost always a slight
151 *          overestimate of the true error.
152 *
153 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
154 *          The componentwise relative backward error of each solution
155 *          vector X(j) (i.e., the smallest relative change in
156 *          any element of A or B that makes X(j) an exact solution).
157 *
158 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
159 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
160 *
161 *  LWORK   (input) INTEGER
162 *          The length of WORK.  LWORK >= max(1,2*N), and for best
163 *          performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
164 *          NB is the optimal blocksize for ZSYTRF.
165 *
166 *          If LWORK = -1, then a workspace query is assumed; the routine
167 *          only calculates the optimal size of the WORK array, returns
168 *          this value as the first entry of the WORK array, and no error
169 *          message related to LWORK is issued by XERBLA.
170 *
171 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
172 *
173 *  INFO    (output) INTEGER
174 *          = 0: successful exit
175 *          < 0: if INFO = -i, the i-th argument had an illegal value
176 *          > 0: if INFO = i, and i is
177 *                <= N:  D(i,i) is exactly zero.  The factorization
178 *                       has been completed but the factor D is exactly
179 *                       singular, so the solution and error bounds could
180 *                       not be computed. RCOND = 0 is returned.
181 *                = N+1: D is nonsingular, but RCOND is less than machine
182 *                       precision, meaning that the matrix is singular
183 *                       to working precision.  Nevertheless, the
184 *                       solution and error bounds are computed because
185 *                       there are a number of situations where the
186 *                       computed solution can be more accurate than the
187 *                       value of RCOND would suggest.
188 *
189 *  =====================================================================
190 *
191 *     .. Parameters ..
192       DOUBLE PRECISION   ZERO
193       PARAMETER          ( ZERO = 0.0D+0 )
194 *     ..
195 *     .. Local Scalars ..
196       LOGICAL            LQUERY, NOFACT
197       INTEGER            LWKOPT, NB
198       DOUBLE PRECISION   ANORM
199 *     ..
200 *     .. External Functions ..
201       LOGICAL            LSAME
202       INTEGER            ILAENV
203       DOUBLE PRECISION   DLAMCH, ZLANSY
204       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANSY
205 *     ..
206 *     .. External Subroutines ..
207       EXTERNAL           XERBLA, ZLACPY, ZSYCON, ZSYRFS, ZSYTRF, ZSYTRS
208 *     ..
209 *     .. Intrinsic Functions ..
210       INTRINSIC          MAX
211 *     ..
212 *     .. Executable Statements ..
213 *
214 *     Test the input parameters.
215 *
216       INFO = 0
217       NOFACT = LSAME( FACT, 'N' )
218       LQUERY = ( LWORK.EQ.-1 )
219       IF.NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
220          INFO = -1
221       ELSE IF.NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
222      $          THEN
223          INFO = -2
224       ELSE IF( N.LT.0 ) THEN
225          INFO = -3
226       ELSE IF( NRHS.LT.0 ) THEN
227          INFO = -4
228       ELSE IF( LDA.LT.MAX1, N ) ) THEN
229          INFO = -6
230       ELSE IF( LDAF.LT.MAX1, N ) ) THEN
231          INFO = -8
232       ELSE IF( LDB.LT.MAX1, N ) ) THEN
233          INFO = -11
234       ELSE IF( LDX.LT.MAX1, N ) ) THEN
235          INFO = -13
236       ELSE IF( LWORK.LT.MAX12*N ) .AND. .NOT.LQUERY ) THEN
237          INFO = -18
238       END IF
239 *
240       IF( INFO.EQ.0 ) THEN
241          LWKOPT = MAX12*N )
242          IF( NOFACT ) THEN
243             NB = ILAENV( 1'ZSYTRF', UPLO, N, -1-1-1 )
244             LWKOPT = MAX( LWKOPT, N*NB )
245          END IF
246          WORK( 1 ) = LWKOPT
247       END IF
248 *
249       IF( INFO.NE.0 ) THEN
250          CALL XERBLA( 'ZSYSVX'-INFO )
251          RETURN
252       ELSE IF( LQUERY ) THEN
253          RETURN
254       END IF
255 *
256       IF( NOFACT ) THEN
257 *
258 *        Compute the factorization A = U*D*U**T or A = L*D*L**T.
259 *
260          CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
261          CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
262 *
263 *        Return if INFO is non-zero.
264 *
265          IF( INFO.GT.0 )THEN
266             RCOND = ZERO
267             RETURN
268          END IF
269       END IF
270 *
271 *     Compute the norm of the matrix A.
272 *
273       ANORM = ZLANSY( 'I', UPLO, N, A, LDA, RWORK )
274 *
275 *     Compute the reciprocal of the condition number of A.
276 *
277       CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
278 *
279 *     Compute the solution vectors X.
280 *
281       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
282       CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
283 *
284 *     Use iterative refinement to improve the computed solutions and
285 *     compute error bounds and backward error estimates for them.
286 *
287       CALL ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
288      $             LDX, FERR, BERR, WORK, RWORK, INFO )
289 *
290 *     Set INFO = N+1 if the matrix is singular to working precision.
291 *
292       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
293      $   INFO = N + 1
294 *
295       WORK( 1 ) = LWKOPT
296 *
297       RETURN
298 *
299 *     End of ZSYSVX
300 *
301       END