00001 /* slasd8.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__1 = 1; 00019 static integer c__0 = 0; 00020 static real c_b8 = 1.f; 00021 00022 /* Subroutine */ int slasd8_(integer *icompq, integer *k, real *d__, real * 00023 z__, real *vf, real *vl, real *difl, real *difr, integer *lddifr, 00024 real *dsigma, real *work, integer *info) 00025 { 00026 /* System generated locals */ 00027 integer difr_dim1, difr_offset, i__1, i__2; 00028 real r__1, r__2; 00029 00030 /* Builtin functions */ 00031 double sqrt(doublereal), r_sign(real *, real *); 00032 00033 /* Local variables */ 00034 integer i__, j; 00035 real dj, rho; 00036 integer iwk1, iwk2, iwk3; 00037 real temp; 00038 extern doublereal sdot_(integer *, real *, integer *, real *, integer *); 00039 integer iwk2i, iwk3i; 00040 extern doublereal snrm2_(integer *, real *, integer *); 00041 real diflj, difrj, dsigj; 00042 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 00043 integer *); 00044 extern doublereal slamc3_(real *, real *); 00045 extern /* Subroutine */ int slasd4_(integer *, integer *, real *, real *, 00046 real *, real *, real *, real *, integer *), xerbla_(char *, 00047 integer *); 00048 real dsigjp; 00049 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 00050 real *, integer *, integer *, real *, integer *, integer *), slaset_(char *, integer *, integer *, real *, real *, 00051 real *, integer *); 00052 00053 00054 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00055 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00056 /* October 2006 */ 00057 00058 /* .. Scalar Arguments .. */ 00059 /* .. */ 00060 /* .. Array Arguments .. */ 00061 /* .. */ 00062 00063 /* Purpose */ 00064 /* ======= */ 00065 00066 /* SLASD8 finds the square roots of the roots of the secular equation, */ 00067 /* as defined by the values in DSIGMA and Z. It makes the appropriate */ 00068 /* calls to SLASD4, and stores, for each element in D, the distance */ 00069 /* to its two nearest poles (elements in DSIGMA). It also updates */ 00070 /* the arrays VF and VL, the first and last components of all the */ 00071 /* right singular vectors of the original bidiagonal matrix. */ 00072 00073 /* SLASD8 is called from SLASD6. */ 00074 00075 /* Arguments */ 00076 /* ========= */ 00077 00078 /* ICOMPQ (input) INTEGER */ 00079 /* Specifies whether singular vectors are to be computed in */ 00080 /* factored form in the calling routine: */ 00081 /* = 0: Compute singular values only. */ 00082 /* = 1: Compute singular vectors in factored form as well. */ 00083 00084 /* K (input) INTEGER */ 00085 /* The number of terms in the rational function to be solved */ 00086 /* by SLASD4. K >= 1. */ 00087 00088 /* D (output) REAL array, dimension ( K ) */ 00089 /* On output, D contains the updated singular values. */ 00090 00091 /* Z (input/output) REAL array, dimension ( K ) */ 00092 /* On entry, the first K elements of this array contain the */ 00093 /* components of the deflation-adjusted updating row vector. */ 00094 /* On exit, Z is updated. */ 00095 00096 /* VF (input/output) REAL array, dimension ( K ) */ 00097 /* On entry, VF contains information passed through DBEDE8. */ 00098 /* On exit, VF contains the first K components of the first */ 00099 /* components of all right singular vectors of the bidiagonal */ 00100 /* matrix. */ 00101 00102 /* VL (input/output) REAL array, dimension ( K ) */ 00103 /* On entry, VL contains information passed through DBEDE8. */ 00104 /* On exit, VL contains the first K components of the last */ 00105 /* components of all right singular vectors of the bidiagonal */ 00106 /* matrix. */ 00107 00108 /* DIFL (output) REAL array, dimension ( K ) */ 00109 /* On exit, DIFL(I) = D(I) - DSIGMA(I). */ 00110 00111 /* DIFR (output) REAL array, */ 00112 /* dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and */ 00113 /* dimension ( K ) if ICOMPQ = 0. */ 00114 /* On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not */ 00115 /* defined and will not be referenced. */ 00116 00117 /* If ICOMPQ = 1, DIFR(1:K,2) is an array containing the */ 00118 /* normalizing factors for the right singular vector matrix. */ 00119 00120 /* LDDIFR (input) INTEGER */ 00121 /* The leading dimension of DIFR, must be at least K. */ 00122 00123 /* DSIGMA (input/output) REAL array, dimension ( K ) */ 00124 /* On entry, the first K elements of this array contain the old */ 00125 /* roots of the deflated updating problem. These are the poles */ 00126 /* of the secular equation. */ 00127 /* On exit, the elements of DSIGMA may be very slightly altered */ 00128 /* in value. */ 00129 00130 /* WORK (workspace) REAL array, dimension at least 3 * K */ 00131 00132 /* INFO (output) INTEGER */ 00133 /* = 0: successful exit. */ 00134 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00135 /* > 0: if INFO = 1, an singular value did not converge */ 00136 00137 /* Further Details */ 00138 /* =============== */ 00139 00140 /* Based on contributions by */ 00141 /* Ming Gu and Huan Ren, Computer Science Division, University of */ 00142 /* California at Berkeley, USA */ 00143 00144 /* ===================================================================== */ 00145 00146 /* .. Parameters .. */ 00147 /* .. */ 00148 /* .. Local Scalars .. */ 00149 /* .. */ 00150 /* .. External Subroutines .. */ 00151 /* .. */ 00152 /* .. External Functions .. */ 00153 /* .. */ 00154 /* .. Intrinsic Functions .. */ 00155 /* .. */ 00156 /* .. Executable Statements .. */ 00157 00158 /* Test the input parameters. */ 00159 00160 /* Parameter adjustments */ 00161 --d__; 00162 --z__; 00163 --vf; 00164 --vl; 00165 --difl; 00166 difr_dim1 = *lddifr; 00167 difr_offset = 1 + difr_dim1; 00168 difr -= difr_offset; 00169 --dsigma; 00170 --work; 00171 00172 /* Function Body */ 00173 *info = 0; 00174 00175 if (*icompq < 0 || *icompq > 1) { 00176 *info = -1; 00177 } else if (*k < 1) { 00178 *info = -2; 00179 } else if (*lddifr < *k) { 00180 *info = -9; 00181 } 00182 if (*info != 0) { 00183 i__1 = -(*info); 00184 xerbla_("SLASD8", &i__1); 00185 return 0; 00186 } 00187 00188 /* Quick return if possible */ 00189 00190 if (*k == 1) { 00191 d__[1] = dabs(z__[1]); 00192 difl[1] = d__[1]; 00193 if (*icompq == 1) { 00194 difl[2] = 1.f; 00195 difr[(difr_dim1 << 1) + 1] = 1.f; 00196 } 00197 return 0; 00198 } 00199 00200 /* Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can */ 00201 /* be computed with high relative accuracy (barring over/underflow). */ 00202 /* This is a problem on machines without a guard digit in */ 00203 /* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */ 00204 /* The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), */ 00205 /* which on any of these machines zeros out the bottommost */ 00206 /* bit of DSIGMA(I) if it is 1; this makes the subsequent */ 00207 /* subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation */ 00208 /* occurs. On binary machines with a guard digit (almost all */ 00209 /* machines) it does not change DSIGMA(I) at all. On hexadecimal */ 00210 /* and decimal machines with a guard digit, it slightly */ 00211 /* changes the bottommost bits of DSIGMA(I). It does not account */ 00212 /* for hexadecimal or decimal machines without guard digits */ 00213 /* (we know of none). We use a subroutine call to compute */ 00214 /* 2*DLAMBDA(I) to prevent optimizing compilers from eliminating */ 00215 /* this code. */ 00216 00217 i__1 = *k; 00218 for (i__ = 1; i__ <= i__1; ++i__) { 00219 dsigma[i__] = slamc3_(&dsigma[i__], &dsigma[i__]) - dsigma[i__]; 00220 /* L10: */ 00221 } 00222 00223 /* Book keeping. */ 00224 00225 iwk1 = 1; 00226 iwk2 = iwk1 + *k; 00227 iwk3 = iwk2 + *k; 00228 iwk2i = iwk2 - 1; 00229 iwk3i = iwk3 - 1; 00230 00231 /* Normalize Z. */ 00232 00233 rho = snrm2_(k, &z__[1], &c__1); 00234 slascl_("G", &c__0, &c__0, &rho, &c_b8, k, &c__1, &z__[1], k, info); 00235 rho *= rho; 00236 00237 /* Initialize WORK(IWK3). */ 00238 00239 slaset_("A", k, &c__1, &c_b8, &c_b8, &work[iwk3], k); 00240 00241 /* Compute the updated singular values, the arrays DIFL, DIFR, */ 00242 /* and the updated Z. */ 00243 00244 i__1 = *k; 00245 for (j = 1; j <= i__1; ++j) { 00246 slasd4_(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[ 00247 iwk2], info); 00248 00249 /* If the root finder fails, the computation is terminated. */ 00250 00251 if (*info != 0) { 00252 return 0; 00253 } 00254 work[iwk3i + j] = work[iwk3i + j] * work[j] * work[iwk2i + j]; 00255 difl[j] = -work[j]; 00256 difr[j + difr_dim1] = -work[j + 1]; 00257 i__2 = j - 1; 00258 for (i__ = 1; i__ <= i__2; ++i__) { 00259 work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i + 00260 i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[ 00261 j]); 00262 /* L20: */ 00263 } 00264 i__2 = *k; 00265 for (i__ = j + 1; i__ <= i__2; ++i__) { 00266 work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i + 00267 i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[ 00268 j]); 00269 /* L30: */ 00270 } 00271 /* L40: */ 00272 } 00273 00274 /* Compute updated Z. */ 00275 00276 i__1 = *k; 00277 for (i__ = 1; i__ <= i__1; ++i__) { 00278 r__2 = sqrt((r__1 = work[iwk3i + i__], dabs(r__1))); 00279 z__[i__] = r_sign(&r__2, &z__[i__]); 00280 /* L50: */ 00281 } 00282 00283 /* Update VF and VL. */ 00284 00285 i__1 = *k; 00286 for (j = 1; j <= i__1; ++j) { 00287 diflj = difl[j]; 00288 dj = d__[j]; 00289 dsigj = -dsigma[j]; 00290 if (j < *k) { 00291 difrj = -difr[j + difr_dim1]; 00292 dsigjp = -dsigma[j + 1]; 00293 } 00294 work[j] = -z__[j] / diflj / (dsigma[j] + dj); 00295 i__2 = j - 1; 00296 for (i__ = 1; i__ <= i__2; ++i__) { 00297 work[i__] = z__[i__] / (slamc3_(&dsigma[i__], &dsigj) - diflj) / ( 00298 dsigma[i__] + dj); 00299 /* L60: */ 00300 } 00301 i__2 = *k; 00302 for (i__ = j + 1; i__ <= i__2; ++i__) { 00303 work[i__] = z__[i__] / (slamc3_(&dsigma[i__], &dsigjp) + difrj) / 00304 (dsigma[i__] + dj); 00305 /* L70: */ 00306 } 00307 temp = snrm2_(k, &work[1], &c__1); 00308 work[iwk2i + j] = sdot_(k, &work[1], &c__1, &vf[1], &c__1) / temp; 00309 work[iwk3i + j] = sdot_(k, &work[1], &c__1, &vl[1], &c__1) / temp; 00310 if (*icompq == 1) { 00311 difr[j + (difr_dim1 << 1)] = temp; 00312 } 00313 /* L80: */ 00314 } 00315 00316 scopy_(k, &work[iwk2], &c__1, &vf[1], &c__1); 00317 scopy_(k, &work[iwk3], &c__1, &vl[1], &c__1); 00318 00319 return 0; 00320 00321 /* End of SLASD8 */ 00322 00323 } /* slasd8_ */