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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c__2 = 2;
   7 static integer c__0 = 0;
   8
   9 /* Subroutine */ int slasq1_(integer *n, real *d__, real *e, real *work,
  10         integer *info)
  11 {
  12     /* System generated locals */
  13     integer i__1, i__2;
  14     real r__1, r__2, r__3;
  15
  16     /* Builtin functions */
  17     double sqrt(doublereal);
  18
  19     /* Local variables */
  20     integer i__;
  21     real eps;
  22     extern /* Subroutine */ int slas2_(real *, real *, real *, real *, real *)
  23             ;
  24     real scale;
  25     integer iinfo;
  26     real sigmn, sigmx;
  27     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  28             integer *), slasq2_(integer *, real *, integer *);
  29     extern doublereal slamch_(char *);
  30     real safmin;
  31     extern /* Subroutine */ int xerbla_(char *, integer *), slascl_(
  32             char *, integer *, integer *, real *, real *, integer *, integer *
  33 , real *, integer *, integer *), slasrt_(char *, integer *
  34 , real *, integer *);
  35
  36
  37 /*  -- LAPACK routine (version 3.1) -- */
  38 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  39 /*     November 2006 */
  40
  41 /*     .. Scalar Arguments .. */
  42 /*     .. */
  43 /*     .. Array Arguments .. */
  44 /*     .. */
  45
  46 /*  Purpose */
  47 /*  ======= */
  48
  49 /*  SLASQ1 computes the singular values of a real N-by-N bidiagonal */
  50 /*  matrix with diagonal D and off-diagonal E. The singular values */
  51 /*  are computed to high relative accuracy, in the absence of */
  52 /*  denormalization, underflow and overflow. The algorithm was first */
  53 /*  presented in */
  54
  55 /*  "Accurate singular values and differential qd algorithms" by K. V. */
  56 /*  Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, */
  57 /*  1994, */
  58
  59 /*  and the present implementation is described in "An implementation of */
  60 /*  the dqds Algorithm (Positive Case)", LAPACK Working Note. */
  61
  62 /*  Arguments */
  63 /*  ========= */
  64
  65 /*  N     (input) INTEGER */
  66 /*        The number of rows and columns in the matrix. N >= 0. */
  67
  68 /*  D     (input/output) REAL array, dimension (N) */
  69 /*        On entry, D contains the diagonal elements of the */
  70 /*        bidiagonal matrix whose SVD is desired. On normal exit, */
  71 /*        D contains the singular values in decreasing order. */
  72
  73 /*  E     (input/output) REAL array, dimension (N) */
  74 /*        On entry, elements E(1:N-1) contain the off-diagonal elements */
  75 /*        of the bidiagonal matrix whose SVD is desired. */
  76 /*        On exit, E is overwritten. */
  77
  78 /*  WORK  (workspace) REAL array, dimension (4*N) */
  79
  80 /*  INFO  (output) INTEGER */
  81 /*        = 0: successful exit */
  82 /*        < 0: if INFO = -i, the i-th argument had an illegal value */
  83 /*        > 0: the algorithm failed */
  84 /*             = 1, a split was marked by a positive value in E */
  85 /*             = 2, current block of Z not diagonalized after 30*N */
  86 /*                  iterations (in inner while loop) */
  87 /*             = 3, termination criterion of outer while loop not met */
  88 /*                  (program created more than N unreduced blocks) */
  89
  90 /*  ===================================================================== */
  91
  92 /*     .. Parameters .. */
  93 /*     .. */
  94 /*     .. Local Scalars .. */
  95 /*     .. */
  96 /*     .. External Subroutines .. */
  97 /*     .. */
  98 /*     .. External Functions .. */
  99 /*     .. */
 100 /*     .. Intrinsic Functions .. */
 101 /*     .. */
 102 /*     .. Executable Statements .. */
 103
 104     /* Parameter adjustments */
 105     --work;
 106     --e;
 107     --d__;
 108
 109     /* Function Body */
 110     *info = 0;
 111     if (*n < 0) {
 112         *info = -2;
 113         i__1 = -(*info);
 114         xerbla_("SLASQ1", &i__1);
 115         return 0;
 116     } else if (*n == 0) {
 117         return 0;
 118     } else if (*n == 1) {
 119         d__[1] = dabs(d__[1]);
 120         return 0;
 121     } else if (*n == 2) {
 122         slas2_(&d__[1], &e[1], &d__[2], &sigmn, &sigmx);
 123         d__[1] = sigmx;
 124         d__[2] = sigmn;
 125         return 0;
 126     }
 127
 128 /*     Estimate the largest singular value. */
 129
 130     sigmx = 0.f;
 131     i__1 = *n - 1;
 132     for (i__ = 1; i__ <= i__1; ++i__) {
 133         d__[i__] = (r__1 = d__[i__], dabs(r__1));
 134 /* Computing MAX */
 135         r__2 = sigmx, r__3 = (r__1 = e[i__], dabs(r__1));
 136         sigmx = dmax(r__2,r__3);
 137 /* L10: */
 138     }
 139     d__[*n] = (r__1 = d__[*n], dabs(r__1));
 140
 141 /*     Early return if SIGMX is zero (matrix is already diagonal). */
 142
 143     if (sigmx == 0.f) {
 144         slasrt_("D", n, &d__[1], &iinfo);
 145         return 0;
 146     }
 147
 148     i__1 = *n;
 149     for (i__ = 1; i__ <= i__1; ++i__) {
 150 /* Computing MAX */
 151         r__1 = sigmx, r__2 = d__[i__];
 152         sigmx = dmax(r__1,r__2);
 153 /* L20: */
 154     }
 155
 156 /*     Copy D and E into WORK (in the Z format) and scale (squaring the */
 157 /*     input data makes scaling by a power of the radix pointless). */
 158
 159     eps = slamch_("Precision");
 160     safmin = slamch_("Safe minimum");
 161     scale = sqrt(eps / safmin);
 162     scopy_(n, &d__[1], &c__1, &work[1], &c__2);
 163     i__1 = *n - 1;
 164     scopy_(&i__1, &e[1], &c__1, &work[2], &c__2);
 165     i__1 = (*n << 1) - 1;
 166     i__2 = (*n << 1) - 1;
 167     slascl_("G", &c__0, &c__0, &sigmx, &scale, &i__1, &c__1, &work[1], &i__2,
 168             &iinfo);
 169
 170 /*     Compute the q's and e's. */
 171
 172     i__1 = (*n << 1) - 1;
 173     for (i__ = 1; i__ <= i__1; ++i__) {
 174 /* Computing 2nd power */
 175         r__1 = work[i__];
 176         work[i__] = r__1 * r__1;
 177 /* L30: */
 178     }
 179     work[*n * 2] = 0.f;
 180
 181     slasq2_(n, &work[1], info);
 182
 183     if (*info == 0) {
 184         i__1 = *n;
 185         for (i__ = 1; i__ <= i__1; ++i__) {
 186             d__[i__] = sqrt(work[i__]);
 187 /* L40: */
 188         }
 189         slascl_("G", &c__0, &c__0, &scale, &sigmx, n, &c__1, &d__[1], n, &
 190                 iinfo);
 191     }
 192
 193     return 0;
 194
 195 /*     End of SLASQ1 */
 196
 197 } /* slasq1_ */