src/gmxlib/gmx_lapack/dlasd8.c

   1 #include <math.h>
   2 #include "gmx_blas.h"
   3 #include "gmx_lapack.h"
   4
   5 void
   6 F77_FUNC(dlasd8,DLASD8)(int *icompq,
   7         int *k,
   8         double *d__,
   9         double *z__,
  10         double *vf,
  11         double *vl,
  12         double *difl,
  13         double *difr,
  14         int *lddifr,
  15         double *dsigma,
  16         double *work,
  17         int *info)
  18 {
  19     int difr_dim1, difr_offset, i__1, i__2;
  20     double d__2;
  21     double *p1,*p2,t1,t2;
  22
  23     int i__, j;
  24     double dj, rho;
  25     int iwk1, iwk2, iwk3;
  26     double temp;
  27     int iwk2i, iwk3i;
  28     double diflj, difrj, dsigj;
  29     double dsigjp;
  30     int c__1 = 1;
  31     int c__0 = 0;
  32     double one = 1.;
  33
  34     /* avoid warnings on high gcc optimization levels */
  35     difrj = dsigjp = 0;
  36
  37      --d__;
  38     --z__;
  39     --vf;
  40     --vl;
  41     --difl;
  42     difr_dim1 = *lddifr;
  43     difr_offset = 1 + difr_dim1;
  44     difr -= difr_offset;
  45     --dsigma;
  46     --work;
  47
  48     *info = 0;
  49
  50     p1 = &t1;
  51     p2 = &t2;
  52
  53     if (*k == 1) {
  54         d__[1] = fabs(z__[1]);
  55         difl[1] = d__[1];
  56         if (*icompq == 1) {
  57             difl[2] = 1.;
  58             difr[(difr_dim1 << 1) + 1] = 1.;
  59         }
  60         return;
  61     }
  62
  63     i__1 = *k;
  64     for (i__ = 1; i__ <= i__1; ++i__) {
  65       t1 = dsigma[i__];
  66       t2 = dsigma[i__];
  67       /* force store and reload from memory */
  68       d__2 = (*p1) + (*p2) - dsigma[i__];
  69     }
  70
  71     iwk1 = 1;
  72     iwk2 = iwk1 + *k;
  73     iwk3 = iwk2 + *k;
  74     iwk2i = iwk2 - 1;
  75     iwk3i = iwk3 - 1;
  76
  77     rho = F77_FUNC(dnrm2,DNRM2)(k, &z__[1], &c__1);
  78     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &rho, &one, k, &c__1, &z__[1], k, info);
  79     rho *= rho;
  80
  81     F77_FUNC(dlaset,DLASET)("A", k, &c__1, &one, &one, &work[iwk3], k);
  82
  83     i__1 = *k;
  84     for (j = 1; j <= i__1; ++j) {
  85         F77_FUNC(dlasd4,DLASD4)(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[
  86                 iwk2], info);
  87
  88         if (*info != 0) {
  89             return;
  90         }
  91         work[iwk3i + j] = work[iwk3i + j] * work[j] * work[iwk2i + j];
  92         difl[j] = -work[j];
  93         difr[j + difr_dim1] = -work[j + 1];
  94         i__2 = j - 1;
  95         for (i__ = 1; i__ <= i__2; ++i__) {
  96             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
  97                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
  98                     j]);
  99         }
 100         i__2 = *k;
 101         for (i__ = j + 1; i__ <= i__2; ++i__) {
 102             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
 103                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
 104                     j]);
 105         }
 106     }
 107
 108     i__1 = *k;
 109     for (i__ = 1; i__ <= i__1; ++i__) {
 110         d__2 = sqrt(fabs(work[iwk3i + i__]));
 111         z__[i__] = (z__[i__] > 0) ? d__2 : -d__2;
 112     }
 113
 114     i__1 = *k;
 115     for (j = 1; j <= i__1; ++j) {
 116         diflj = difl[j];
 117         dj = d__[j];
 118         dsigj = -dsigma[j];
 119         if (j < *k) {
 120             difrj = -difr[j + difr_dim1];
 121             dsigjp = -dsigma[j + 1];
 122         }
 123         work[j] = -z__[j] / diflj / (dsigma[j] + dj);
 124         i__2 = j - 1;
 125         for (i__ = 1; i__ <= i__2; ++i__) {
 126           t1 = dsigma[i__];
 127           t2 = dsigj;
 128           /* force store and reload from memory */
 129           t1 = (*p1) + (*p2) - diflj;
 130           work[i__] = z__[i__] / t1 / ( dsigma[i__] + dj);
 131         }
 132         i__2 = *k;
 133         for (i__ = j + 1; i__ <= i__2; ++i__) {
 134           t1 = dsigma[i__];
 135           t2 = dsigjp;
 136           /* force store and reload from memory */
 137           t1 = (*p1) + (*p2) - difrj;
 138             work[i__] = z__[i__] / t1 / (dsigma[i__] + dj);
 139         }
 140         temp = F77_FUNC(dnrm2,DNRM2)(k, &work[1], &c__1);
 141         work[iwk2i + j] = F77_FUNC(ddot,DDOT)(k, &work[1], &c__1, &vf[1], &c__1) / temp;
 142         work[iwk3i + j] = F77_FUNC(ddot,DDOT)(k, &work[1], &c__1, &vl[1], &c__1) / temp;
 143         if (*icompq == 1) {
 144             difr[j + (difr_dim1 << 1)] = temp;
 145         }
 146     }
 147
 148     F77_FUNC(dcopy,DCOPY)(k, &work[iwk2], &c__1, &vf[1], &c__1);
 149     F77_FUNC(dcopy,DCOPY)(k, &work[iwk3], &c__1, &vl[1], &c__1);
 150
 151     return;
 152
 153 }