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