slasd8.c
Go to the documentation of this file.
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_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:11