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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 /* Subroutine */ int dlasdq_(char *uplo, integer *sqre, integer *n, integer *
   8         ncvt, integer *nru, integer *ncc, doublereal *d__, doublereal *e,
   9         doublereal *vt, integer *ldvt, doublereal *u, integer *ldu,
  10         doublereal *c__, integer *ldc, doublereal *work, integer *info)
  11 {
  12     /* System generated locals */
  13     integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  14             i__2;
  15
  16     /* Local variables */
  17     integer i__, j;
  18     doublereal r__, cs, sn;
  19     integer np1, isub;
  20     doublereal smin;
  21     integer sqre1;
  22     extern logical lsame_(char *, char *);
  23     extern /* Subroutine */ int dlasr_(char *, char *, char *, integer *,
  24             integer *, doublereal *, doublereal *, doublereal *, integer *), dswap_(integer *, doublereal *, integer *
  25 , doublereal *, integer *);
  26     integer iuplo;
  27     extern /* Subroutine */ int dlartg_(doublereal *, doublereal *,
  28             doublereal *, doublereal *, doublereal *), xerbla_(char *,
  29             integer *), dbdsqr_(char *, integer *, integer *, integer
  30             *, integer *, doublereal *, doublereal *, doublereal *, integer *,
  31              doublereal *, integer *, doublereal *, integer *, doublereal *,
  32             integer *);
  33     logical rotate;
  34
  35
  36 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  37 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  38 /*     November 2006 */
  39
  40 /*     .. Scalar Arguments .. */
  41 /*     .. */
  42 /*     .. Array Arguments .. */
  43 /*     .. */
  44
  45 /*  Purpose */
  46 /*  ======= */
  47
  48 /*  DLASDQ computes the singular value decomposition (SVD) of a real */
  49 /*  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal */
  50 /*  E, accumulating the transformations if desired. Letting B denote */
  51 /*  the input bidiagonal matrix, the algorithm computes orthogonal */
  52 /*  matrices Q and P such that B = Q * S * P' (P' denotes the transpose */
  53 /*  of P). The singular values S are overwritten on D. */
  54
  55 /*  The input matrix U  is changed to U  * Q  if desired. */
  56 /*  The input matrix VT is changed to P' * VT if desired. */
  57 /*  The input matrix C  is changed to Q' * C  if desired. */
  58
  59 /*  See "Computing  Small Singular Values of Bidiagonal Matrices With */
  60 /*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
  61 /*  LAPACK Working Note #3, for a detailed description of the algorithm. */
  62
  63 /*  Arguments */
  64 /*  ========= */
  65
  66 /*  UPLO  (input) CHARACTER*1 */
  67 /*        On entry, UPLO specifies whether the input bidiagonal matrix */
  68 /*        is upper or lower bidiagonal, and wether it is square are */
  69 /*        not. */
  70 /*           UPLO = 'U' or 'u'   B is upper bidiagonal. */
  71 /*           UPLO = 'L' or 'l'   B is lower bidiagonal. */
  72
  73 /*  SQRE  (input) INTEGER */
  74 /*        = 0: then the input matrix is N-by-N. */
  75 /*        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and */
  76 /*             (N+1)-by-N if UPLU = 'L'. */
  77
  78 /*        The bidiagonal matrix has */
  79 /*        N = NL + NR + 1 rows and */
  80 /*        M = N + SQRE >= N columns. */
  81
  82 /*  N     (input) INTEGER */
  83 /*        On entry, N specifies the number of rows and columns */
  84 /*        in the matrix. N must be at least 0. */
  85
  86 /*  NCVT  (input) INTEGER */
  87 /*        On entry, NCVT specifies the number of columns of */
  88 /*        the matrix VT. NCVT must be at least 0. */
  89
  90 /*  NRU   (input) INTEGER */
  91 /*        On entry, NRU specifies the number of rows of */
  92 /*        the matrix U. NRU must be at least 0. */
  93
  94 /*  NCC   (input) INTEGER */
  95 /*        On entry, NCC specifies the number of columns of */
  96 /*        the matrix C. NCC must be at least 0. */
  97
  98 /*  D     (input/output) DOUBLE PRECISION array, dimension (N) */
  99 /*        On entry, D contains the diagonal entries of the */
 100 /*        bidiagonal matrix whose SVD is desired. On normal exit, */
 101 /*        D contains the singular values in ascending order. */
 102
 103 /*  E     (input/output) DOUBLE PRECISION array. */
 104 /*        dimension is (N-1) if SQRE = 0 and N if SQRE = 1. */
 105 /*        On entry, the entries of E contain the offdiagonal entries */
 106 /*        of the bidiagonal matrix whose SVD is desired. On normal */
 107 /*        exit, E will contain 0. If the algorithm does not converge, */
 108 /*        D and E will contain the diagonal and superdiagonal entries */
 109 /*        of a bidiagonal matrix orthogonally equivalent to the one */
 110 /*        given as input. */
 111
 112 /*  VT    (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) */
 113 /*        On entry, contains a matrix which on exit has been */
 114 /*        premultiplied by P', dimension N-by-NCVT if SQRE = 0 */
 115 /*        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). */
 116
 117 /*  LDVT  (input) INTEGER */
 118 /*        On entry, LDVT specifies the leading dimension of VT as */
 119 /*        declared in the calling (sub) program. LDVT must be at */
 120 /*        least 1. If NCVT is nonzero LDVT must also be at least N. */
 121
 122 /*  U     (input/output) DOUBLE PRECISION array, dimension (LDU, N) */
 123 /*        On entry, contains a  matrix which on exit has been */
 124 /*        postmultiplied by Q, dimension NRU-by-N if SQRE = 0 */
 125 /*        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). */
 126
 127 /*  LDU   (input) INTEGER */
 128 /*        On entry, LDU  specifies the leading dimension of U as */
 129 /*        declared in the calling (sub) program. LDU must be at */
 130 /*        least max( 1, NRU ) . */
 131
 132 /*  C     (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) */
 133 /*        On entry, contains an N-by-NCC matrix which on exit */
 134 /*        has been premultiplied by Q'  dimension N-by-NCC if SQRE = 0 */
 135 /*        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). */
 136
 137 /*  LDC   (input) INTEGER */
 138 /*        On entry, LDC  specifies the leading dimension of C as */
 139 /*        declared in the calling (sub) program. LDC must be at */
 140 /*        least 1. If NCC is nonzero, LDC must also be at least N. */
 141
 142 /*  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N) */
 143 /*        Workspace. Only referenced if one of NCVT, NRU, or NCC is */
 144 /*        nonzero, and if N is at least 2. */
 145
 146 /*  INFO  (output) INTEGER */
 147 /*        On exit, a value of 0 indicates a successful exit. */
 148 /*        If INFO < 0, argument number -INFO is illegal. */
 149 /*        If INFO > 0, the algorithm did not converge, and INFO */
 150 /*        specifies how many superdiagonals did not converge. */
 151
 152 /*  Further Details */
 153 /*  =============== */
 154
 155 /*  Based on contributions by */
 156 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 157 /*     California at Berkeley, USA */
 158
 159 /*  ===================================================================== */
 160
 161 /*     .. Parameters .. */
 162 /*     .. */
 163 /*     .. Local Scalars .. */
 164 /*     .. */
 165 /*     .. External Subroutines .. */
 166 /*     .. */
 167 /*     .. External Functions .. */
 168 /*     .. */
 169 /*     .. Intrinsic Functions .. */
 170 /*     .. */
 171 /*     .. Executable Statements .. */
 172
 173 /*     Test the input parameters. */
 174
 175     /* Parameter adjustments */
 176     --d__;
 177     --e;
 178     vt_dim1 = *ldvt;
 179     vt_offset = 1 + vt_dim1;
 180     vt -= vt_offset;
 181     u_dim1 = *ldu;
 182     u_offset = 1 + u_dim1;
 183     u -= u_offset;
 184     c_dim1 = *ldc;
 185     c_offset = 1 + c_dim1;
 186     c__ -= c_offset;
 187     --work;
 188
 189     /* Function Body */
 190     *info = 0;
 191     iuplo = 0;
 192     if (lsame_(uplo, "U")) {
 193         iuplo = 1;
 194     }
 195     if (lsame_(uplo, "L")) {
 196         iuplo = 2;
 197     }
 198     if (iuplo == 0) {
 199         *info = -1;
 200     } else if (*sqre < 0 || *sqre > 1) {
 201         *info = -2;
 202     } else if (*n < 0) {
 203         *info = -3;
 204     } else if (*ncvt < 0) {
 205         *info = -4;
 206     } else if (*nru < 0) {
 207         *info = -5;
 208     } else if (*ncc < 0) {
 209         *info = -6;
 210     } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) {
 211         *info = -10;
 212     } else if (*ldu < max(1,*nru)) {
 213         *info = -12;
 214     } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) {
 215         *info = -14;
 216     }
 217     if (*info != 0) {
 218         i__1 = -(*info);
 219         xerbla_("DLASDQ", &i__1);
 220         return 0;
 221     }
 222     if (*n == 0) {
 223         return 0;
 224     }
 225
 226 /*     ROTATE is true if any singular vectors desired, false otherwise */
 227
 228     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
 229     np1 = *n + 1;
 230     sqre1 = *sqre;
 231
 232 /*     If matrix non-square upper bidiagonal, rotate to be lower */
 233 /*     bidiagonal.  The rotations are on the right. */
 234
 235     if (iuplo == 1 && sqre1 == 1) {
 236         i__1 = *n - 1;
 237         for (i__ = 1; i__ <= i__1; ++i__) {
 238             dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 239             d__[i__] = r__;
 240             e[i__] = sn * d__[i__ + 1];
 241             d__[i__ + 1] = cs * d__[i__ + 1];
 242             if (rotate) {
 243                 work[i__] = cs;
 244                 work[*n + i__] = sn;
 245             }
 246 /* L10: */
 247         }
 248         dlartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
 249         d__[*n] = r__;
 250         e[*n] = 0.;
 251         if (rotate) {
 252             work[*n] = cs;
 253             work[*n + *n] = sn;
 254         }
 255         iuplo = 2;
 256         sqre1 = 0;
 257
 258 /*        Update singular vectors if desired. */
 259
 260         if (*ncvt > 0) {
 261             dlasr_("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[
 262                     vt_offset], ldvt);
 263         }
 264     }
 265
 266 /*     If matrix lower bidiagonal, rotate to be upper bidiagonal */
 267 /*     by applying Givens rotations on the left. */
 268
 269     if (iuplo == 2) {
 270         i__1 = *n - 1;
 271         for (i__ = 1; i__ <= i__1; ++i__) {
 272             dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 273             d__[i__] = r__;
 274             e[i__] = sn * d__[i__ + 1];
 275             d__[i__ + 1] = cs * d__[i__ + 1];
 276             if (rotate) {
 277                 work[i__] = cs;
 278                 work[*n + i__] = sn;
 279             }
 280 /* L20: */
 281         }
 282
 283 /*        If matrix (N+1)-by-N lower bidiagonal, one additional */
 284 /*        rotation is needed. */
 285
 286         if (sqre1 == 1) {
 287             dlartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
 288             d__[*n] = r__;
 289             if (rotate) {
 290                 work[*n] = cs;
 291                 work[*n + *n] = sn;
 292             }
 293         }
 294
 295 /*        Update singular vectors if desired. */
 296
 297         if (*nru > 0) {
 298             if (sqre1 == 0) {
 299                 dlasr_("R", "V", "F", nru, n, &work[1], &work[np1], &u[
 300                         u_offset], ldu);
 301             } else {
 302                 dlasr_("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[
 303                         u_offset], ldu);
 304             }
 305         }
 306         if (*ncc > 0) {
 307             if (sqre1 == 0) {
 308                 dlasr_("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[
 309                         c_offset], ldc);
 310             } else {
 311                 dlasr_("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[
 312                         c_offset], ldc);
 313             }
 314         }
 315     }
 316
 317 /*     Call DBDSQR to compute the SVD of the reduced real */
 318 /*     N-by-N upper bidiagonal matrix. */
 319
 320     dbdsqr_("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[
 321             u_offset], ldu, &c__[c_offset], ldc, &work[1], info);
 322
 323 /*     Sort the singular values into ascending order (insertion sort on */
 324 /*     singular values, but only one transposition per singular vector) */
 325
 326     i__1 = *n;
 327     for (i__ = 1; i__ <= i__1; ++i__) {
 328
 329 /*        Scan for smallest D(I). */
 330
 331         isub = i__;
 332         smin = d__[i__];
 333         i__2 = *n;
 334         for (j = i__ + 1; j <= i__2; ++j) {
 335             if (d__[j] < smin) {
 336                 isub = j;
 337                 smin = d__[j];
 338             }
 339 /* L30: */
 340         }
 341         if (isub != i__) {
 342
 343 /*           Swap singular values and vectors. */
 344
 345             d__[isub] = d__[i__];
 346             d__[i__] = smin;
 347             if (*ncvt > 0) {
 348                 dswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1],
 349                         ldvt);
 350             }
 351             if (*nru > 0) {
 352                 dswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1]
 353 , &c__1);
 354             }
 355             if (*ncc > 0) {
 356                 dswap_(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc)
 357                         ;
 358             }
 359         }
 360 /* L40: */
 361     }
 362
 363     return 0;
 364
 365 /*     End of DLASDQ */
 366
 367 } /* dlasdq_ */