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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:31