00001 SUBROUTINE ZLAQHP( UPLO, N, AP, S, SCOND, AMAX, EQUED ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER EQUED, UPLO 00010 INTEGER N 00011 DOUBLE PRECISION AMAX, SCOND 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION S( * ) 00015 COMPLEX*16 AP( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZLAQHP equilibrates a Hermitian matrix A using the scaling factors 00022 * in the vector S. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * Specifies whether the upper or lower triangular part of the 00029 * Hermitian matrix A is stored. 00030 * = 'U': Upper triangular 00031 * = 'L': Lower triangular 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. N >= 0. 00035 * 00036 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 00037 * On entry, the upper or lower triangle of the Hermitian matrix 00038 * A, packed columnwise in a linear array. The j-th column of A 00039 * is stored in the array AP as follows: 00040 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00041 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00042 * 00043 * On exit, the equilibrated matrix: diag(S) * A * diag(S), in 00044 * the same storage format as A. 00045 * 00046 * S (input) DOUBLE PRECISION array, dimension (N) 00047 * The scale factors for A. 00048 * 00049 * SCOND (input) DOUBLE PRECISION 00050 * Ratio of the smallest S(i) to the largest S(i). 00051 * 00052 * AMAX (input) DOUBLE PRECISION 00053 * Absolute value of largest matrix entry. 00054 * 00055 * EQUED (output) CHARACTER*1 00056 * Specifies whether or not equilibration was done. 00057 * = 'N': No equilibration. 00058 * = 'Y': Equilibration was done, i.e., A has been replaced by 00059 * diag(S) * A * diag(S). 00060 * 00061 * Internal Parameters 00062 * =================== 00063 * 00064 * THRESH is a threshold value used to decide if scaling should be done 00065 * based on the ratio of the scaling factors. If SCOND < THRESH, 00066 * scaling is done. 00067 * 00068 * LARGE and SMALL are threshold values used to decide if scaling should 00069 * be done based on the absolute size of the largest matrix element. 00070 * If AMAX > LARGE or AMAX < SMALL, scaling is done. 00071 * 00072 * ===================================================================== 00073 * 00074 * .. Parameters .. 00075 DOUBLE PRECISION ONE, THRESH 00076 PARAMETER ( ONE = 1.0D+0, THRESH = 0.1D+0 ) 00077 * .. 00078 * .. Local Scalars .. 00079 INTEGER I, J, JC 00080 DOUBLE PRECISION CJ, LARGE, SMALL 00081 * .. 00082 * .. External Functions .. 00083 LOGICAL LSAME 00084 DOUBLE PRECISION DLAMCH 00085 EXTERNAL LSAME, DLAMCH 00086 * .. 00087 * .. Intrinsic Functions .. 00088 INTRINSIC DBLE 00089 * .. 00090 * .. Executable Statements .. 00091 * 00092 * Quick return if possible 00093 * 00094 IF( N.LE.0 ) THEN 00095 EQUED = 'N' 00096 RETURN 00097 END IF 00098 * 00099 * Initialize LARGE and SMALL. 00100 * 00101 SMALL = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) 00102 LARGE = ONE / SMALL 00103 * 00104 IF( SCOND.GE.THRESH .AND. AMAX.GE.SMALL .AND. AMAX.LE.LARGE ) THEN 00105 * 00106 * No equilibration 00107 * 00108 EQUED = 'N' 00109 ELSE 00110 * 00111 * Replace A by diag(S) * A * diag(S). 00112 * 00113 IF( LSAME( UPLO, 'U' ) ) THEN 00114 * 00115 * Upper triangle of A is stored. 00116 * 00117 JC = 1 00118 DO 20 J = 1, N 00119 CJ = S( J ) 00120 DO 10 I = 1, J - 1 00121 AP( JC+I-1 ) = CJ*S( I )*AP( JC+I-1 ) 00122 10 CONTINUE 00123 AP( JC+J-1 ) = CJ*CJ*DBLE( AP( JC+J-1 ) ) 00124 JC = JC + J 00125 20 CONTINUE 00126 ELSE 00127 * 00128 * Lower triangle of A is stored. 00129 * 00130 JC = 1 00131 DO 40 J = 1, N 00132 CJ = S( J ) 00133 AP( JC ) = CJ*CJ*DBLE( AP( JC ) ) 00134 DO 30 I = J + 1, N 00135 AP( JC+I-J ) = CJ*S( I )*AP( JC+I-J ) 00136 30 CONTINUE 00137 JC = JC + N - J + 1 00138 40 CONTINUE 00139 END IF 00140 EQUED = 'Y' 00141 END IF 00142 * 00143 RETURN 00144 * 00145 * End of ZLAQHP 00146 * 00147 END
1.7.2