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


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