00001 /* slasd1.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__0 = 0; 00019 static real c_b7 = 1.f; 00020 static integer c__1 = 1; 00021 static integer c_n1 = -1; 00022 00023 /* Subroutine */ int slasd1_(integer *nl, integer *nr, integer *sqre, real * 00024 d__, real *alpha, real *beta, real *u, integer *ldu, real *vt, 00025 integer *ldvt, integer *idxq, integer *iwork, real *work, integer * 00026 info) 00027 { 00028 /* System generated locals */ 00029 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1; 00030 real r__1, r__2; 00031 00032 /* Local variables */ 00033 integer i__, k, m, n, n1, n2, iq, iz, iu2, ldq, idx, ldu2, ivt2, idxc, 00034 idxp, ldvt2; 00035 extern /* Subroutine */ int slasd2_(integer *, integer *, integer *, 00036 integer *, real *, real *, real *, real *, real *, integer *, 00037 real *, integer *, real *, real *, integer *, real *, integer *, 00038 integer *, integer *, integer *, integer *, integer *, integer *), 00039 slasd3_(integer *, integer *, integer *, integer *, real *, real 00040 *, integer *, real *, real *, integer *, real *, integer *, real * 00041 , integer *, real *, integer *, integer *, integer *, real *, 00042 integer *); 00043 integer isigma; 00044 extern /* Subroutine */ int xerbla_(char *, integer *), slascl_( 00045 char *, integer *, integer *, real *, real *, integer *, integer * 00046 , real *, integer *, integer *), slamrg_(integer *, 00047 integer *, real *, integer *, integer *, integer *); 00048 real orgnrm; 00049 integer coltyp; 00050 00051 00052 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00053 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00054 /* November 2006 */ 00055 00056 /* .. Scalar Arguments .. */ 00057 /* .. */ 00058 /* .. Array Arguments .. */ 00059 /* .. */ 00060 00061 /* Purpose */ 00062 /* ======= */ 00063 00064 /* SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, */ 00065 /* where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0. */ 00066 00067 /* A related subroutine SLASD7 handles the case in which the singular */ 00068 /* values (and the singular vectors in factored form) are desired. */ 00069 00070 /* SLASD1 computes the SVD as follows: */ 00071 00072 /* ( D1(in) 0 0 0 ) */ 00073 /* B = U(in) * ( Z1' a Z2' b ) * VT(in) */ 00074 /* ( 0 0 D2(in) 0 ) */ 00075 00076 /* = U(out) * ( D(out) 0) * VT(out) */ 00077 00078 /* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M */ 00079 /* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros */ 00080 /* elsewhere; and the entry b is empty if SQRE = 0. */ 00081 00082 /* The left singular vectors of the original matrix are stored in U, and */ 00083 /* the transpose of the right singular vectors are stored in VT, and the */ 00084 /* singular values are in D. The algorithm consists of three stages: */ 00085 00086 /* The first stage consists of deflating the size of the problem */ 00087 /* when there are multiple singular values or when there are zeros in */ 00088 /* the Z vector. For each such occurence the dimension of the */ 00089 /* secular equation problem is reduced by one. This stage is */ 00090 /* performed by the routine SLASD2. */ 00091 00092 /* The second stage consists of calculating the updated */ 00093 /* singular values. This is done by finding the square roots of the */ 00094 /* roots of the secular equation via the routine SLASD4 (as called */ 00095 /* by SLASD3). This routine also calculates the singular vectors of */ 00096 /* the current problem. */ 00097 00098 /* The final stage consists of computing the updated singular vectors */ 00099 /* directly using the updated singular values. The singular vectors */ 00100 /* for the current problem are multiplied with the singular vectors */ 00101 /* from the overall problem. */ 00102 00103 /* Arguments */ 00104 /* ========= */ 00105 00106 /* NL (input) INTEGER */ 00107 /* The row dimension of the upper block. NL >= 1. */ 00108 00109 /* NR (input) INTEGER */ 00110 /* The row dimension of the lower block. NR >= 1. */ 00111 00112 /* SQRE (input) INTEGER */ 00113 /* = 0: the lower block is an NR-by-NR square matrix. */ 00114 /* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */ 00115 00116 /* The bidiagonal matrix has row dimension N = NL + NR + 1, */ 00117 /* and column dimension M = N + SQRE. */ 00118 00119 /* D (input/output) REAL array, dimension (NL+NR+1). */ 00120 /* N = NL+NR+1 */ 00121 /* On entry D(1:NL,1:NL) contains the singular values of the */ 00122 /* upper block; and D(NL+2:N) contains the singular values of */ 00123 /* the lower block. On exit D(1:N) contains the singular values */ 00124 /* of the modified matrix. */ 00125 00126 /* ALPHA (input/output) REAL */ 00127 /* Contains the diagonal element associated with the added row. */ 00128 00129 /* BETA (input/output) REAL */ 00130 /* Contains the off-diagonal element associated with the added */ 00131 /* row. */ 00132 00133 /* U (input/output) REAL array, dimension (LDU,N) */ 00134 /* On entry U(1:NL, 1:NL) contains the left singular vectors of */ 00135 /* the upper block; U(NL+2:N, NL+2:N) contains the left singular */ 00136 /* vectors of the lower block. On exit U contains the left */ 00137 /* singular vectors of the bidiagonal matrix. */ 00138 00139 /* LDU (input) INTEGER */ 00140 /* The leading dimension of the array U. LDU >= max( 1, N ). */ 00141 00142 /* VT (input/output) REAL array, dimension (LDVT,M) */ 00143 /* where M = N + SQRE. */ 00144 /* On entry VT(1:NL+1, 1:NL+1)' contains the right singular */ 00145 /* vectors of the upper block; VT(NL+2:M, NL+2:M)' contains */ 00146 /* the right singular vectors of the lower block. On exit */ 00147 /* VT' contains the right singular vectors of the */ 00148 /* bidiagonal matrix. */ 00149 00150 /* LDVT (input) INTEGER */ 00151 /* The leading dimension of the array VT. LDVT >= max( 1, M ). */ 00152 00153 /* IDXQ (output) INTEGER array, dimension (N) */ 00154 /* This contains the permutation which will reintegrate the */ 00155 /* subproblem just solved back into sorted order, i.e. */ 00156 /* D( IDXQ( I = 1, N ) ) will be in ascending order. */ 00157 00158 /* IWORK (workspace) INTEGER array, dimension (4*N) */ 00159 00160 /* WORK (workspace) REAL array, dimension (3*M**2+2*M) */ 00161 00162 /* INFO (output) INTEGER */ 00163 /* = 0: successful exit. */ 00164 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00165 /* > 0: if INFO = 1, an singular value did not converge */ 00166 00167 /* Further Details */ 00168 /* =============== */ 00169 00170 /* Based on contributions by */ 00171 /* Ming Gu and Huan Ren, Computer Science Division, University of */ 00172 /* California at Berkeley, USA */ 00173 00174 /* ===================================================================== */ 00175 00176 /* .. Parameters .. */ 00177 00178 /* .. */ 00179 /* .. Local Scalars .. */ 00180 /* .. */ 00181 /* .. External Subroutines .. */ 00182 /* .. */ 00183 /* .. Intrinsic Functions .. */ 00184 /* .. */ 00185 /* .. Executable Statements .. */ 00186 00187 /* Test the input parameters. */ 00188 00189 /* Parameter adjustments */ 00190 --d__; 00191 u_dim1 = *ldu; 00192 u_offset = 1 + u_dim1; 00193 u -= u_offset; 00194 vt_dim1 = *ldvt; 00195 vt_offset = 1 + vt_dim1; 00196 vt -= vt_offset; 00197 --idxq; 00198 --iwork; 00199 --work; 00200 00201 /* Function Body */ 00202 *info = 0; 00203 00204 if (*nl < 1) { 00205 *info = -1; 00206 } else if (*nr < 1) { 00207 *info = -2; 00208 } else if (*sqre < 0 || *sqre > 1) { 00209 *info = -3; 00210 } 00211 if (*info != 0) { 00212 i__1 = -(*info); 00213 xerbla_("SLASD1", &i__1); 00214 return 0; 00215 } 00216 00217 n = *nl + *nr + 1; 00218 m = n + *sqre; 00219 00220 /* The following values are for bookkeeping purposes only. They are */ 00221 /* integer pointers which indicate the portion of the workspace */ 00222 /* used by a particular array in SLASD2 and SLASD3. */ 00223 00224 ldu2 = n; 00225 ldvt2 = m; 00226 00227 iz = 1; 00228 isigma = iz + m; 00229 iu2 = isigma + n; 00230 ivt2 = iu2 + ldu2 * n; 00231 iq = ivt2 + ldvt2 * m; 00232 00233 idx = 1; 00234 idxc = idx + n; 00235 coltyp = idxc + n; 00236 idxp = coltyp + n; 00237 00238 /* Scale. */ 00239 00240 /* Computing MAX */ 00241 r__1 = dabs(*alpha), r__2 = dabs(*beta); 00242 orgnrm = dmax(r__1,r__2); 00243 d__[*nl + 1] = 0.f; 00244 i__1 = n; 00245 for (i__ = 1; i__ <= i__1; ++i__) { 00246 if ((r__1 = d__[i__], dabs(r__1)) > orgnrm) { 00247 orgnrm = (r__1 = d__[i__], dabs(r__1)); 00248 } 00249 /* L10: */ 00250 } 00251 slascl_("G", &c__0, &c__0, &orgnrm, &c_b7, &n, &c__1, &d__[1], &n, info); 00252 *alpha /= orgnrm; 00253 *beta /= orgnrm; 00254 00255 /* Deflate singular values. */ 00256 00257 slasd2_(nl, nr, sqre, &k, &d__[1], &work[iz], alpha, beta, &u[u_offset], 00258 ldu, &vt[vt_offset], ldvt, &work[isigma], &work[iu2], &ldu2, & 00259 work[ivt2], &ldvt2, &iwork[idxp], &iwork[idx], &iwork[idxc], & 00260 idxq[1], &iwork[coltyp], info); 00261 00262 /* Solve Secular Equation and update singular vectors. */ 00263 00264 ldq = k; 00265 slasd3_(nl, nr, sqre, &k, &d__[1], &work[iq], &ldq, &work[isigma], &u[ 00266 u_offset], ldu, &work[iu2], &ldu2, &vt[vt_offset], ldvt, &work[ 00267 ivt2], &ldvt2, &iwork[idxc], &iwork[coltyp], &work[iz], info); 00268 if (*info != 0) { 00269 return 0; 00270 } 00271 00272 /* Unscale. */ 00273 00274 slascl_("G", &c__0, &c__0, &c_b7, &orgnrm, &n, &c__1, &d__[1], &n, info); 00275 00276 /* Prepare the IDXQ sorting permutation. */ 00277 00278 n1 = k; 00279 n2 = n - k; 00280 slamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]); 00281 00282 return 0; 00283 00284 /* End of SLASD1 */ 00285 00286 } /* slasd1_ */