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