LAPACK  3.6.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
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.EQ.0 and ITER.GE.0, see description below), then A is
97 *> unchanged, if double precision factorization has been used
98 *> (INFO.EQ.0 and ITER.LT.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.EQ.0 and ITER.GE.0) or the double precision
116 *> factorization (if INFO.EQ.0 and ITER.LT.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 *> \date June 2016
197 *
198 *> \ingroup complex16GEsolve
199 *
200 * =====================================================================
201  SUBROUTINE zcgesv( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
202  \$ swork, rwork, iter, info )
203 *
204 * -- LAPACK driver routine (version 3.6.1) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * June 2016
208 *
209 * .. Scalar Arguments ..
210  INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
211 * ..
212 * .. Array Arguments ..
213  INTEGER IPIV( * )
214  DOUBLE PRECISION RWORK( * )
215  COMPLEX SWORK( * )
216  COMPLEX*16 A( lda, * ), B( ldb, * ), WORK( n, * ),
217  \$ x( ldx, * )
218 * ..
219 *
220 * =====================================================================
221 *
222 * .. Parameters ..
223  LOGICAL DOITREF
224  parameter ( doitref = .true. )
225 *
226  INTEGER ITERMAX
227  parameter ( itermax = 30 )
228 *
229  DOUBLE PRECISION BWDMAX
230  parameter ( bwdmax = 1.0e+00 )
231 *
232  COMPLEX*16 NEGONE, ONE
233  parameter ( negone = ( -1.0d+00, 0.0d+00 ),
234  \$ one = ( 1.0d+00, 0.0d+00 ) )
235 *
236 * .. Local Scalars ..
237  INTEGER I, IITER, PTSA, PTSX
238  DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
239  COMPLEX*16 ZDUM
240 *
241 * .. External Subroutines ..
242  EXTERNAL cgetrs, cgetrf, clag2z, xerbla, zaxpy, zgemm,
243  \$ zlacpy, zlag2c
244 * ..
245 * .. External Functions ..
246  INTEGER IZAMAX
247  DOUBLE PRECISION DLAMCH, ZLANGE
248  EXTERNAL izamax, dlamch, zlange
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC abs, dble, max, sqrt
252 * ..
253 * .. Statement Functions ..
254  DOUBLE PRECISION CABS1
255 * ..
256 * .. Statement Function definitions ..
257  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
258 * ..
259 * .. Executable Statements ..
260 *
261  info = 0
262  iter = 0
263 *
264 * Test the input parameters.
265 *
266  IF( n.LT.0 ) THEN
267  info = -1
268  ELSE IF( nrhs.LT.0 ) THEN
269  info = -2
270  ELSE IF( lda.LT.max( 1, n ) ) THEN
271  info = -4
272  ELSE IF( ldb.LT.max( 1, n ) ) THEN
273  info = -7
274  ELSE IF( ldx.LT.max( 1, n ) ) THEN
275  info = -9
276  END IF
277  IF( info.NE.0 ) THEN
278  CALL xerbla( 'ZCGESV', -info )
279  RETURN
280  END IF
281 *
282 * Quick return if (N.EQ.0).
283 *
284  IF( n.EQ.0 )
285  \$ RETURN
286 *
287 * Skip single precision iterative refinement if a priori slower
288 * than double precision factorization.
289 *
290  IF( .NOT.doitref ) THEN
291  iter = -1
292  GO TO 40
293  END IF
294 *
295 * Compute some constants.
296 *
297  anrm = zlange( 'I', n, n, a, lda, rwork )
298  eps = dlamch( 'Epsilon' )
299  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
300 *
301 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
302 *
303  ptsa = 1
304  ptsx = ptsa + n*n
305 *
306 * Convert B from double precision to single precision and store the
307 * result in SX.
308 *
309  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
310 *
311  IF( info.NE.0 ) THEN
312  iter = -2
313  GO TO 40
314  END IF
315 *
316 * Convert A from double precision to single precision and store the
317 * result in SA.
318 *
319  CALL zlag2c( n, n, a, lda, swork( ptsa ), n, info )
320 *
321  IF( info.NE.0 ) THEN
322  iter = -2
323  GO TO 40
324  END IF
325 *
326 * Compute the LU factorization of SA.
327 *
328  CALL cgetrf( n, n, swork( ptsa ), n, ipiv, info )
329 *
330  IF( info.NE.0 ) THEN
331  iter = -3
332  GO TO 40
333  END IF
334 *
335 * Solve the system SA*SX = SB.
336 *
337  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
338  \$ swork( ptsx ), n, info )
339 *
340 * Convert SX back to double precision
341 *
342  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
343 *
344 * Compute R = B - AX (R is WORK).
345 *
346  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
347 *
348  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
349  \$ lda, x, ldx, one, work, n )
350 *
351 * Check whether the NRHS normwise backward errors satisfy the
352 * stopping criterion. If yes, set ITER=0 and return.
353 *
354  DO i = 1, nrhs
355  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
356  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
357  IF( rnrm.GT.xnrm*cte )
358  \$ GO TO 10
359  END DO
360 *
361 * If we are here, the NRHS normwise backward errors satisfy the
362 * stopping criterion. We are good to exit.
363 *
364  iter = 0
365  RETURN
366 *
367  10 CONTINUE
368 *
369  DO 30 iiter = 1, itermax
370 *
371 * Convert R (in WORK) from double precision to single precision
372 * and store the result in SX.
373 *
374  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
375 *
376  IF( info.NE.0 ) THEN
377  iter = -2
378  GO TO 40
379  END IF
380 *
381 * Solve the system SA*SX = SR.
382 *
383  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
384  \$ swork( ptsx ), n, info )
385 *
386 * Convert SX back to double precision and update the current
387 * iterate.
388 *
389  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
390 *
391  DO i = 1, nrhs
392  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
393  END DO
394 *
395 * Compute R = B - AX (R is WORK).
396 *
397  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
398 *
399  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
400  \$ a, lda, x, ldx, one, work, n )
401 *
402 * Check whether the NRHS normwise backward errors satisfy the
403 * stopping criterion. If yes, set ITER=IITER>0 and return.
404 *
405  DO i = 1, nrhs
406  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
407  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
408  IF( rnrm.GT.xnrm*cte )
409  \$ GO TO 20
410  END DO
411 *
412 * If we are here, the NRHS normwise backward errors satisfy the
413 * stopping criterion, we are good to exit.
414 *
415  iter = iiter
416 *
417  RETURN
418 *
419  20 CONTINUE
420 *
421  30 CONTINUE
422 *
423 * If we are at this place of the code, this is because we have
424 * performed ITER=ITERMAX iterations and never satisified the stopping
425 * criterion, set up the ITER flag accordingly and follow up on double
426 * precision routine.
427 *
428  iter = -itermax - 1
429 *
430  40 CONTINUE
431 *
432 * Single-precision iterative refinement failed to converge to a
433 * satisfactory solution, so we resort to double precision.
434 *
435  CALL zgetrf( n, n, a, lda, ipiv, info )
436 *
437  IF( info.NE.0 )
438  \$ RETURN
439 *
440  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
441  CALL zgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
442  \$ info )
443 *
444  RETURN
445 *
446 * End of ZCGESV.
447 *
448  END
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
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:109
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:123
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
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:105
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 precisio...
Definition: zcgesv.f:203
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53