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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slagts_(integer *job, integer *n, real *a, real *b, real
   4         *c__, real *d__, integer *in, real *y, real *tol, integer *info)
   5 {
   6     /* System generated locals */
   7     integer i__1;
   8     real r__1, r__2, r__3, r__4, r__5;
   9
  10     /* Builtin functions */
  11     double r_sign(real *, real *);
  12
  13     /* Local variables */
  14     integer k;
  15     real ak, eps, temp, pert, absak, sfmin;
  16     extern doublereal slamch_(char *);
  17     extern /* Subroutine */ int xerbla_(char *, integer *);
  18     real bignum;
  19
  20
  21 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  22 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  23 /*     November 2006 */
  24
  25 /*     .. Scalar Arguments .. */
  26 /*     .. */
  27 /*     .. Array Arguments .. */
  28 /*     .. */
  29
  30 /*  Purpose */
  31 /*  ======= */
  32
  33 /*  SLAGTS may be used to solve one of the systems of equations */
  34
  35 /*     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y, */
  36
  37 /*  where T is an n by n tridiagonal matrix, for x, following the */
  38 /*  factorization of (T - lambda*I) as */
  39
  40 /*     (T - lambda*I) = P*L*U , */
  41
  42 /*  by routine SLAGTF. The choice of equation to be solved is */
  43 /*  controlled by the argument JOB, and in each case there is an option */
  44 /*  to perturb zero or very small diagonal elements of U, this option */
  45 /*  being intended for use in applications such as inverse iteration. */
  46
  47 /*  Arguments */
  48 /*  ========= */
  49
  50 /*  JOB     (input) INTEGER */
  51 /*          Specifies the job to be performed by SLAGTS as follows: */
  52 /*          =  1: The equations  (T - lambda*I)x = y  are to be solved, */
  53 /*                but diagonal elements of U are not to be perturbed. */
  54 /*          = -1: The equations  (T - lambda*I)x = y  are to be solved */
  55 /*                and, if overflow would otherwise occur, the diagonal */
  56 /*                elements of U are to be perturbed. See argument TOL */
  57 /*                below. */
  58 /*          =  2: The equations  (T - lambda*I)'x = y  are to be solved, */
  59 /*                but diagonal elements of U are not to be perturbed. */
  60 /*          = -2: The equations  (T - lambda*I)'x = y  are to be solved */
  61 /*                and, if overflow would otherwise occur, the diagonal */
  62 /*                elements of U are to be perturbed. See argument TOL */
  63 /*                below. */
  64
  65 /*  N       (input) INTEGER */
  66 /*          The order of the matrix T. */
  67
  68 /*  A       (input) REAL array, dimension (N) */
  69 /*          On entry, A must contain the diagonal elements of U as */
  70 /*          returned from SLAGTF. */
  71
  72 /*  B       (input) REAL array, dimension (N-1) */
  73 /*          On entry, B must contain the first super-diagonal elements of */
  74 /*          U as returned from SLAGTF. */
  75
  76 /*  C       (input) REAL array, dimension (N-1) */
  77 /*          On entry, C must contain the sub-diagonal elements of L as */
  78 /*          returned from SLAGTF. */
  79
  80 /*  D       (input) REAL array, dimension (N-2) */
  81 /*          On entry, D must contain the second super-diagonal elements */
  82 /*          of U as returned from SLAGTF. */
  83
  84 /*  IN      (input) INTEGER array, dimension (N) */
  85 /*          On entry, IN must contain details of the matrix P as returned */
  86 /*          from SLAGTF. */
  87
  88 /*  Y       (input/output) REAL array, dimension (N) */
  89 /*          On entry, the right hand side vector y. */
  90 /*          On exit, Y is overwritten by the solution vector x. */
  91
  92 /*  TOL     (input/output) REAL */
  93 /*          On entry, with  JOB .lt. 0, TOL should be the minimum */
  94 /*          perturbation to be made to very small diagonal elements of U. */
  95 /*          TOL should normally be chosen as about eps*norm(U), where eps */
  96 /*          is the relative machine precision, but if TOL is supplied as */
  97 /*          non-positive, then it is reset to eps*max( abs( u(i,j) ) ). */
  98 /*          If  JOB .gt. 0  then TOL is not referenced. */
  99
 100 /*          On exit, TOL is changed as described above, only if TOL is */
 101 /*          non-positive on entry. Otherwise TOL is unchanged. */
 102
 103 /*  INFO    (output) INTEGER */
 104 /*          = 0   : successful exit */
 105 /*          .lt. 0: if INFO = -i, the i-th argument had an illegal value */
 106 /*          .gt. 0: overflow would occur when computing the INFO(th) */
 107 /*                  element of the solution vector x. This can only occur */
 108 /*                  when JOB is supplied as positive and either means */
 109 /*                  that a diagonal element of U is very small, or that */
 110 /*                  the elements of the right-hand side vector y are very */
 111 /*                  large. */
 112
 113 /*  ===================================================================== */
 114
 115 /*     .. Parameters .. */
 116 /*     .. */
 117 /*     .. Local Scalars .. */
 118 /*     .. */
 119 /*     .. Intrinsic Functions .. */
 120 /*     .. */
 121 /*     .. External Functions .. */
 122 /*     .. */
 123 /*     .. External Subroutines .. */
 124 /*     .. */
 125 /*     .. Executable Statements .. */
 126
 127     /* Parameter adjustments */
 128     --y;
 129     --in;
 130     --d__;
 131     --c__;
 132     --b;
 133     --a;
 134
 135     /* Function Body */
 136     *info = 0;
 137     if (abs(*job) > 2 || *job == 0) {
 138         *info = -1;
 139     } else if (*n < 0) {
 140         *info = -2;
 141     }
 142     if (*info != 0) {
 143         i__1 = -(*info);
 144         xerbla_("SLAGTS", &i__1);
 145         return 0;
 146     }
 147
 148     if (*n == 0) {
 149         return 0;
 150     }
 151
 152     eps = slamch_("Epsilon");
 153     sfmin = slamch_("Safe minimum");
 154     bignum = 1.f / sfmin;
 155
 156     if (*job < 0) {
 157         if (*tol <= 0.f) {
 158             *tol = dabs(a[1]);
 159             if (*n > 1) {
 160 /* Computing MAX */
 161                 r__1 = *tol, r__2 = dabs(a[2]), r__1 = max(r__1,r__2), r__2 =
 162                         dabs(b[1]);
 163                 *tol = dmax(r__1,r__2);
 164             }
 165             i__1 = *n;
 166             for (k = 3; k <= i__1; ++k) {
 167 /* Computing MAX */
 168                 r__4 = *tol, r__5 = (r__1 = a[k], dabs(r__1)), r__4 = max(
 169                         r__4,r__5), r__5 = (r__2 = b[k - 1], dabs(r__2)),
 170                         r__4 = max(r__4,r__5), r__5 = (r__3 = d__[k - 2],
 171                         dabs(r__3));
 172                 *tol = dmax(r__4,r__5);
 173 /* L10: */
 174             }
 175             *tol *= eps;
 176             if (*tol == 0.f) {
 177                 *tol = eps;
 178             }
 179         }
 180     }
 181
 182     if (abs(*job) == 1) {
 183         i__1 = *n;
 184         for (k = 2; k <= i__1; ++k) {
 185             if (in[k - 1] == 0) {
 186                 y[k] -= c__[k - 1] * y[k - 1];
 187             } else {
 188                 temp = y[k - 1];
 189                 y[k - 1] = y[k];
 190                 y[k] = temp - c__[k - 1] * y[k];
 191             }
 192 /* L20: */
 193         }
 194         if (*job == 1) {
 195             for (k = *n; k >= 1; --k) {
 196                 if (k <= *n - 2) {
 197                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
 198                 } else if (k == *n - 1) {
 199                     temp = y[k] - b[k] * y[k + 1];
 200                 } else {
 201                     temp = y[k];
 202                 }
 203                 ak = a[k];
 204                 absak = dabs(ak);
 205                 if (absak < 1.f) {
 206                     if (absak < sfmin) {
 207                         if (absak == 0.f || dabs(temp) * sfmin > absak) {
 208                             *info = k;
 209                             return 0;
 210                         } else {
 211                             temp *= bignum;
 212                             ak *= bignum;
 213                         }
 214                     } else if (dabs(temp) > absak * bignum) {
 215                         *info = k;
 216                         return 0;
 217                     }
 218                 }
 219                 y[k] = temp / ak;
 220 /* L30: */
 221             }
 222         } else {
 223             for (k = *n; k >= 1; --k) {
 224                 if (k <= *n - 2) {
 225                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
 226                 } else if (k == *n - 1) {
 227                     temp = y[k] - b[k] * y[k + 1];
 228                 } else {
 229                     temp = y[k];
 230                 }
 231                 ak = a[k];
 232                 pert = r_sign(tol, &ak);
 233 L40:
 234                 absak = dabs(ak);
 235                 if (absak < 1.f) {
 236                     if (absak < sfmin) {
 237                         if (absak == 0.f || dabs(temp) * sfmin > absak) {
 238                             ak += pert;
 239                             pert *= 2;
 240                             goto L40;
 241                         } else {
 242                             temp *= bignum;
 243                             ak *= bignum;
 244                         }
 245                     } else if (dabs(temp) > absak * bignum) {
 246                         ak += pert;
 247                         pert *= 2;
 248                         goto L40;
 249                     }
 250                 }
 251                 y[k] = temp / ak;
 252 /* L50: */
 253             }
 254         }
 255     } else {
 256
 257 /*        Come to here if  JOB = 2 or -2 */
 258
 259         if (*job == 2) {
 260             i__1 = *n;
 261             for (k = 1; k <= i__1; ++k) {
 262                 if (k >= 3) {
 263                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 264                 } else if (k == 2) {
 265                     temp = y[k] - b[k - 1] * y[k - 1];
 266                 } else {
 267                     temp = y[k];
 268                 }
 269                 ak = a[k];
 270                 absak = dabs(ak);
 271                 if (absak < 1.f) {
 272                     if (absak < sfmin) {
 273                         if (absak == 0.f || dabs(temp) * sfmin > absak) {
 274                             *info = k;
 275                             return 0;
 276                         } else {
 277                             temp *= bignum;
 278                             ak *= bignum;
 279                         }
 280                     } else if (dabs(temp) > absak * bignum) {
 281                         *info = k;
 282                         return 0;
 283                     }
 284                 }
 285                 y[k] = temp / ak;
 286 /* L60: */
 287             }
 288         } else {
 289             i__1 = *n;
 290             for (k = 1; k <= i__1; ++k) {
 291                 if (k >= 3) {
 292                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 293                 } else if (k == 2) {
 294                     temp = y[k] - b[k - 1] * y[k - 1];
 295                 } else {
 296                     temp = y[k];
 297                 }
 298                 ak = a[k];
 299                 pert = r_sign(tol, &ak);
 300 L70:
 301                 absak = dabs(ak);
 302                 if (absak < 1.f) {
 303                     if (absak < sfmin) {
 304                         if (absak == 0.f || dabs(temp) * sfmin > absak) {
 305                             ak += pert;
 306                             pert *= 2;
 307                             goto L70;
 308                         } else {
 309                             temp *= bignum;
 310                             ak *= bignum;
 311                         }
 312                     } else if (dabs(temp) > absak * bignum) {
 313                         ak += pert;
 314                         pert *= 2;
 315                         goto L70;
 316                     }
 317                 }
 318                 y[k] = temp / ak;
 319 /* L80: */
 320             }
 321         }
 322
 323         for (k = *n; k >= 2; --k) {
 324             if (in[k - 1] == 0) {
 325                 y[k - 1] -= c__[k - 1] * y[k];
 326             } else {
 327                 temp = y[k - 1];
 328                 y[k - 1] = y[k];
 329                 y[k] = temp - c__[k - 1] * y[k];
 330             }
 331 /* L90: */
 332         }
 333     }
 334
 335 /*     End of SLAGTS */
 336
 337     return 0;
 338 } /* slagts_ */