LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 principal minor of order i
                of (COMPLEX*16) A is not positive, 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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 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 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 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
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zpotrf(uplo, n, a, lda, info)
ZPOTRF
Definition zpotrf.f:107
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 zpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
ZPOTRS
Definition zpotrs.f:110
Here is the call graph for this function:
Here is the caller graph for this function: