git.maemo.org Git - opencv/blob - 3rdparty/lapack/slascl.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slascl_(char *type__, integer *kl, integer *ku, real *
   4         cfrom, real *cto, integer *m, integer *n, real *a, integer *lda,
   5         integer *info)
   6 {
   7     /* System generated locals */
   8     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
   9
  10     /* Local variables */
  11     integer i__, j, k1, k2, k3, k4;
  12     real mul, cto1;
  13     logical done;
  14     real ctoc;
  15     extern logical lsame_(char *, char *);
  16     integer itype;
  17     real cfrom1;
  18     extern doublereal slamch_(char *);
  19     real cfromc;
  20     extern /* Subroutine */ int xerbla_(char *, integer *);
  21     real bignum, smlnum;
  22
  23
  24 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  25 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  26 /*     November 2006 */
  27
  28 /*     .. Scalar Arguments .. */
  29 /*     .. */
  30 /*     .. Array Arguments .. */
  31 /*     .. */
  32
  33 /*  Purpose */
  34 /*  ======= */
  35
  36 /*  SLASCL multiplies the M by N real matrix A by the real scalar */
  37 /*  CTO/CFROM.  This is done without over/underflow as long as the final */
  38 /*  result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that */
  39 /*  A may be full, upper triangular, lower triangular, upper Hessenberg, */
  40 /*  or banded. */
  41
  42 /*  Arguments */
  43 /*  ========= */
  44
  45 /*  TYPE    (input) CHARACTER*1 */
  46 /*          TYPE indices the storage type of the input matrix. */
  47 /*          = 'G':  A is a full matrix. */
  48 /*          = 'L':  A is a lower triangular matrix. */
  49 /*          = 'U':  A is an upper triangular matrix. */
  50 /*          = 'H':  A is an upper Hessenberg matrix. */
  51 /*          = 'B':  A is a symmetric band matrix with lower bandwidth KL */
  52 /*                  and upper bandwidth KU and with the only the lower */
  53 /*                  half stored. */
  54 /*          = 'Q':  A is a symmetric band matrix with lower bandwidth KL */
  55 /*                  and upper bandwidth KU and with the only the upper */
  56 /*                  half stored. */
  57 /*          = 'Z':  A is a band matrix with lower bandwidth KL and upper */
  58 /*                  bandwidth KU. */
  59
  60 /*  KL      (input) INTEGER */
  61 /*          The lower bandwidth of A.  Referenced only if TYPE = 'B', */
  62 /*          'Q' or 'Z'. */
  63
  64 /*  KU      (input) INTEGER */
  65 /*          The upper bandwidth of A.  Referenced only if TYPE = 'B', */
  66 /*          'Q' or 'Z'. */
  67
  68 /*  CFROM   (input) REAL */
  69 /*  CTO     (input) REAL */
  70 /*          The matrix A is multiplied by CTO/CFROM. A(I,J) is computed */
  71 /*          without over/underflow if the final result CTO*A(I,J)/CFROM */
  72 /*          can be represented without over/underflow.  CFROM must be */
  73 /*          nonzero. */
  74
  75 /*  M       (input) INTEGER */
  76 /*          The number of rows of the matrix A.  M >= 0. */
  77
  78 /*  N       (input) INTEGER */
  79 /*          The number of columns of the matrix A.  N >= 0. */
  80
  81 /*  A       (input/output) REAL array, dimension (LDA,N) */
  82 /*          The matrix to be multiplied by CTO/CFROM.  See TYPE for the */
  83 /*          storage type. */
  84
  85 /*  LDA     (input) INTEGER */
  86 /*          The leading dimension of the array A.  LDA >= max(1,M). */
  87
  88 /*  INFO    (output) INTEGER */
  89 /*          0  - successful exit */
  90 /*          <0 - if INFO = -i, the i-th argument had an illegal value. */
  91
  92 /*  ===================================================================== */
  93
  94 /*     .. Parameters .. */
  95 /*     .. */
  96 /*     .. Local Scalars .. */
  97 /*     .. */
  98 /*     .. External Functions .. */
  99 /*     .. */
 100 /*     .. Intrinsic Functions .. */
 101 /*     .. */
 102 /*     .. External Subroutines .. */
 103 /*     .. */
 104 /*     .. Executable Statements .. */
 105
 106 /*     Test the input arguments */
 107
 108     /* Parameter adjustments */
 109     a_dim1 = *lda;
 110     a_offset = 1 + a_dim1;
 111     a -= a_offset;
 112
 113     /* Function Body */
 114     *info = 0;
 115
 116     if (lsame_(type__, "G")) {
 117         itype = 0;
 118     } else if (lsame_(type__, "L")) {
 119         itype = 1;
 120     } else if (lsame_(type__, "U")) {
 121         itype = 2;
 122     } else if (lsame_(type__, "H")) {
 123         itype = 3;
 124     } else if (lsame_(type__, "B")) {
 125         itype = 4;
 126     } else if (lsame_(type__, "Q")) {
 127         itype = 5;
 128     } else if (lsame_(type__, "Z")) {
 129         itype = 6;
 130     } else {
 131         itype = -1;
 132     }
 133
 134     if (itype == -1) {
 135         *info = -1;
 136     } else if (*cfrom == 0.f) {
 137         *info = -4;
 138     } else if (*m < 0) {
 139         *info = -6;
 140     } else if (*n < 0 || itype == 4 && *n != *m || itype == 5 && *n != *m) {
 141         *info = -7;
 142     } else if (itype <= 3 && *lda < max(1,*m)) {
 143         *info = -9;
 144     } else if (itype >= 4) {
 145 /* Computing MAX */
 146         i__1 = *m - 1;
 147         if (*kl < 0 || *kl > max(i__1,0)) {
 148             *info = -2;
 149         } else /* if(complicated condition) */ {
 150 /* Computing MAX */
 151             i__1 = *n - 1;
 152             if (*ku < 0 || *ku > max(i__1,0) || (itype == 4 || itype == 5) &&
 153                     *kl != *ku) {
 154                 *info = -3;
 155             } else if (itype == 4 && *lda < *kl + 1 || itype == 5 && *lda < *
 156                     ku + 1 || itype == 6 && *lda < (*kl << 1) + *ku + 1) {
 157                 *info = -9;
 158             }
 159         }
 160     }
 161
 162     if (*info != 0) {
 163         i__1 = -(*info);
 164         xerbla_("SLASCL", &i__1);
 165         return 0;
 166     }
 167
 168 /*     Quick return if possible */
 169
 170     if (*n == 0 || *m == 0) {
 171         return 0;
 172     }
 173
 174 /*     Get machine parameters */
 175
 176     smlnum = slamch_("S");
 177     bignum = 1.f / smlnum;
 178
 179     cfromc = *cfrom;
 180     ctoc = *cto;
 181
 182 L10:
 183     cfrom1 = cfromc * smlnum;
 184     cto1 = ctoc / bignum;
 185     if (dabs(cfrom1) > dabs(ctoc) && ctoc != 0.f) {
 186         mul = smlnum;
 187         done = FALSE_;
 188         cfromc = cfrom1;
 189     } else if (dabs(cto1) > dabs(cfromc)) {
 190         mul = bignum;
 191         done = FALSE_;
 192         ctoc = cto1;
 193     } else {
 194         mul = ctoc / cfromc;
 195         done = TRUE_;
 196     }
 197
 198     if (itype == 0) {
 199
 200 /*        Full matrix */
 201
 202         i__1 = *n;
 203         for (j = 1; j <= i__1; ++j) {
 204             i__2 = *m;
 205             for (i__ = 1; i__ <= i__2; ++i__) {
 206                 a[i__ + j * a_dim1] *= mul;
 207 /* L20: */
 208             }
 209 /* L30: */
 210         }
 211
 212     } else if (itype == 1) {
 213
 214 /*        Lower triangular matrix */
 215
 216         i__1 = *n;
 217         for (j = 1; j <= i__1; ++j) {
 218             i__2 = *m;
 219             for (i__ = j; i__ <= i__2; ++i__) {
 220                 a[i__ + j * a_dim1] *= mul;
 221 /* L40: */
 222             }
 223 /* L50: */
 224         }
 225
 226     } else if (itype == 2) {
 227
 228 /*        Upper triangular matrix */
 229
 230         i__1 = *n;
 231         for (j = 1; j <= i__1; ++j) {
 232             i__2 = min(j,*m);
 233             for (i__ = 1; i__ <= i__2; ++i__) {
 234                 a[i__ + j * a_dim1] *= mul;
 235 /* L60: */
 236             }
 237 /* L70: */
 238         }
 239
 240     } else if (itype == 3) {
 241
 242 /*        Upper Hessenberg matrix */
 243
 244         i__1 = *n;
 245         for (j = 1; j <= i__1; ++j) {
 246 /* Computing MIN */
 247             i__3 = j + 1;
 248             i__2 = min(i__3,*m);
 249             for (i__ = 1; i__ <= i__2; ++i__) {
 250                 a[i__ + j * a_dim1] *= mul;
 251 /* L80: */
 252             }
 253 /* L90: */
 254         }
 255
 256     } else if (itype == 4) {
 257
 258 /*        Lower half of a symmetric band matrix */
 259
 260         k3 = *kl + 1;
 261         k4 = *n + 1;
 262         i__1 = *n;
 263         for (j = 1; j <= i__1; ++j) {
 264 /* Computing MIN */
 265             i__3 = k3, i__4 = k4 - j;
 266             i__2 = min(i__3,i__4);
 267             for (i__ = 1; i__ <= i__2; ++i__) {
 268                 a[i__ + j * a_dim1] *= mul;
 269 /* L100: */
 270             }
 271 /* L110: */
 272         }
 273
 274     } else if (itype == 5) {
 275
 276 /*        Upper half of a symmetric band matrix */
 277
 278         k1 = *ku + 2;
 279         k3 = *ku + 1;
 280         i__1 = *n;
 281         for (j = 1; j <= i__1; ++j) {
 282 /* Computing MAX */
 283             i__2 = k1 - j;
 284             i__3 = k3;
 285             for (i__ = max(i__2,1); i__ <= i__3; ++i__) {
 286                 a[i__ + j * a_dim1] *= mul;
 287 /* L120: */
 288             }
 289 /* L130: */
 290         }
 291
 292     } else if (itype == 6) {
 293
 294 /*        Band matrix */
 295
 296         k1 = *kl + *ku + 2;
 297         k2 = *kl + 1;
 298         k3 = (*kl << 1) + *ku + 1;
 299         k4 = *kl + *ku + 1 + *m;
 300         i__1 = *n;
 301         for (j = 1; j <= i__1; ++j) {
 302 /* Computing MAX */
 303             i__3 = k1 - j;
 304 /* Computing MIN */
 305             i__4 = k3, i__5 = k4 - j;
 306             i__2 = min(i__4,i__5);
 307             for (i__ = max(i__3,k2); i__ <= i__2; ++i__) {
 308                 a[i__ + j * a_dim1] *= mul;
 309 /* L140: */
 310             }
 311 /* L150: */
 312         }
 313
 314     }
 315
 316     if (! done) {
 317         goto L10;
 318     }
 319
 320     return 0;
 321
 322 /*     End of SLASCL */
 323
 324 } /* slascl_ */