src/gmxlib/gmx_lapack/dlasd8.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38
  39 void
  40 F77_FUNC(dlasd8,DLASD8)(int *icompq,
  41         int *k,
  42         double *d__,
  43         double *z__,
  44         double *vf,
  45         double *vl,
  46         double *difl,
  47         double *difr,
  48         int *lddifr,
  49         double *dsigma,
  50         double *work,
  51         int *info)
  52 {
  53     int difr_dim1, difr_offset, i__1, i__2;
  54     double d__2;
  55     double *p1,*p2,t1,t2;
  56
  57     int i__, j;
  58     double dj, rho;
  59     int iwk1, iwk2, iwk3;
  60     double temp;
  61     int iwk2i, iwk3i;
  62     double diflj, difrj, dsigj;
  63     double dsigjp;
  64     int c__1 = 1;
  65     int c__0 = 0;
  66     double one = 1.;
  67
  68     /* avoid warnings on high gcc optimization levels */
  69     difrj = dsigjp = 0;
  70
  71      --d__;
  72     --z__;
  73     --vf;
  74     --vl;
  75     --difl;
  76     difr_dim1 = *lddifr;
  77     difr_offset = 1 + difr_dim1;
  78     difr -= difr_offset;
  79     --dsigma;
  80     --work;
  81
  82     *info = 0;
  83
  84     p1 = &t1;
  85     p2 = &t2;
  86
  87     if (*k == 1) {
  88         d__[1] = fabs(z__[1]);
  89         difl[1] = d__[1];
  90         if (*icompq == 1) {
  91             difl[2] = 1.;
  92             difr[(difr_dim1 << 1) + 1] = 1.;
  93         }
  94         return;
  95     }
  96
  97     i__1 = *k;
  98     for (i__ = 1; i__ <= i__1; ++i__) {
  99       t1 = dsigma[i__];
 100       t2 = dsigma[i__];
 101       /* force store and reload from memory */
 102       d__2 = (*p1) + (*p2) - dsigma[i__];
 103     }
 104
 105     iwk1 = 1;
 106     iwk2 = iwk1 + *k;
 107     iwk3 = iwk2 + *k;
 108     iwk2i = iwk2 - 1;
 109     iwk3i = iwk3 - 1;
 110
 111     rho = F77_FUNC(dnrm2,DNRM2)(k, &z__[1], &c__1);
 112     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &rho, &one, k, &c__1, &z__[1], k, info);
 113     rho *= rho;
 114
 115     F77_FUNC(dlaset,DLASET)("A", k, &c__1, &one, &one, &work[iwk3], k);
 116
 117     i__1 = *k;
 118     for (j = 1; j <= i__1; ++j) {
 119         F77_FUNC(dlasd4,DLASD4)(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[
 120                 iwk2], info);
 121
 122         if (*info != 0) {
 123             return;
 124         }
 125         work[iwk3i + j] = work[iwk3i + j] * work[j] * work[iwk2i + j];
 126         difl[j] = -work[j];
 127         difr[j + difr_dim1] = -work[j + 1];
 128         i__2 = j - 1;
 129         for (i__ = 1; i__ <= i__2; ++i__) {
 130             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
 131                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
 132                     j]);
 133         }
 134         i__2 = *k;
 135         for (i__ = j + 1; i__ <= i__2; ++i__) {
 136             work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
 137                     i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
 138                     j]);
 139         }
 140     }
 141
 142     i__1 = *k;
 143     for (i__ = 1; i__ <= i__1; ++i__) {
 144         d__2 = sqrt(fabs(work[iwk3i + i__]));
 145         z__[i__] = (z__[i__] > 0) ? d__2 : -d__2;
 146     }
 147
 148     i__1 = *k;
 149     for (j = 1; j <= i__1; ++j) {
 150         diflj = difl[j];
 151         dj = d__[j];
 152         dsigj = -dsigma[j];
 153         if (j < *k) {
 154             difrj = -difr[j + difr_dim1];
 155             dsigjp = -dsigma[j + 1];
 156         }
 157         work[j] = -z__[j] / diflj / (dsigma[j] + dj);
 158         i__2 = j - 1;
 159         for (i__ = 1; i__ <= i__2; ++i__) {
 160           t1 = dsigma[i__];
 161           t2 = dsigj;
 162           /* force store and reload from memory */
 163           t1 = (*p1) + (*p2) - diflj;
 164           work[i__] = z__[i__] / t1 / ( dsigma[i__] + dj);
 165         }
 166         i__2 = *k;
 167         for (i__ = j + 1; i__ <= i__2; ++i__) {
 168           t1 = dsigma[i__];
 169           t2 = dsigjp;
 170           /* force store and reload from memory */
 171           t1 = (*p1) + (*p2) - difrj;
 172             work[i__] = z__[i__] / t1 / (dsigma[i__] + dj);
 173         }
 174         temp = F77_FUNC(dnrm2,DNRM2)(k, &work[1], &c__1);
 175         work[iwk2i + j] = F77_FUNC(ddot,DDOT)(k, &work[1], &c__1, &vf[1], &c__1) / temp;
 176         work[iwk3i + j] = F77_FUNC(ddot,DDOT)(k, &work[1], &c__1, &vl[1], &c__1) / temp;
 177         if (*icompq == 1) {
 178             difr[j + (difr_dim1 << 1)] = temp;
 179         }
 180     }
 181
 182     F77_FUNC(dcopy,DCOPY)(k, &work[iwk2], &c__1, &vf[1], &c__1);
 183     F77_FUNC(dcopy,DCOPY)(k, &work[iwk3], &c__1, &vl[1], &c__1);
 184
 185     return;
 186
 187 }