claqsb.c
Go to the documentation of this file.
00001 /* claqsb.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 claqsb_(char *uplo, integer *n, integer *kd, complex *ab, 
00017          integer *ldab, real *s, real *scond, real *amax, char *equed)
00018 {
00019     /* System generated locals */
00020     integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4;
00021     real r__1;
00022     complex q__1;
00023 
00024     /* Local variables */
00025     integer i__, j;
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 /*  CLAQSB equilibrates a symmetric band matrix A using the scaling */
00045 /*  factors 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 /*  KD      (input) INTEGER */
00060 /*          The number of super-diagonals of the matrix A if UPLO = 'U', */
00061 /*          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0. */
00062 
00063 /*  AB      (input/output) COMPLEX array, dimension (LDAB,N) */
00064 /*          On entry, the upper or lower triangle of the symmetric band */
00065 /*          matrix A, stored in the first KD+1 rows of the array.  The */
00066 /*          j-th column of A is stored in the j-th column of the array AB */
00067 /*          as follows: */
00068 /*          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
00069 /*          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd). */
00070 
00071 /*          On exit, if INFO = 0, the triangular factor U or L from the */
00072 /*          Cholesky factorization A = U'*U or A = L*L' of the band */
00073 /*          matrix A, in the same storage format as A. */
00074 
00075 /*  LDAB    (input) INTEGER */
00076 /*          The leading dimension of the array AB.  LDAB >= KD+1. */
00077 
00078 /*  S       (input) REAL array, dimension (N) */
00079 /*          The scale factors for A. */
00080 
00081 /*  SCOND   (input) REAL */
00082 /*          Ratio of the smallest S(i) to the largest S(i). */
00083 
00084 /*  AMAX    (input) REAL */
00085 /*          Absolute value of largest matrix entry. */
00086 
00087 /*  EQUED   (output) CHARACTER*1 */
00088 /*          Specifies whether or not equilibration was done. */
00089 /*          = 'N':  No equilibration. */
00090 /*          = 'Y':  Equilibration was done, i.e., A has been replaced by */
00091 /*                  diag(S) * A * diag(S). */
00092 
00093 /*  Internal Parameters */
00094 /*  =================== */
00095 
00096 /*  THRESH is a threshold value used to decide if scaling should be done */
00097 /*  based on the ratio of the scaling factors.  If SCOND < THRESH, */
00098 /*  scaling is done. */
00099 
00100 /*  LARGE and SMALL are threshold values used to decide if scaling should */
00101 /*  be done based on the absolute size of the largest matrix element. */
00102 /*  If AMAX > LARGE or AMAX < SMALL, scaling is done. */
00103 
00104 /*  ===================================================================== */
00105 
00106 /*     .. Parameters .. */
00107 /*     .. */
00108 /*     .. Local Scalars .. */
00109 /*     .. */
00110 /*     .. External Functions .. */
00111 /*     .. */
00112 /*     .. Intrinsic Functions .. */
00113 /*     .. */
00114 /*     .. Executable Statements .. */
00115 
00116 /*     Quick return if possible */
00117 
00118     /* Parameter adjustments */
00119     ab_dim1 = *ldab;
00120     ab_offset = 1 + ab_dim1;
00121     ab -= ab_offset;
00122     --s;
00123 
00124     /* Function Body */
00125     if (*n <= 0) {
00126         *(unsigned char *)equed = 'N';
00127         return 0;
00128     }
00129 
00130 /*     Initialize LARGE and SMALL. */
00131 
00132     small = slamch_("Safe minimum") / slamch_("Precision");
00133     large = 1.f / small;
00134 
00135     if (*scond >= .1f && *amax >= small && *amax <= large) {
00136 
00137 /*        No equilibration */
00138 
00139         *(unsigned char *)equed = 'N';
00140     } else {
00141 
00142 /*        Replace A by diag(S) * A * diag(S). */
00143 
00144         if (lsame_(uplo, "U")) {
00145 
00146 /*           Upper triangle of A is stored in band format. */
00147 
00148             i__1 = *n;
00149             for (j = 1; j <= i__1; ++j) {
00150                 cj = s[j];
00151 /* Computing MAX */
00152                 i__2 = 1, i__3 = j - *kd;
00153                 i__4 = j;
00154                 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
00155                     i__2 = *kd + 1 + i__ - j + j * ab_dim1;
00156                     r__1 = cj * s[i__];
00157                     i__3 = *kd + 1 + i__ - j + j * ab_dim1;
00158                     q__1.r = r__1 * ab[i__3].r, q__1.i = r__1 * ab[i__3].i;
00159                     ab[i__2].r = q__1.r, ab[i__2].i = q__1.i;
00160 /* L10: */
00161                 }
00162 /* L20: */
00163             }
00164         } else {
00165 
00166 /*           Lower triangle of A is stored. */
00167 
00168             i__1 = *n;
00169             for (j = 1; j <= i__1; ++j) {
00170                 cj = s[j];
00171 /* Computing MIN */
00172                 i__2 = *n, i__3 = j + *kd;
00173                 i__4 = min(i__2,i__3);
00174                 for (i__ = j; i__ <= i__4; ++i__) {
00175                     i__2 = i__ + 1 - j + j * ab_dim1;
00176                     r__1 = cj * s[i__];
00177                     i__3 = i__ + 1 - j + j * ab_dim1;
00178                     q__1.r = r__1 * ab[i__3].r, q__1.i = r__1 * ab[i__3].i;
00179                     ab[i__2].r = q__1.r, ab[i__2].i = q__1.i;
00180 /* L30: */
00181                 }
00182 /* L40: */
00183             }
00184         }
00185         *(unsigned char *)equed = 'Y';
00186     }
00187 
00188     return 0;
00189 
00190 /*     End of CLAQSB */
00191 
00192 } /* claqsb_ */


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