--- rpl/lapack/lapack/zcposv.f 2010/12/21 13:48:06 1.4 +++ rpl/lapack/lapack/zcposv.f 2023/08/07 08:39:15 1.19 @@ -1,13 +1,216 @@ - SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, - + SWORK, RWORK, ITER, INFO ) +*> \brief ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * -* -- LAPACK PROTOTYPE driver routine (version 3.3.0) -- +*> \htmlonly +*> Download ZCPOSV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, +* SWORK, RWORK, ITER, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS +* .. +* .. Array Arguments .. +* DOUBLE PRECISION RWORK( * ) +* COMPLEX SWORK( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ), +* $ X( LDX, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> 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. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,NRHS) +*> The N-by-NRHS right hand side matrix B. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (LDX,NRHS) +*> If INFO = 0, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDX +*> \verbatim +*> LDX is INTEGER +*> The leading dimension of the array X. LDX >= max(1,N). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (N,NRHS) +*> This array is used to hold the residual vectors. +*> \endverbatim +*> +*> \param[out] SWORK +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] ITER +*> \verbatim +*> 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 +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> 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. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. * -* November 2010 +*> \ingroup complex16POsolve +* +* ===================================================================== + SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, + $ SWORK, RWORK, ITER, INFO ) * +* -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* .. +* * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS @@ -16,125 +219,10 @@ DOUBLE PRECISION RWORK( * ) COMPLEX SWORK( * ) COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ), - + X( LDX, * ) + $ X( LDX, * ) * .. * -* 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. -* -* Arguments -* ========= -* -* UPLO (input) CHARACTER*1 -* = 'U': Upper triangle of A is stored; -* = 'L': Lower triangle of A is stored. -* -* N (input) INTEGER -* The number of linear equations, i.e., the order of the -* matrix A. N >= 0. -* -* NRHS (input) INTEGER -* The number of right hand sides, i.e., the number of columns -* of the matrix B. NRHS >= 0. -* -* A (input/output) 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. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* B (input) COMPLEX*16 array, dimension (LDB,NRHS) -* The N-by-NRHS right hand side matrix B. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,N). -* -* X (output) COMPLEX*16 array, dimension (LDX,NRHS) -* If INFO = 0, the N-by-NRHS solution matrix X. -* -* LDX (input) INTEGER -* The leading dimension of the array X. LDX >= max(1,N). -* -* WORK (workspace) COMPLEX*16 array, dimension (N*NRHS) -* This array is used to hold the residual vectors. -* -* SWORK (workspace) 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. -* -* RWORK (workspace) DOUBLE PRECISION array, dimension (N) -* -* ITER (output) 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 sucessfully used. -* Returns the number of iterations -* -* INFO (output) 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. -* -* ========= +* ===================================================================== * * .. Parameters .. LOGICAL DOITREF @@ -148,7 +236,7 @@ * COMPLEX*16 NEGONE, ONE PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ), - + ONE = ( 1.0D+00, 0.0D+00 ) ) + $ ONE = ( 1.0D+00, 0.0D+00 ) ) * * .. Local Scalars .. INTEGER I, IITER, PTSA, PTSX @@ -157,7 +245,7 @@ * * .. External Subroutines .. EXTERNAL ZAXPY, ZHEMM, ZLACPY, ZLAT2C, ZLAG2C, CLAG2Z, - + CPOTRF, CPOTRS, XERBLA + $ CPOTRF, CPOTRS, XERBLA, ZPOTRF, ZPOTRS * .. * .. External Functions .. INTEGER IZAMAX @@ -201,7 +289,7 @@ * Quick return if (N.EQ.0). * IF( N.EQ.0 ) - + RETURN + $ RETURN * * Skip single precision iterative refinement if a priori slower * than double precision factorization. @@ -254,7 +342,7 @@ * Solve the system SA*SX = SB. * CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, - + INFO ) + $ INFO ) * * Convert SX back to COMPLEX*16 * @@ -265,7 +353,7 @@ CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N ) * CALL ZHEMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, - + WORK, N ) + $ WORK, N ) * * Check whether the NRHS normwise backward errors satisfy the * stopping criterion. If yes, set ITER=0 and return. @@ -274,7 +362,7 @@ XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) ) RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) ) IF( RNRM.GT.XNRM*CTE ) - + GO TO 10 + $ GO TO 10 END DO * * If we are here, the NRHS normwise backward errors satisfy the @@ -300,7 +388,7 @@ * Solve the system SA*SX = SR. * CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, - + INFO ) + $ INFO ) * * Convert SX back to double precision and update the current * iterate. @@ -316,7 +404,7 @@ CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N ) * CALL ZHEMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, - + WORK, N ) + $ WORK, N ) * * Check whether the NRHS normwise backward errors satisfy the * stopping criterion. If yes, set ITER=IITER>0 and return. @@ -325,7 +413,7 @@ XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) ) RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) ) IF( RNRM.GT.XNRM*CTE ) - + GO TO 20 + $ GO TO 20 END DO * * If we are here, the NRHS normwise backward errors satisfy the @@ -340,7 +428,7 @@ 30 CONTINUE * * If we are at this place of the code, this is because we have -* performed ITER=ITERMAX iterations and never satisified the +* performed ITER=ITERMAX iterations and never satisfied the * stopping criterion, set up the ITER flag accordingly and follow * up on double precision routine. * @@ -354,13 +442,13 @@ CALL ZPOTRF( UPLO, N, A, LDA, INFO ) * IF( INFO.NE.0 ) - + RETURN + $ RETURN * CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX ) CALL ZPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO ) * RETURN * -* End of ZCPOSV. +* End of ZCPOSV * END