src/gmxlib/gmx_lapack/sbdsdc.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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 <ctype.h>
  36 #include <math.h>
  37 #include "gmx_blas.h"
  38 #include "gmx_lapack.h"
  39 #include "lapack_limits.h"
  40
  41 #include <types/simple.h>
  42
  43 void
  44 F77_FUNC(sbdsdc,SBDSDC)(const char *uplo,
  45         const char *compq,
  46         int *n,
  47         float *d__,
  48         float *e,
  49         float *u,
  50         int *ldu,
  51         float *vt,
  52         int *ldvt,
  53         float *q,
  54         int *iq,
  55         float *work,
  56         int *iwork,
  57         int *info)
  58 {
  59     int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  60     int i__, j, k;
  61     float p, r__;
  62     int z__, ic, ii, kk;
  63     float cs;
  64     int is, iu;
  65     float sn;
  66     int nm1;
  67     float eps;
  68     int ivt, difl, difr, ierr, perm, mlvl, sqre;
  69     int poles, iuplo, nsize, start;
  70     int givcol;
  71     int icompq;
  72     float orgnrm;
  73     int givnum, givptr, qstart, smlsiz, wstart, smlszp;
  74     float zero = 0.0;
  75     float one = 1.0;
  76     int c_0 = 0;
  77     int c_1 = 1;
  78
  79     --d__;
  80     --e;
  81     u_dim1 = *ldu;
  82     u_offset = 1 + u_dim1;
  83     u -= u_offset;
  84     vt_dim1 = *ldvt;
  85     vt_offset = 1 + vt_dim1;
  86     vt -= vt_offset;
  87     --q;
  88     --iq;
  89     --work;
  90     --iwork;
  91
  92     k = iu = z__ = ic = is = ivt = difl = difr = perm = 0;
  93     poles = givnum = givptr = givcol = 0;
  94
  95     smlsiz = DBDSDC_SMALLSIZE;
  96     *info = 0;
  97
  98     iuplo = (*uplo=='U' || *uplo=='u') ? 1 : 2;
  99
 100     switch(*compq) {
 101     case 'n':
 102     case 'N':
 103       icompq = 0;
 104       break;
 105     case 'p':
 106     case 'P':
 107       icompq = 1;
 108       break;
 109     case 'i':
 110     case 'I':
 111       icompq = 2;
 112       break;
 113     default:
 114       return;
 115     }
 116
 117     if (*n <= 0)
 118         return;
 119
 120     if (*n == 1) {
 121         if (icompq == 1) {
 122           q[1] = (d__[1]>0) ? 1.0 : -1.0;
 123           q[smlsiz * *n + 1] = 1.;
 124         } else if (icompq == 2) {
 125           u[u_dim1 + 1] = (d__[1]>0) ? 1.0 : -1.0;
 126           vt[vt_dim1 + 1] = 1.;
 127         }
 128         d__[1] = fabs(d__[1]);
 129         return;
 130     }
 131     nm1 = *n - 1;
 132     wstart = 1;
 133     qstart = 3;
 134     if (icompq == 1) {
 135         F77_FUNC(scopy,SCOPY)(n, &d__[1], &c_1, &q[1], &c_1);
 136         i__1 = *n - 1;
 137         F77_FUNC(scopy,SCOPY)(&i__1, &e[1], &c_1, &q[*n + 1], &c_1);
 138     }
 139     if (iuplo == 2) {
 140         qstart = 5;
 141         wstart = (*n << 1) - 1;
 142         i__1 = *n - 1;
 143         for (i__ = 1; i__ <= i__1; ++i__) {
 144             F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 145             d__[i__] = r__;
 146             e[i__] = sn * d__[i__ + 1];
 147             d__[i__ + 1] = cs * d__[i__ + 1];
 148             if (icompq == 1) {
 149                 q[i__ + (*n << 1)] = cs;
 150                 q[i__ + *n * 3] = sn;
 151             } else if (icompq == 2) {
 152                 work[i__] = cs;
 153                 work[nm1 + i__] = -sn;
 154             }
 155         }
 156     }
 157     if (icompq == 0) {
 158       F77_FUNC(slasdq,SLASDQ)("U",&c_0,n,&c_0,&c_0,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
 159               &u[u_offset], ldu, &u[u_offset], ldu, &work[wstart], info);
 160         goto L40;
 161     }
 162     if (*n <= smlsiz) {
 163         if (icompq == 2) {
 164             F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
 165             F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
 166             F77_FUNC(slasdq,SLASDQ)("U",&c_0,n,n,n,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
 167                     &u[u_offset],ldu,&u[u_offset],ldu,&work[wstart],info);
 168         } else if (icompq == 1) {
 169             iu = 1;
 170             ivt = iu + *n;
 171             F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &q[iu + (qstart - 1) * *n], n);
 172             F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &q[ivt + (qstart - 1) * *n], n);
 173             F77_FUNC(slasdq,SLASDQ)("U", &c_0, n, n, n, &c_0, &d__[1], &e[1],
 174                     &q[ivt + (qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n],
 175                     n, &q[iu + (qstart - 1) * *n], n, &work[wstart], info);
 176         }
 177         goto L40;
 178     }
 179
 180     if (icompq == 2) {
 181         F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
 182         F77_FUNC(slaset,SLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
 183     }
 184
 185     orgnrm = F77_FUNC(slanst,SLANST)("M", n, &d__[1], &e[1]);
 186     if ( fabs(orgnrm)<GMX_FLOAT_MIN) {
 187         return;
 188     }
 189     F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &orgnrm, &one, n, &c_1, &d__[1], n, &ierr);
 190     F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &orgnrm, &one, &nm1, &c_1, &e[1], &nm1, &ierr);
 191
 192     eps = GMX_FLOAT_EPS;
 193
 194     mlvl = (int) (log((float) (*n) / (float) (smlsiz + 1)) /
 195             log(2.)) + 1;
 196     smlszp = smlsiz + 1;
 197
 198     if (icompq == 1) {
 199         iu = 1;
 200         ivt = smlsiz + 1;
 201         difl = ivt + smlszp;
 202         difr = difl + mlvl;
 203         z__ = difr + (mlvl << 1);
 204         ic = z__ + mlvl;
 205         is = ic + 1;
 206         poles = is + 1;
 207         givnum = poles + (mlvl << 1);
 208
 209         k = 1;
 210         givptr = 2;
 211         perm = 3;
 212         givcol = perm + mlvl;
 213     }
 214
 215     i__1 = *n;
 216     for (i__ = 1; i__ <= i__1; ++i__) {
 217         if (fabs(d__[i__]) < eps)
 218             d__[i__] = (d__[i__]>0) ? eps : -eps;
 219     }
 220
 221     start = 1;
 222     sqre = 0;
 223
 224     i__1 = nm1;
 225     for (i__ = 1; i__ <= i__1; ++i__) {
 226         if (fabs(e[i__]) < eps || i__ == nm1) {
 227             if (i__ < nm1) {
 228                 nsize = i__ - start + 1;
 229             } else if (fabs(e[i__]) >= eps) {
 230                 nsize = *n - start + 1;
 231             } else {
 232                 nsize = i__ - start + 1;
 233                 if (icompq == 2) {
 234                     u[*n + *n * u_dim1] = (d__[*n]>0) ? 1.0 : -1.0;
 235                     vt[*n + *n * vt_dim1] = 1.;
 236                 } else if (icompq == 1) {
 237                     q[*n + (qstart - 1) * *n] = (d__[*n]>0) ? 1.0 : -1.0;
 238                     q[*n + (smlsiz + qstart - 1) * *n] = 1.;
 239                 }
 240                 d__[*n] = fabs(d__[*n]);
 241             }
 242             if (icompq == 2) {
 243                 F77_FUNC(slasd0,SLASD0)(&nsize, &sqre, &d__[start], &e[start],
 244                         &u[start + start * u_dim1], ldu,
 245                         &vt[start + start * vt_dim1],
 246                         ldvt, &smlsiz, &iwork[1], &work[wstart], info);
 247             } else {
 248                 F77_FUNC(slasda,SLASDA)(&icompq, &smlsiz, &nsize, &sqre, &d__[start],
 249                         &e[start], &q[start + (iu + qstart - 2) * *n], n,
 250                         &q[start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
 251                         &q[start + (difl + qstart - 2) * *n],
 252                         &q[start + (difr + qstart - 2) * *n],
 253                         &q[start + (z__ + qstart - 2) * *n],
 254                         &q[start + (poles + qstart - 2) * *n],
 255                         &iq[start + givptr * *n], &iq[start + givcol * *n], n,
 256                         &iq[start + perm * *n],
 257                         &q[start + (givnum + qstart - 2) * *n],
 258                         &q[start + (ic + qstart - 2) * *n],
 259                         &q[start + (is + qstart - 2) * *n], &work[wstart],
 260                         &iwork[1], info);
 261                 if (*info != 0) {
 262                     return;
 263                 }
 264             }
 265             start = i__ + 1;
 266         }
 267     }
 268     F77_FUNC(slascl,SLASCL)("G", &c_0, &c_0, &one, &orgnrm, n, &c_1, &d__[1], n, &ierr);
 269 L40:
 270     i__1 = *n;
 271     for (ii = 2; ii <= i__1; ++ii) {
 272         i__ = ii - 1;
 273         kk = i__;
 274         p = d__[i__];
 275         i__2 = *n;
 276         for (j = ii; j <= i__2; ++j) {
 277             if (d__[j] > p) {
 278                 kk = j;
 279                 p = d__[j];
 280             }
 281         }
 282         if (kk != i__) {
 283             d__[kk] = d__[i__];
 284             d__[i__] = p;
 285             if (icompq == 1) {
 286                 iq[i__] = kk;
 287             } else if (icompq == 2) {
 288                 F77_FUNC(sswap,SSWAP)(n, &u[i__ * u_dim1 + 1],&c_1,&u[kk*u_dim1+1],&c_1);
 289                 F77_FUNC(sswap,SSWAP)(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
 290             }
 291         } else if (icompq == 1) {
 292             iq[i__] = i__;
 293         }
 294     }
 295     if (icompq == 1) {
 296         if (iuplo == 1) {
 297             iq[*n] = 1;
 298         } else {
 299             iq[*n] = 0;
 300         }
 301     }
 302     if (iuplo == 2 && icompq == 2) {
 303         F77_FUNC(slasr,SLASR)("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu);
 304     }
 305
 306     return;
 307 }