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

Download ZCPOSV + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: