 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ zcposv()

 subroutine zcposv ( character UPLO, integer N, integer NRHS, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldx, * ) X, integer LDX, complex*16, dimension( n, * ) WORK, complex, dimension( * ) SWORK, double precision, dimension( * ) RWORK, integer ITER, integer INFO )

ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices

Purpose:
``` ZCPOSV computes the solution to a complex system of linear equations
A * X = B,
where A is an N-by-N Hermitian positive definite matrix and X and B
are N-by-NRHS matrices.

ZCPOSV first attempts to factorize the matrix in COMPLEX and use this
factorization within an iterative refinement procedure to produce a
solution with COMPLEX*16 normwise backward error quality (see below).
If the approach fails the method switches to a COMPLEX*16
factorization and solve.

The iterative refinement is not going to be a winning strategy if
the ratio COMPLEX performance over COMPLEX*16 performance is too
small. A reasonable strategy should take the number of right-hand
sides and the size of the matrix into account. This might be done
with a call to ILAENV in the future. Up to now, we always try
iterative refinement.

The iterative refinement process is stopped if
ITER > ITERMAX
or for all the RHS we have:
RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
where
o ITER is the number of the current iteration in the iterative
refinement process
o RNRM is the infinity-norm of the residual
o XNRM is the infinity-norm of the solution
o ANRM is the infinity-operator-norm of the matrix A
o EPS is the machine epsilon returned by DLAMCH('Epsilon')
The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
respectively.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.``` [in,out] A ``` A is COMPLEX*16 array, dimension (LDA,N) On entry, the Hermitian matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. Note that the imaginary parts of the diagonal elements need not be set and are assumed to be zero. On exit, if iterative refinement has been successfully used (INFO = 0 and ITER >= 0, see description below), then A is unchanged, if double precision factorization has been used (INFO = 0 and ITER < 0, see description below), then the array A contains the factor U or L from the Cholesky factorization A = U**H*U or A = L*L**H.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in] B ``` B is COMPLEX*16 array, dimension (LDB,NRHS) The N-by-NRHS right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [out] X ``` X is COMPLEX*16 array, dimension (LDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] WORK ``` WORK is COMPLEX*16 array, dimension (N,NRHS) This array is used to hold the residual vectors.``` [out] SWORK ``` SWORK is COMPLEX array, dimension (N*(N+NRHS)) This array is used to use the single precision matrix and the right-hand sides or solutions in single precision.``` [out] RWORK ` RWORK is DOUBLE PRECISION array, dimension (N)` [out] ITER ``` ITER is INTEGER < 0: iterative refinement has failed, COMPLEX*16 factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of CPOTRF -31: stop the iterative refinement after the 30th iterations > 0: iterative refinement has been successfully used. Returns the number of iterations``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i of (COMPLEX*16) A is not positive definite, so the factorization could not be completed, and the solution has not been computed.```

Definition at line 207 of file zcposv.f.

209 *
210 * -- LAPACK driver routine --
211 * -- LAPACK is a software package provided by Univ. of Tennessee, --
212 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 *
214 * .. Scalar Arguments ..
215  CHARACTER UPLO
216  INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
217 * ..
218 * .. Array Arguments ..
219  DOUBLE PRECISION RWORK( * )
220  COMPLEX SWORK( * )
221  COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
222  \$ X( LDX, * )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228  LOGICAL DOITREF
229  parameter( doitref = .true. )
230 *
231  INTEGER ITERMAX
232  parameter( itermax = 30 )
233 *
234  DOUBLE PRECISION BWDMAX
235  parameter( bwdmax = 1.0e+00 )
236 *
237  COMPLEX*16 NEGONE, ONE
238  parameter( negone = ( -1.0d+00, 0.0d+00 ),
239  \$ one = ( 1.0d+00, 0.0d+00 ) )
240 *
241 * .. Local Scalars ..
242  INTEGER I, IITER, PTSA, PTSX
243  DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
244  COMPLEX*16 ZDUM
245 *
246 * .. External Subroutines ..
247  EXTERNAL zaxpy, zhemm, zlacpy, zlat2c, zlag2c, clag2z,
249 * ..
250 * .. External Functions ..
251  INTEGER IZAMAX
252  DOUBLE PRECISION DLAMCH, ZLANHE
253  LOGICAL LSAME
254  EXTERNAL izamax, dlamch, zlanhe, lsame
255 * ..
256 * .. Intrinsic Functions ..
257  INTRINSIC abs, dble, max, sqrt
258 * .. Statement Functions ..
259  DOUBLE PRECISION CABS1
260 * ..
261 * .. Statement Function definitions ..
262  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
263 * ..
264 * .. Executable Statements ..
265 *
266  info = 0
267  iter = 0
268 *
269 * Test the input parameters.
270 *
271  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
272  info = -1
273  ELSE IF( n.LT.0 ) THEN
274  info = -2
275  ELSE IF( nrhs.LT.0 ) THEN
276  info = -3
277  ELSE IF( lda.LT.max( 1, n ) ) THEN
278  info = -5
279  ELSE IF( ldb.LT.max( 1, n ) ) THEN
280  info = -7
281  ELSE IF( ldx.LT.max( 1, n ) ) THEN
282  info = -9
283  END IF
284  IF( info.NE.0 ) THEN
285  CALL xerbla( 'ZCPOSV', -info )
286  RETURN
287  END IF
288 *
289 * Quick return if (N.EQ.0).
290 *
291  IF( n.EQ.0 )
292  \$ RETURN
293 *
294 * Skip single precision iterative refinement if a priori slower
295 * than double precision factorization.
296 *
297  IF( .NOT.doitref ) THEN
298  iter = -1
299  GO TO 40
300  END IF
301 *
302 * Compute some constants.
303 *
304  anrm = zlanhe( 'I', uplo, n, a, lda, rwork )
305  eps = dlamch( 'Epsilon' )
306  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
307 *
308 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
309 *
310  ptsa = 1
311  ptsx = ptsa + n*n
312 *
313 * Convert B from double precision to single precision and store the
314 * result in SX.
315 *
316  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
317 *
318  IF( info.NE.0 ) THEN
319  iter = -2
320  GO TO 40
321  END IF
322 *
323 * Convert A from double precision to single precision and store the
324 * result in SA.
325 *
326  CALL zlat2c( uplo, n, a, lda, swork( ptsa ), n, info )
327 *
328  IF( info.NE.0 ) THEN
329  iter = -2
330  GO TO 40
331  END IF
332 *
333 * Compute the Cholesky factorization of SA.
334 *
335  CALL cpotrf( uplo, n, swork( ptsa ), n, info )
336 *
337  IF( info.NE.0 ) THEN
338  iter = -3
339  GO TO 40
340  END IF
341 *
342 * Solve the system SA*SX = SB.
343 *
344  CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
345  \$ info )
346 *
347 * Convert SX back to COMPLEX*16
348 *
349  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
350 *
351 * Compute R = B - AX (R is WORK).
352 *
353  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
354 *
355  CALL zhemm( 'Left', uplo, n, nrhs, negone, a, lda, x, ldx, one,
356  \$ work, n )
357 *
358 * Check whether the NRHS normwise backward errors satisfy the
359 * stopping criterion. If yes, set ITER=0 and return.
360 *
361  DO i = 1, nrhs
362  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
363  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
364  IF( rnrm.GT.xnrm*cte )
365  \$ GO TO 10
366  END DO
367 *
368 * If we are here, the NRHS normwise backward errors satisfy the
369 * stopping criterion. We are good to exit.
370 *
371  iter = 0
372  RETURN
373 *
374  10 CONTINUE
375 *
376  DO 30 iiter = 1, itermax
377 *
378 * Convert R (in WORK) from double precision to single precision
379 * and store the result in SX.
380 *
381  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
382 *
383  IF( info.NE.0 ) THEN
384  iter = -2
385  GO TO 40
386  END IF
387 *
388 * Solve the system SA*SX = SR.
389 *
390  CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
391  \$ info )
392 *
393 * Convert SX back to double precision and update the current
394 * iterate.
395 *
396  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
397 *
398  DO i = 1, nrhs
399  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
400  END DO
401 *
402 * Compute R = B - AX (R is WORK).
403 *
404  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
405 *
406  CALL zhemm( 'L', uplo, n, nrhs, negone, a, lda, x, ldx, one,
407  \$ work, n )
408 *
409 * Check whether the NRHS normwise backward errors satisfy the
410 * stopping criterion. If yes, set ITER=IITER>0 and return.
411 *
412  DO i = 1, nrhs
413  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
414  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
415  IF( rnrm.GT.xnrm*cte )
416  \$ GO TO 20
417  END DO
418 *
419 * If we are here, the NRHS normwise backward errors satisfy the
420 * stopping criterion, we are good to exit.
421 *
422  iter = iiter
423 *
424  RETURN
425 *
426  20 CONTINUE
427 *
428  30 CONTINUE
429 *
430 * If we are at this place of the code, this is because we have
431 * performed ITER=ITERMAX iterations and never satisfied the
432 * stopping criterion, set up the ITER flag accordingly and follow
433 * up on double precision routine.
434 *
435  iter = -itermax - 1
436 *
437  40 CONTINUE
438 *
439 * Single-precision iterative refinement failed to converge to a
440 * satisfactory solution, so we resort to double precision.
441 *
442  CALL zpotrf( uplo, n, a, lda, info )
443 *
444  IF( info.NE.0 )
445  \$ RETURN
446 *
447  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
448  CALL zpotrs( uplo, n, nrhs, a, lda, x, ldx, info )
449 *
450  RETURN
451 *
452 * End of ZCPOSV
453 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zhemm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZHEMM
Definition: zhemm.f:191
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: zlanhe.f:124
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 zlat2c(UPLO, N, A, LDA, SA, LDSA, INFO)
ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.
Definition: zlat2c.f:111
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 zpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
ZPOTRS
Definition: zpotrs.f:110
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:107
subroutine cpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
CPOTRS
Definition: cpotrs.f:110
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
Here is the call graph for this function:
Here is the caller graph for this function: