vcs.maemo.org Git - opencv/blob - 3rdparty/lapack/slasd1.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__0 = 0;
   6 static real c_b7 = 1.f;
   7 static integer c__1 = 1;
   8 static integer c_n1 = -1;
   9
  10 /* Subroutine */ int slasd1_(integer *nl, integer *nr, integer *sqre, real *
  11         d__, real *alpha, real *beta, real *u, integer *ldu, real *vt,
  12         integer *ldvt, integer *idxq, integer *iwork, real *work, integer *
  13         info)
  14 {
  15     /* System generated locals */
  16     integer u_dim1, u_offset, vt_dim1, vt_offset, i__1;
  17     real r__1, r__2;
  18
  19     /* Local variables */
  20     integer i__, k, m, n, n1, n2, iq, iz, iu2, ldq, idx, ldu2, ivt2, idxc,
  21             idxp, ldvt2;
  22     extern /* Subroutine */ int slasd2_(integer *, integer *, integer *,
  23             integer *, real *, real *, real *, real *, real *, integer *,
  24             real *, integer *, real *, real *, integer *, real *, integer *,
  25             integer *, integer *, integer *, integer *, integer *, integer *),
  26              slasd3_(integer *, integer *, integer *, integer *, real *, real
  27             *, integer *, real *, real *, integer *, real *, integer *, real *
  28 , integer *, real *, integer *, integer *, integer *, real *,
  29             integer *);
  30     integer isigma;
  31     extern /* Subroutine */ int xerbla_(char *, integer *), slascl_(
  32             char *, integer *, integer *, real *, real *, integer *, integer *
  33 , real *, integer *, integer *), slamrg_(integer *,
  34             integer *, real *, integer *, integer *, integer *);
  35     real orgnrm;
  36     integer coltyp;
  37
  38
  39 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  40 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  41 /*     November 2006 */
  42
  43 /*     .. Scalar Arguments .. */
  44 /*     .. */
  45 /*     .. Array Arguments .. */
  46 /*     .. */
  47
  48 /*  Purpose */
  49 /*  ======= */
  50
  51 /*  SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, */
  52 /*  where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0. */
  53
  54 /*  A related subroutine SLASD7 handles the case in which the singular */
  55 /*  values (and the singular vectors in factored form) are desired. */
  56
  57 /*  SLASD1 computes the SVD as follows: */
  58
  59 /*                ( D1(in)  0    0     0 ) */
  60 /*    B = U(in) * (   Z1'   a   Z2'    b ) * VT(in) */
  61 /*                (   0     0   D2(in) 0 ) */
  62
  63 /*      = U(out) * ( D(out) 0) * VT(out) */
  64
  65 /*  where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M */
  66 /*  with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros */
  67 /*  elsewhere; and the entry b is empty if SQRE = 0. */
  68
  69 /*  The left singular vectors of the original matrix are stored in U, and */
  70 /*  the transpose of the right singular vectors are stored in VT, and the */
  71 /*  singular values are in D.  The algorithm consists of three stages: */
  72
  73 /*     The first stage consists of deflating the size of the problem */
  74 /*     when there are multiple singular values or when there are zeros in */
  75 /*     the Z vector.  For each such occurence the dimension of the */
  76 /*     secular equation problem is reduced by one.  This stage is */
  77 /*     performed by the routine SLASD2. */
  78
  79 /*     The second stage consists of calculating the updated */
  80 /*     singular values. This is done by finding the square roots of the */
  81 /*     roots of the secular equation via the routine SLASD4 (as called */
  82 /*     by SLASD3). This routine also calculates the singular vectors of */
  83 /*     the current problem. */
  84
  85 /*     The final stage consists of computing the updated singular vectors */
  86 /*     directly using the updated singular values.  The singular vectors */
  87 /*     for the current problem are multiplied with the singular vectors */
  88 /*     from the overall problem. */
  89
  90 /*  Arguments */
  91 /*  ========= */
  92
  93 /*  NL     (input) INTEGER */
  94 /*         The row dimension of the upper block.  NL >= 1. */
  95
  96 /*  NR     (input) INTEGER */
  97 /*         The row dimension of the lower block.  NR >= 1. */
  98
  99 /*  SQRE   (input) INTEGER */
 100 /*         = 0: the lower block is an NR-by-NR square matrix. */
 101 /*         = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
 102
 103 /*         The bidiagonal matrix has row dimension N = NL + NR + 1, */
 104 /*         and column dimension M = N + SQRE. */
 105
 106 /*  D      (input/output) REAL array, dimension (NL+NR+1). */
 107 /*         N = NL+NR+1 */
 108 /*         On entry D(1:NL,1:NL) contains the singular values of the */
 109 /*         upper block; and D(NL+2:N) contains the singular values of */
 110 /*         the lower block. On exit D(1:N) contains the singular values */
 111 /*         of the modified matrix. */
 112
 113 /*  ALPHA  (input/output) REAL */
 114 /*         Contains the diagonal element associated with the added row. */
 115
 116 /*  BETA   (input/output) REAL */
 117 /*         Contains the off-diagonal element associated with the added */
 118 /*         row. */
 119
 120 /*  U      (input/output) REAL array, dimension (LDU,N) */
 121 /*         On entry U(1:NL, 1:NL) contains the left singular vectors of */
 122 /*         the upper block; U(NL+2:N, NL+2:N) contains the left singular */
 123 /*         vectors of the lower block. On exit U contains the left */
 124 /*         singular vectors of the bidiagonal matrix. */
 125
 126 /*  LDU    (input) INTEGER */
 127 /*         The leading dimension of the array U.  LDU >= max( 1, N ). */
 128
 129 /*  VT     (input/output) REAL array, dimension (LDVT,M) */
 130 /*         where M = N + SQRE. */
 131 /*         On entry VT(1:NL+1, 1:NL+1)' contains the right singular */
 132 /*         vectors of the upper block; VT(NL+2:M, NL+2:M)' contains */
 133 /*         the right singular vectors of the lower block. On exit */
 134 /*         VT' contains the right singular vectors of the */
 135 /*         bidiagonal matrix. */
 136
 137 /*  LDVT   (input) INTEGER */
 138 /*         The leading dimension of the array VT.  LDVT >= max( 1, M ). */
 139
 140 /*  IDXQ  (output) INTEGER array, dimension (N) */
 141 /*         This contains the permutation which will reintegrate the */
 142 /*         subproblem just solved back into sorted order, i.e. */
 143 /*         D( IDXQ( I = 1, N ) ) will be in ascending order. */
 144
 145 /*  IWORK  (workspace) INTEGER array, dimension (4*N) */
 146
 147 /*  WORK   (workspace) REAL array, dimension (3*M**2+2*M) */
 148
 149 /*  INFO   (output) INTEGER */
 150 /*          = 0:  successful exit. */
 151 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 152 /*          > 0:  if INFO = 1, an singular value did not converge */
 153
 154 /*  Further Details */
 155 /*  =============== */
 156
 157 /*  Based on contributions by */
 158 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
 159 /*     California at Berkeley, USA */
 160
 161 /*  ===================================================================== */
 162
 163 /*     .. Parameters .. */
 164
 165 /*     .. */
 166 /*     .. Local Scalars .. */
 167 /*     .. */
 168 /*     .. External Subroutines .. */
 169 /*     .. */
 170 /*     .. Intrinsic Functions .. */
 171 /*     .. */
 172 /*     .. Executable Statements .. */
 173
 174 /*     Test the input parameters. */
 175
 176     /* Parameter adjustments */
 177     --d__;
 178     u_dim1 = *ldu;
 179     u_offset = 1 + u_dim1;
 180     u -= u_offset;
 181     vt_dim1 = *ldvt;
 182     vt_offset = 1 + vt_dim1;
 183     vt -= vt_offset;
 184     --idxq;
 185     --iwork;
 186     --work;
 187
 188     /* Function Body */
 189     *info = 0;
 190
 191     if (*nl < 1) {
 192         *info = -1;
 193     } else if (*nr < 1) {
 194         *info = -2;
 195     } else if (*sqre < 0 || *sqre > 1) {
 196         *info = -3;
 197     }
 198     if (*info != 0) {
 199         i__1 = -(*info);
 200         xerbla_("SLASD1", &i__1);
 201         return 0;
 202     }
 203
 204     n = *nl + *nr + 1;
 205     m = n + *sqre;
 206
 207 /*     The following values are for bookkeeping purposes only.  They are */
 208 /*     integer pointers which indicate the portion of the workspace */
 209 /*     used by a particular array in SLASD2 and SLASD3. */
 210
 211     ldu2 = n;
 212     ldvt2 = m;
 213
 214     iz = 1;
 215     isigma = iz + m;
 216     iu2 = isigma + n;
 217     ivt2 = iu2 + ldu2 * n;
 218     iq = ivt2 + ldvt2 * m;
 219
 220     idx = 1;
 221     idxc = idx + n;
 222     coltyp = idxc + n;
 223     idxp = coltyp + n;
 224
 225 /*     Scale. */
 226
 227 /* Computing MAX */
 228     r__1 = dabs(*alpha), r__2 = dabs(*beta);
 229     orgnrm = dmax(r__1,r__2);
 230     d__[*nl + 1] = 0.f;
 231     i__1 = n;
 232     for (i__ = 1; i__ <= i__1; ++i__) {
 233         if ((r__1 = d__[i__], dabs(r__1)) > orgnrm) {
 234             orgnrm = (r__1 = d__[i__], dabs(r__1));
 235         }
 236 /* L10: */
 237     }
 238     slascl_("G", &c__0, &c__0, &orgnrm, &c_b7, &n, &c__1, &d__[1], &n, info);
 239     *alpha /= orgnrm;
 240     *beta /= orgnrm;
 241
 242 /*     Deflate singular values. */
 243
 244     slasd2_(nl, nr, sqre, &k, &d__[1], &work[iz], alpha, beta, &u[u_offset],
 245             ldu, &vt[vt_offset], ldvt, &work[isigma], &work[iu2], &ldu2, &
 246             work[ivt2], &ldvt2, &iwork[idxp], &iwork[idx], &iwork[idxc], &
 247             idxq[1], &iwork[coltyp], info);
 248
 249 /*     Solve Secular Equation and update singular vectors. */
 250
 251     ldq = k;
 252     slasd3_(nl, nr, sqre, &k, &d__[1], &work[iq], &ldq, &work[isigma], &u[
 253             u_offset], ldu, &work[iu2], &ldu2, &vt[vt_offset], ldvt, &work[
 254             ivt2], &ldvt2, &iwork[idxc], &iwork[coltyp], &work[iz], info);
 255     if (*info != 0) {
 256         return 0;
 257     }
 258
 259 /*     Unscale. */
 260
 261     slascl_("G", &c__0, &c__0, &c_b7, &orgnrm, &n, &c__1, &d__[1], &n, info);
 262
 263 /*     Prepare the IDXQ sorting permutation. */
 264
 265     n1 = k;
 266     n2 = n - k;
 267     slamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);
 268
 269     return 0;
 270
 271 /*     End of SLASD1 */
 272
 273 } /* slasd1_ */