src/gmxlib/gmx_lapack/slasd3.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(slasd3,SLASD3)(int *nl,
  41         int *nr,
  42         int *sqre,
  43         int *k,
  44         float *d__,
  45         float *q,
  46         int *ldq,
  47         float *dsigma,
  48         float *u,
  49         int *ldu,
  50         float *u2,
  51         int *ldu2,
  52         float *vt,
  53         int *ldvt,
  54         float *vt2,
  55         int *ldvt2,
  56         int *idxc,
  57         int *ctot,
  58         float *z__,
  59         int *info)
  60 {
  61     int q_dim1, q_offset, u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1,
  62             vt_offset, vt2_dim1, vt2_offset, i__1, i__2;
  63     float d__2;
  64
  65     int i__, j, m, n, jc;
  66     float rho;
  67     int nlp1, nlp2, nrp1;
  68     float temp;
  69     int ctemp;
  70     int ktemp;
  71     int c__1 = 1;
  72     int c__0 = 0;
  73     float zero = 0.0;
  74     float one = 1.0;
  75     float *p1,*p2,t1,t2;
  76
  77     p1 = &t1;
  78     p2 = &t2;
  79
  80     --d__;
  81     q_dim1 = *ldq;
  82     q_offset = 1 + q_dim1;
  83     q -= q_offset;
  84     --dsigma;
  85     u_dim1 = *ldu;
  86     u_offset = 1 + u_dim1;
  87     u -= u_offset;
  88     u2_dim1 = *ldu2;
  89     u2_offset = 1 + u2_dim1;
  90     u2 -= u2_offset;
  91     vt_dim1 = *ldvt;
  92     vt_offset = 1 + vt_dim1;
  93     vt -= vt_offset;
  94     vt2_dim1 = *ldvt2;
  95     vt2_offset = 1 + vt2_dim1;
  96     vt2 -= vt2_offset;
  97     --idxc;
  98     --ctot;
  99     --z__;
 100
 101     /* Function Body */
 102     *info = 0;
 103
 104     if (*nl < 1) {
 105         *info = -1;
 106     } else if (*nr < 1) {
 107         *info = -2;
 108     } else if (*sqre != 1 && *sqre != 0) {
 109         *info = -3;
 110     }
 111
 112     n = *nl + *nr + 1;
 113     m = n + *sqre;
 114     nlp1 = *nl + 1;
 115     nlp2 = *nl + 2;
 116
 117     if (*k == 1) {
 118         d__[1] = fabs(z__[1]);
 119         F77_FUNC(scopy,SCOPY)(&m, &vt2[vt2_dim1 + 1], ldvt2, &vt[vt_dim1 + 1], ldvt);
 120         if (z__[1] > 0.) {
 121             F77_FUNC(scopy,SCOPY)(&n, &u2[u2_dim1 + 1], &c__1, &u[u_dim1 + 1], &c__1);
 122         } else {
 123             i__1 = n;
 124             for (i__ = 1; i__ <= i__1; ++i__) {
 125                 u[i__ + u_dim1] = -u2[i__ + u2_dim1];
 126             }
 127         }
 128         return;
 129     }
 130
 131     i__1 = *k;
 132     for (i__ = 1; i__ <= i__1; ++i__) {
 133       t1 = dsigma[i__];
 134       t2 = dsigma[i__];
 135       /* force store and reload from memory */
 136       t1 = (*p1) + (*p2) - dsigma[i__];
 137
 138       dsigma[i__] = t1;
 139     }
 140
 141     F77_FUNC(scopy,SCOPY)(k, &z__[1], &c__1, &q[q_offset], &c__1);
 142
 143     rho = F77_FUNC(snrm2,SNRM2)(k, &z__[1], &c__1);
 144     F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &rho, &one, k, &c__1, &z__[1], k, info);
 145     rho *= rho;
 146
 147
 148     i__1 = *k;
 149     for (j = 1; j <= i__1; ++j) {
 150         F77_FUNC(slasd4,SLASD4)(k, &j, &dsigma[1], &z__[1], &u[j * u_dim1 + 1], &rho, &d__[j],
 151                  &vt[j * vt_dim1 + 1], info);
 152
 153         if (*info != 0) {
 154             return;
 155         }
 156     }
 157
 158     i__1 = *k;
 159     for (i__ = 1; i__ <= i__1; ++i__) {
 160         z__[i__] = u[i__ + *k * u_dim1] * vt[i__ + *k * vt_dim1];
 161         i__2 = i__ - 1;
 162         for (j = 1; j <= i__2; ++j) {
 163             z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
 164                     i__] - dsigma[j]) / (dsigma[i__] + dsigma[j]);
 165         }
 166         i__2 = *k - 1;
 167         for (j = i__; j <= i__2; ++j) {
 168             z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
 169                     i__] - dsigma[j + 1]) / (dsigma[i__] + dsigma[j + 1]);
 170         }
 171         d__2 = sqrt(fabs(z__[i__]));
 172         z__[i__] = (q[i__ + q_dim1] > 0) ? d__2 : -d__2;
 173     }
 174
 175     i__1 = *k;
 176     for (i__ = 1; i__ <= i__1; ++i__) {
 177         vt[i__ * vt_dim1 + 1] = z__[1] / u[i__ * u_dim1 + 1] / vt[i__ *
 178                 vt_dim1 + 1];
 179         u[i__ * u_dim1 + 1] = -1.;
 180         i__2 = *k;
 181         for (j = 2; j <= i__2; ++j) {
 182             vt[j + i__ * vt_dim1] = z__[j] / u[j + i__ * u_dim1] / vt[j + i__
 183                     * vt_dim1];
 184             u[j + i__ * u_dim1] = dsigma[j] * vt[j + i__ * vt_dim1];
 185         }
 186         temp = F77_FUNC(snrm2,SNRM2)(k, &u[i__ * u_dim1 + 1], &c__1);
 187         q[i__ * q_dim1 + 1] = u[i__ * u_dim1 + 1] / temp;
 188         i__2 = *k;
 189         for (j = 2; j <= i__2; ++j) {
 190             jc = idxc[j];
 191             q[j + i__ * q_dim1] = u[jc + i__ * u_dim1] / temp;
 192         }
 193     }
 194
 195     if (*k == 2) {
 196         F77_FUNC(sgemm,SGEMM)("N", "N", &n, k, k, &one, &u2[u2_offset], ldu2, &q[q_offset],
 197                  ldq, &zero, &u[u_offset], ldu);
 198         goto L100;
 199     }
 200     if (ctot[1] > 0) {
 201         F77_FUNC(sgemm,SGEMM)("N", "N", nl, k, &ctot[1], &one, &u2[(u2_dim1 << 1) + 1],
 202                 ldu2, &q[q_dim1 + 2], ldq, &zero, &u[u_dim1 + 1], ldu);
 203         if (ctot[3] > 0) {
 204             ktemp = ctot[1] + 2 + ctot[2];
 205             F77_FUNC(sgemm,SGEMM)("N", "N", nl, k, &ctot[3], &one, &u2[ktemp * u2_dim1 + 1]
 206                     , ldu2, &q[ktemp + q_dim1], ldq, &one, &u[u_dim1 + 1],
 207                     ldu);
 208         }
 209     } else if (ctot[3] > 0) {
 210         ktemp = ctot[1] + 2 + ctot[2];
 211         F77_FUNC(sgemm,SGEMM)("N", "N", nl, k, &ctot[3], &one, &u2[ktemp * u2_dim1 + 1],
 212                 ldu2, &q[ktemp + q_dim1], ldq, &zero, &u[u_dim1 + 1], ldu);
 213     } else {
 214         F77_FUNC(slacpy,SLACPY)("F", nl, k, &u2[u2_offset], ldu2, &u[u_offset], ldu);
 215     }
 216     F77_FUNC(scopy,SCOPY)(k, &q[q_dim1 + 1], ldq, &u[nlp1 + u_dim1], ldu);
 217     ktemp = ctot[1] + 2;
 218     ctemp = ctot[2] + ctot[3];
 219     F77_FUNC(sgemm,SGEMM)("N", "N", nr, k, &ctemp, &one, &u2[nlp2 + ktemp * u2_dim1], ldu2,
 220              &q[ktemp + q_dim1], ldq, &zero, &u[nlp2 + u_dim1], ldu);
 221
 222 L100:
 223     i__1 = *k;
 224     for (i__ = 1; i__ <= i__1; ++i__) {
 225         temp = F77_FUNC(snrm2,SNRM2)(k, &vt[i__ * vt_dim1 + 1], &c__1);
 226         q[i__ + q_dim1] = vt[i__ * vt_dim1 + 1] / temp;
 227         i__2 = *k;
 228         for (j = 2; j <= i__2; ++j) {
 229             jc = idxc[j];
 230             q[i__ + j * q_dim1] = vt[jc + i__ * vt_dim1] / temp;
 231         }
 232     }
 233
 234     if (*k == 2) {
 235         F77_FUNC(sgemm,SGEMM)("N", "N", k, &m, k, &one, &q[q_offset], ldq, &vt2[vt2_offset]
 236                 , ldvt2, &zero, &vt[vt_offset], ldvt);
 237         return;
 238     }
 239     ktemp = ctot[1] + 1;
 240     F77_FUNC(sgemm,SGEMM)("N", "N", k, &nlp1, &ktemp, &one, &q[q_dim1 + 1], ldq, &vt2[
 241             vt2_dim1 + 1], ldvt2, &zero, &vt[vt_dim1 + 1], ldvt);
 242     ktemp = ctot[1] + 2 + ctot[2];
 243     if (ktemp <= *ldvt2) {
 244         F77_FUNC(sgemm,SGEMM)("N", "N", k, &nlp1, &ctot[3], &one, &q[ktemp * q_dim1 + 1],
 245                 ldq, &vt2[ktemp + vt2_dim1], ldvt2, &one, &vt[vt_dim1 + 1],
 246                 ldvt);
 247     }
 248
 249     ktemp = ctot[1] + 1;
 250     nrp1 = *nr + *sqre;
 251     if (ktemp > 1) {
 252         i__1 = *k;
 253         for (i__ = 1; i__ <= i__1; ++i__) {
 254             q[i__ + ktemp * q_dim1] = q[i__ + q_dim1];
 255         }
 256         i__1 = m;
 257         for (i__ = nlp2; i__ <= i__1; ++i__) {
 258             vt2[ktemp + i__ * vt2_dim1] = vt2[i__ * vt2_dim1 + 1];
 259         }
 260     }
 261     ctemp = ctot[2] + 1 + ctot[3];
 262     F77_FUNC(sgemm,SGEMM)("N", "N", k, &nrp1, &ctemp, &one, &q[ktemp * q_dim1 + 1], ldq, &
 263             vt2[ktemp + nlp2 * vt2_dim1], ldvt2, &zero, &vt[nlp2 * vt_dim1 +
 264             1], ldvt);
 265
 266     return;
 267
 268
 269 }
 270
 271