slaqsp.c
Go to the documentation of this file.
00001 /* slaqsp.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 slaqsp_(char *uplo, integer *n, real *ap, real *s, real *
00017         scond, real *amax, char *equed)
00018 {
00019     /* System generated locals */
00020     integer i__1, i__2;
00021 
00022     /* Local variables */
00023     integer i__, j, jc;
00024     real cj, large;
00025     extern logical lsame_(char *, char *);
00026     real small;
00027     extern doublereal slamch_(char *);
00028 
00029 
00030 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00031 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00032 /*     November 2006 */
00033 
00034 /*     .. Scalar Arguments .. */
00035 /*     .. */
00036 /*     .. Array Arguments .. */
00037 /*     .. */
00038 
00039 /*  Purpose */
00040 /*  ======= */
00041 
00042 /*  SLAQSP equilibrates a symmetric matrix A using the scaling factors */
00043 /*  in the vector S. */
00044 
00045 /*  Arguments */
00046 /*  ========= */
00047 
00048 /*  UPLO    (input) CHARACTER*1 */
00049 /*          Specifies whether the upper or lower triangular part of the */
00050 /*          symmetric matrix A is stored. */
00051 /*          = 'U':  Upper triangular */
00052 /*          = 'L':  Lower triangular */
00053 
00054 /*  N       (input) INTEGER */
00055 /*          The order of the matrix A.  N >= 0. */
00056 
00057 /*  AP      (input/output) REAL array, dimension (N*(N+1)/2) */
00058 /*          On entry, the upper or lower triangle of the symmetric matrix */
00059 /*          A, packed columnwise in a linear array.  The j-th column of A */
00060 /*          is stored in the array AP as follows: */
00061 /*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
00062 /*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
00063 
00064 /*          On exit, the equilibrated matrix:  diag(S) * A * diag(S), in */
00065 /*          the same storage format as A. */
00066 
00067 /*  S       (input) REAL array, dimension (N) */
00068 /*          The scale factors for A. */
00069 
00070 /*  SCOND   (input) REAL */
00071 /*          Ratio of the smallest S(i) to the largest S(i). */
00072 
00073 /*  AMAX    (input) REAL */
00074 /*          Absolute value of largest matrix entry. */
00075 
00076 /*  EQUED   (output) CHARACTER*1 */
00077 /*          Specifies whether or not equilibration was done. */
00078 /*          = 'N':  No equilibration. */
00079 /*          = 'Y':  Equilibration was done, i.e., A has been replaced by */
00080 /*                  diag(S) * A * diag(S). */
00081 
00082 /*  Internal Parameters */
00083 /*  =================== */
00084 
00085 /*  THRESH is a threshold value used to decide if scaling should be done */
00086 /*  based on the ratio of the scaling factors.  If SCOND < THRESH, */
00087 /*  scaling is done. */
00088 
00089 /*  LARGE and SMALL are threshold values used to decide if scaling should */
00090 /*  be done based on the absolute size of the largest matrix element. */
00091 /*  If AMAX > LARGE or AMAX < SMALL, scaling is done. */
00092 
00093 /*  ===================================================================== */
00094 
00095 /*     .. Parameters .. */
00096 /*     .. */
00097 /*     .. Local Scalars .. */
00098 /*     .. */
00099 /*     .. External Functions .. */
00100 /*     .. */
00101 /*     .. Executable Statements .. */
00102 
00103 /*     Quick return if possible */
00104 
00105     /* Parameter adjustments */
00106     --s;
00107     --ap;
00108 
00109     /* Function Body */
00110     if (*n <= 0) {
00111         *(unsigned char *)equed = 'N';
00112         return 0;
00113     }
00114 
00115 /*     Initialize LARGE and SMALL. */
00116 
00117     small = slamch_("Safe minimum") / slamch_("Precision");
00118     large = 1.f / small;
00119 
00120     if (*scond >= .1f && *amax >= small && *amax <= large) {
00121 
00122 /*        No equilibration */
00123 
00124         *(unsigned char *)equed = 'N';
00125     } else {
00126 
00127 /*        Replace A by diag(S) * A * diag(S). */
00128 
00129         if (lsame_(uplo, "U")) {
00130 
00131 /*           Upper triangle of A is stored. */
00132 
00133             jc = 1;
00134             i__1 = *n;
00135             for (j = 1; j <= i__1; ++j) {
00136                 cj = s[j];
00137                 i__2 = j;
00138                 for (i__ = 1; i__ <= i__2; ++i__) {
00139                     ap[jc + i__ - 1] = cj * s[i__] * ap[jc + i__ - 1];
00140 /* L10: */
00141                 }
00142                 jc += j;
00143 /* L20: */
00144             }
00145         } else {
00146 
00147 /*           Lower triangle of A is stored. */
00148 
00149             jc = 1;
00150             i__1 = *n;
00151             for (j = 1; j <= i__1; ++j) {
00152                 cj = s[j];
00153                 i__2 = *n;
00154                 for (i__ = j; i__ <= i__2; ++i__) {
00155                     ap[jc + i__ - j] = cj * s[i__] * ap[jc + i__ - j];
00156 /* L30: */
00157                 }
00158                 jc = jc + *n - j + 1;
00159 /* L40: */
00160             }
00161         }
00162         *(unsigned char *)equed = 'Y';
00163     }
00164 
00165     return 0;
00166 
00167 /*     End of SLAQSP */
00168 
00169 } /* slaqsp_ */


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