sppequ.c
Go to the documentation of this file.
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_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:13