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