00001 /* sppequ.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int sppequ_(char *uplo, integer *n, real *ap, real *s, real * 00017 scond, real *amax, integer *info) 00018 { 00019 /* System generated locals */ 00020 integer i__1; 00021 real r__1, r__2; 00022 00023 /* Builtin functions */ 00024 double sqrt(doublereal); 00025 00026 /* Local variables */ 00027 integer i__, jj; 00028 real smin; 00029 extern logical lsame_(char *, char *); 00030 logical upper; 00031 extern /* Subroutine */ int xerbla_(char *, integer *); 00032 00033 00034 /* -- LAPACK routine (version 3.2) -- */ 00035 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00036 /* November 2006 */ 00037 00038 /* .. Scalar Arguments .. */ 00039 /* .. */ 00040 /* .. Array Arguments .. */ 00041 /* .. */ 00042 00043 /* Purpose */ 00044 /* ======= */ 00045 00046 /* SPPEQU computes row and column scalings intended to equilibrate a */ 00047 /* symmetric positive definite matrix A in packed storage and reduce */ 00048 /* its condition number (with respect to the two-norm). S contains the */ 00049 /* scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix */ 00050 /* B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal. */ 00051 /* This choice of S puts the condition number of B within a factor N of */ 00052 /* the smallest possible condition number over all possible diagonal */ 00053 /* scalings. */ 00054 00055 /* Arguments */ 00056 /* ========= */ 00057 00058 /* UPLO (input) CHARACTER*1 */ 00059 /* = 'U': Upper triangle of A is stored; */ 00060 /* = 'L': Lower triangle of A is stored. */ 00061 00062 /* N (input) INTEGER */ 00063 /* The order of the matrix A. N >= 0. */ 00064 00065 /* AP (input) REAL array, dimension (N*(N+1)/2) */ 00066 /* The upper or lower triangle of the symmetric matrix A, packed */ 00067 /* columnwise in a linear array. The j-th column of A is stored */ 00068 /* in the array AP as follows: */ 00069 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */ 00070 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */ 00071 00072 /* S (output) REAL array, dimension (N) */ 00073 /* If INFO = 0, S contains the scale factors for A. */ 00074 00075 /* SCOND (output) REAL */ 00076 /* If INFO = 0, S contains the ratio of the smallest S(i) to */ 00077 /* the largest S(i). If SCOND >= 0.1 and AMAX is neither too */ 00078 /* large nor too small, it is not worth scaling by S. */ 00079 00080 /* AMAX (output) REAL */ 00081 /* Absolute value of largest matrix element. If AMAX is very */ 00082 /* close to overflow or very close to underflow, the matrix */ 00083 /* should be scaled. */ 00084 00085 /* INFO (output) INTEGER */ 00086 /* = 0: successful exit */ 00087 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00088 /* > 0: if INFO = i, the i-th diagonal element is nonpositive. */ 00089 00090 /* ===================================================================== */ 00091 00092 /* .. Parameters .. */ 00093 /* .. */ 00094 /* .. Local Scalars .. */ 00095 /* .. */ 00096 /* .. External Functions .. */ 00097 /* .. */ 00098 /* .. External Subroutines .. */ 00099 /* .. */ 00100 /* .. Intrinsic Functions .. */ 00101 /* .. */ 00102 /* .. Executable Statements .. */ 00103 00104 /* Test the input parameters. */ 00105 00106 /* Parameter adjustments */ 00107 --s; 00108 --ap; 00109 00110 /* Function Body */ 00111 *info = 0; 00112 upper = lsame_(uplo, "U"); 00113 if (! upper && ! lsame_(uplo, "L")) { 00114 *info = -1; 00115 } else if (*n < 0) { 00116 *info = -2; 00117 } 00118 if (*info != 0) { 00119 i__1 = -(*info); 00120 xerbla_("SPPEQU", &i__1); 00121 return 0; 00122 } 00123 00124 /* Quick return if possible */ 00125 00126 if (*n == 0) { 00127 *scond = 1.f; 00128 *amax = 0.f; 00129 return 0; 00130 } 00131 00132 /* Initialize SMIN and AMAX. */ 00133 00134 s[1] = ap[1]; 00135 smin = s[1]; 00136 *amax = s[1]; 00137 00138 if (upper) { 00139 00140 /* UPLO = 'U': Upper triangle of A is stored. */ 00141 /* Find the minimum and maximum diagonal elements. */ 00142 00143 jj = 1; 00144 i__1 = *n; 00145 for (i__ = 2; i__ <= i__1; ++i__) { 00146 jj += i__; 00147 s[i__] = ap[jj]; 00148 /* Computing MIN */ 00149 r__1 = smin, r__2 = s[i__]; 00150 smin = dmin(r__1,r__2); 00151 /* Computing MAX */ 00152 r__1 = *amax, r__2 = s[i__]; 00153 *amax = dmax(r__1,r__2); 00154 /* L10: */ 00155 } 00156 00157 } else { 00158 00159 /* UPLO = 'L': Lower triangle of A is stored. */ 00160 /* Find the minimum and maximum diagonal elements. */ 00161 00162 jj = 1; 00163 i__1 = *n; 00164 for (i__ = 2; i__ <= i__1; ++i__) { 00165 jj = jj + *n - i__ + 2; 00166 s[i__] = ap[jj]; 00167 /* Computing MIN */ 00168 r__1 = smin, r__2 = s[i__]; 00169 smin = dmin(r__1,r__2); 00170 /* Computing MAX */ 00171 r__1 = *amax, r__2 = s[i__]; 00172 *amax = dmax(r__1,r__2); 00173 /* L20: */ 00174 } 00175 } 00176 00177 if (smin <= 0.f) { 00178 00179 /* Find the first non-positive diagonal element and return. */ 00180 00181 i__1 = *n; 00182 for (i__ = 1; i__ <= i__1; ++i__) { 00183 if (s[i__] <= 0.f) { 00184 *info = i__; 00185 return 0; 00186 } 00187 /* L30: */ 00188 } 00189 } else { 00190 00191 /* Set the scale factors to the reciprocals */ 00192 /* of the diagonal elements. */ 00193 00194 i__1 = *n; 00195 for (i__ = 1; i__ <= i__1; ++i__) { 00196 s[i__] = 1.f / sqrt(s[i__]); 00197 /* L40: */ 00198 } 00199 00200 /* Compute SCOND = min(S(I)) / max(S(I)) */ 00201 00202 *scond = sqrt(smin) / sqrt(*amax); 00203 } 00204 return 0; 00205 00206 /* End of SPPEQU */ 00207 00208 } /* sppequ_ */