LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zcgesv.f
Go to the documentation of this file.
1 *> \brief <b> ZCGESV computes the solution to system of linear equations A * X = B for GE matrices</b> (mixed precision with iterative refinement)
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZCGESV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zcgesv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zcgesv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zcgesv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
22 * SWORK, RWORK, ITER, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * DOUBLE PRECISION RWORK( * )
30 * COMPLEX SWORK( * )
31 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
32 * $ X( LDX, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZCGESV computes the solution to a complex system of linear equations
42 *> A * X = B,
43 *> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
44 *>
45 *> ZCGESV first attempts to factorize the matrix in COMPLEX and use this
46 *> factorization within an iterative refinement procedure to produce a
47 *> solution with COMPLEX*16 normwise backward error quality (see below).
48 *> If the approach fails the method switches to a COMPLEX*16
49 *> factorization and solve.
50 *>
51 *> The iterative refinement is not going to be a winning strategy if
52 *> the ratio COMPLEX performance over COMPLEX*16 performance is too
53 *> small. A reasonable strategy should take the number of right-hand
54 *> sides and the size of the matrix into account. This might be done
55 *> with a call to ILAENV in the future. Up to now, we always try
56 *> iterative refinement.
57 *>
58 *> The iterative refinement process is stopped if
59 *> ITER > ITERMAX
60 *> or for all the RHS we have:
61 *> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
62 *> where
63 *> o ITER is the number of the current iteration in the iterative
64 *> refinement process
65 *> o RNRM is the infinity-norm of the residual
66 *> o XNRM is the infinity-norm of the solution
67 *> o ANRM is the infinity-operator-norm of the matrix A
68 *> o EPS is the machine epsilon returned by DLAMCH('Epsilon')
69 *> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
70 *> respectively.
71 *> \endverbatim
72 *
73 * Arguments:
74 * ==========
75 *
76 *> \param[in] N
77 *> \verbatim
78 *> N is INTEGER
79 *> The number of linear equations, i.e., the order of the
80 *> matrix A. N >= 0.
81 *> \endverbatim
82 *>
83 *> \param[in] NRHS
84 *> \verbatim
85 *> NRHS is INTEGER
86 *> The number of right hand sides, i.e., the number of columns
87 *> of the matrix B. NRHS >= 0.
88 *> \endverbatim
89 *>
90 *> \param[in,out] A
91 *> \verbatim
92 *> A is COMPLEX*16 array,
93 *> dimension (LDA,N)
94 *> On entry, the N-by-N coefficient matrix A.
95 *> On exit, if iterative refinement has been successfully used
96 *> (INFO = 0 and ITER >= 0, see description below), then A is
97 *> unchanged, if double precision factorization has been used
98 *> (INFO = 0 and ITER < 0, see description below), then the
99 *> array A contains the factors L and U from the factorization
100 *> A = P*L*U; the unit diagonal elements of L are not stored.
101 *> \endverbatim
102 *>
103 *> \param[in] LDA
104 *> \verbatim
105 *> LDA is INTEGER
106 *> The leading dimension of the array A. LDA >= max(1,N).
107 *> \endverbatim
108 *>
109 *> \param[out] IPIV
110 *> \verbatim
111 *> IPIV is INTEGER array, dimension (N)
112 *> The pivot indices that define the permutation matrix P;
113 *> row i of the matrix was interchanged with row IPIV(i).
114 *> Corresponds either to the single precision factorization
115 *> (if INFO = 0 and ITER >= 0) or the double precision
116 *> factorization (if INFO = 0 and ITER < 0).
117 *> \endverbatim
118 *>
119 *> \param[in] B
120 *> \verbatim
121 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
122 *> The N-by-NRHS right hand side matrix B.
123 *> \endverbatim
124 *>
125 *> \param[in] LDB
126 *> \verbatim
127 *> LDB is INTEGER
128 *> The leading dimension of the array B. LDB >= max(1,N).
129 *> \endverbatim
130 *>
131 *> \param[out] X
132 *> \verbatim
133 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
134 *> If INFO = 0, the N-by-NRHS solution matrix X.
135 *> \endverbatim
136 *>
137 *> \param[in] LDX
138 *> \verbatim
139 *> LDX is INTEGER
140 *> The leading dimension of the array X. LDX >= max(1,N).
141 *> \endverbatim
142 *>
143 *> \param[out] WORK
144 *> \verbatim
145 *> WORK is COMPLEX*16 array, dimension (N,NRHS)
146 *> This array is used to hold the residual vectors.
147 *> \endverbatim
148 *>
149 *> \param[out] SWORK
150 *> \verbatim
151 *> SWORK is COMPLEX array, dimension (N*(N+NRHS))
152 *> This array is used to use the single precision matrix and the
153 *> right-hand sides or solutions in single precision.
154 *> \endverbatim
155 *>
156 *> \param[out] RWORK
157 *> \verbatim
158 *> RWORK is DOUBLE PRECISION array, dimension (N)
159 *> \endverbatim
160 *>
161 *> \param[out] ITER
162 *> \verbatim
163 *> ITER is INTEGER
164 *> < 0: iterative refinement has failed, COMPLEX*16
165 *> factorization has been performed
166 *> -1 : the routine fell back to full precision for
167 *> implementation- or machine-specific reasons
168 *> -2 : narrowing the precision induced an overflow,
169 *> the routine fell back to full precision
170 *> -3 : failure of CGETRF
171 *> -31: stop the iterative refinement after the 30th
172 *> iterations
173 *> > 0: iterative refinement has been successfully used.
174 *> Returns the number of iterations
175 *> \endverbatim
176 *>
177 *> \param[out] INFO
178 *> \verbatim
179 *> INFO is INTEGER
180 *> = 0: successful exit
181 *> < 0: if INFO = -i, the i-th argument had an illegal value
182 *> > 0: if INFO = i, U(i,i) computed in COMPLEX*16 is exactly
183 *> zero. The factorization has been completed, but the
184 *> factor U is exactly singular, so the solution
185 *> could not be computed.
186 *> \endverbatim
187 *
188 * Authors:
189 * ========
190 *
191 *> \author Univ. of Tennessee
192 *> \author Univ. of California Berkeley
193 *> \author Univ. of Colorado Denver
194 *> \author NAG Ltd.
195 *
196 *> \ingroup complex16GEsolve
197 *
198 * =====================================================================
199  SUBROUTINE zcgesv( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
200  $ SWORK, RWORK, ITER, INFO )
201 *
202 * -- LAPACK driver routine --
203 * -- LAPACK is a software package provided by Univ. of Tennessee, --
204 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205 *
206 * .. Scalar Arguments ..
207  INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
208 * ..
209 * .. Array Arguments ..
210  INTEGER IPIV( * )
211  DOUBLE PRECISION RWORK( * )
212  COMPLEX SWORK( * )
213  COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
214  $ x( ldx, * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  LOGICAL DOITREF
221  parameter( doitref = .true. )
222 *
223  INTEGER ITERMAX
224  parameter( itermax = 30 )
225 *
226  DOUBLE PRECISION BWDMAX
227  parameter( bwdmax = 1.0e+00 )
228 *
229  COMPLEX*16 NEGONE, ONE
230  parameter( negone = ( -1.0d+00, 0.0d+00 ),
231  $ one = ( 1.0d+00, 0.0d+00 ) )
232 *
233 * .. Local Scalars ..
234  INTEGER I, IITER, PTSA, PTSX
235  DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
236  COMPLEX*16 ZDUM
237 *
238 * .. External Subroutines ..
239  EXTERNAL cgetrs, cgetrf, clag2z, xerbla, zaxpy, zgemm,
241 * ..
242 * .. External Functions ..
243  INTEGER IZAMAX
244  DOUBLE PRECISION DLAMCH, ZLANGE
245  EXTERNAL izamax, dlamch, zlange
246 * ..
247 * .. Intrinsic Functions ..
248  INTRINSIC abs, dble, max, sqrt
249 * ..
250 * .. Statement Functions ..
251  DOUBLE PRECISION CABS1
252 * ..
253 * .. Statement Function definitions ..
254  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
255 * ..
256 * .. Executable Statements ..
257 *
258  info = 0
259  iter = 0
260 *
261 * Test the input parameters.
262 *
263  IF( n.LT.0 ) THEN
264  info = -1
265  ELSE IF( nrhs.LT.0 ) THEN
266  info = -2
267  ELSE IF( lda.LT.max( 1, n ) ) THEN
268  info = -4
269  ELSE IF( ldb.LT.max( 1, n ) ) THEN
270  info = -7
271  ELSE IF( ldx.LT.max( 1, n ) ) THEN
272  info = -9
273  END IF
274  IF( info.NE.0 ) THEN
275  CALL xerbla( 'ZCGESV', -info )
276  RETURN
277  END IF
278 *
279 * Quick return if (N.EQ.0).
280 *
281  IF( n.EQ.0 )
282  $ RETURN
283 *
284 * Skip single precision iterative refinement if a priori slower
285 * than double precision factorization.
286 *
287  IF( .NOT.doitref ) THEN
288  iter = -1
289  GO TO 40
290  END IF
291 *
292 * Compute some constants.
293 *
294  anrm = zlange( 'I', n, n, a, lda, rwork )
295  eps = dlamch( 'Epsilon' )
296  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
297 *
298 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
299 *
300  ptsa = 1
301  ptsx = ptsa + n*n
302 *
303 * Convert B from double precision to single precision and store the
304 * result in SX.
305 *
306  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
307 *
308  IF( info.NE.0 ) THEN
309  iter = -2
310  GO TO 40
311  END IF
312 *
313 * Convert A from double precision to single precision and store the
314 * result in SA.
315 *
316  CALL zlag2c( n, n, a, lda, swork( ptsa ), n, info )
317 *
318  IF( info.NE.0 ) THEN
319  iter = -2
320  GO TO 40
321  END IF
322 *
323 * Compute the LU factorization of SA.
324 *
325  CALL cgetrf( n, n, swork( ptsa ), n, ipiv, info )
326 *
327  IF( info.NE.0 ) THEN
328  iter = -3
329  GO TO 40
330  END IF
331 *
332 * Solve the system SA*SX = SB.
333 *
334  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
335  $ swork( ptsx ), n, info )
336 *
337 * Convert SX back to double precision
338 *
339  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
340 *
341 * Compute R = B - AX (R is WORK).
342 *
343  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
344 *
345  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
346  $ lda, x, ldx, one, work, n )
347 *
348 * Check whether the NRHS normwise backward errors satisfy the
349 * stopping criterion. If yes, set ITER=0 and return.
350 *
351  DO i = 1, nrhs
352  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
353  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
354  IF( rnrm.GT.xnrm*cte )
355  $ GO TO 10
356  END DO
357 *
358 * If we are here, the NRHS normwise backward errors satisfy the
359 * stopping criterion. We are good to exit.
360 *
361  iter = 0
362  RETURN
363 *
364  10 CONTINUE
365 *
366  DO 30 iiter = 1, itermax
367 *
368 * Convert R (in WORK) from double precision to single precision
369 * and store the result in SX.
370 *
371  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
372 *
373  IF( info.NE.0 ) THEN
374  iter = -2
375  GO TO 40
376  END IF
377 *
378 * Solve the system SA*SX = SR.
379 *
380  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
381  $ swork( ptsx ), n, info )
382 *
383 * Convert SX back to double precision and update the current
384 * iterate.
385 *
386  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
387 *
388  DO i = 1, nrhs
389  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
390  END DO
391 *
392 * Compute R = B - AX (R is WORK).
393 *
394  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
395 *
396  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
397  $ a, lda, x, ldx, one, work, n )
398 *
399 * Check whether the NRHS normwise backward errors satisfy the
400 * stopping criterion. If yes, set ITER=IITER>0 and return.
401 *
402  DO i = 1, nrhs
403  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
404  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
405  IF( rnrm.GT.xnrm*cte )
406  $ GO TO 20
407  END DO
408 *
409 * If we are here, the NRHS normwise backward errors satisfy the
410 * stopping criterion, we are good to exit.
411 *
412  iter = iiter
413 *
414  RETURN
415 *
416  20 CONTINUE
417 *
418  30 CONTINUE
419 *
420 * If we are at this place of the code, this is because we have
421 * performed ITER=ITERMAX iterations and never satisfied the stopping
422 * criterion, set up the ITER flag accordingly and follow up on double
423 * precision routine.
424 *
425  iter = -itermax - 1
426 *
427  40 CONTINUE
428 *
429 * Single-precision iterative refinement failed to converge to a
430 * satisfactory solution, so we resort to double precision.
431 *
432  CALL zgetrf( n, n, a, lda, ipiv, info )
433 *
434  IF( info.NE.0 )
435  $ RETURN
436 *
437  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
438  CALL zgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
439  $ info )
440 *
441  RETURN
442 *
443 * End of ZCGESV
444 *
445  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:121
subroutine zcgesv(N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK, SWORK, RWORK, ITER, INFO)
ZCGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision...
Definition: zcgesv.f:201
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlag2c(M, N, A, LDA, SA, LDSA, INFO)
ZLAG2C converts a complex double precision matrix to a complex single precision matrix.
Definition: zlag2c.f:107
subroutine clag2z(M, N, SA, LDSA, A, LDA, INFO)
CLAG2Z converts a complex single precision matrix to a complex double precision matrix.
Definition: clag2z.f:103
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:121
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:108
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102