LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.EQ.0 and ITER.GE.0, see description below), then A is
          unchanged, if double precision factorization has been used
          (INFO.EQ.0 and ITER.LT.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.
Date
June 2016

Definition at line 211 of file zcposv.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: