src/gromacs/linearalgebra/gmx_lapack/dbdsdc.c

   1 #include <ctype.h>
   2 #include <math.h>
   3 #include "../gmx_blas.h"
   4 #include "../gmx_lapack.h"
   5 #include "lapack_limits.h"
   6
   7 #include "types/simple.h"
   8
   9 void
  10 F77_FUNC(dbdsdc,DBDSDC)(const char *uplo,
  11         const char *compq,
  12         int *n,
  13         double *d__,
  14         double *e,
  15         double *u,
  16         int *ldu,
  17         double *vt,
  18         int *ldvt,
  19         double *q,
  20         int *iq,
  21         double *work,
  22         int *iwork,
  23         int *info)
  24 {
  25     int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  26     int i__, j, k;
  27     double p, r__;
  28     int z__, ic, ii, kk;
  29     double cs;
  30     int is, iu;
  31     double sn;
  32     int nm1;
  33     double eps;
  34     int ivt, difl, difr, ierr, perm, mlvl, sqre;
  35     int poles, iuplo, nsize, start;
  36     int givcol;
  37     int icompq;
  38     double orgnrm;
  39     int givnum, givptr, qstart, smlsiz, wstart, smlszp;
  40     double zero = 0.0;
  41     double one = 1.0;
  42     int c_0 = 0;
  43     int c_1 = 1;
  44
  45     --d__;
  46     --e;
  47     u_dim1 = *ldu;
  48     u_offset = 1 + u_dim1;
  49     u -= u_offset;
  50     vt_dim1 = *ldvt;
  51     vt_offset = 1 + vt_dim1;
  52     vt -= vt_offset;
  53     --q;
  54     --iq;
  55     --work;
  56     --iwork;
  57
  58     k = iu = z__ = ic = is = ivt = difl = difr = perm = 0;
  59     poles = givnum = givptr = givcol = 0;
  60
  61     smlsiz = DBDSDC_SMALLSIZE;
  62     *info = 0;
  63
  64     iuplo = (*uplo=='U' || *uplo=='u') ? 1 : 2;
  65
  66     switch(*compq) {
  67     case 'n':
  68     case 'N':
  69       icompq = 0;
  70       break;
  71     case 'p':
  72     case 'P':
  73       icompq = 1;
  74       break;
  75     case 'i':
  76     case 'I':
  77       icompq = 2;
  78       break;
  79     default:
  80       return;
  81     }
  82
  83     if (*n <= 0)
  84         return;
  85
  86     if (*n == 1) {
  87         if (icompq == 1) {
  88           q[1] = (d__[1]>0) ? 1.0 : -1.0;
  89           q[smlsiz * *n + 1] = 1.;
  90         } else if (icompq == 2) {
  91           u[u_dim1 + 1] = (d__[1]>0) ? 1.0 : -1.0;
  92           vt[vt_dim1 + 1] = 1.;
  93         }
  94         d__[1] = fabs(d__[1]);
  95         return;
  96     }
  97     nm1 = *n - 1;
  98     wstart = 1;
  99     qstart = 3;
 100     if (icompq == 1) {
 101         F77_FUNC(dcopy,DCOPY)(n, &d__[1], &c_1, &q[1], &c_1);
 102         i__1 = *n - 1;
 103         F77_FUNC(dcopy,DCOPY)(&i__1, &e[1], &c_1, &q[*n + 1], &c_1);
 104     }
 105     if (iuplo == 2) {
 106         qstart = 5;
 107         wstart = (*n << 1) - 1;
 108         i__1 = *n - 1;
 109         for (i__ = 1; i__ <= i__1; ++i__) {
 110             F77_FUNC(dlartg,DLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 111             d__[i__] = r__;
 112             e[i__] = sn * d__[i__ + 1];
 113             d__[i__ + 1] = cs * d__[i__ + 1];
 114             if (icompq == 1) {
 115                 q[i__ + (*n << 1)] = cs;
 116                 q[i__ + *n * 3] = sn;
 117             } else if (icompq == 2) {
 118                 work[i__] = cs;
 119                 work[nm1 + i__] = -sn;
 120             }
 121         }
 122     }
 123     if (icompq == 0) {
 124       F77_FUNC(dlasdq,DLASDQ)("U",&c_0,n,&c_0,&c_0,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
 125               &u[u_offset], ldu, &u[u_offset], ldu, &work[wstart], info);
 126         goto L40;
 127     }
 128     if (*n <= smlsiz) {
 129         if (icompq == 2) {
 130             F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
 131             F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
 132             F77_FUNC(dlasdq,DLASDQ)("U",&c_0,n,n,n,&c_0,&d__[1],&e[1],&vt[vt_offset],ldvt,
 133                     &u[u_offset],ldu,&u[u_offset],ldu,&work[wstart],info);
 134         } else if (icompq == 1) {
 135             iu = 1;
 136             ivt = iu + *n;
 137             F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &q[iu + (qstart - 1) * *n], n);
 138             F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &q[ivt + (qstart - 1) * *n], n);
 139             F77_FUNC(dlasdq,DLASDQ)("U", &c_0, n, n, n, &c_0, &d__[1], &e[1],
 140                     &q[ivt + (qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n],
 141                     n, &q[iu + (qstart - 1) * *n], n, &work[wstart], info);
 142         }
 143         goto L40;
 144     }
 145
 146     if (icompq == 2) {
 147         F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &u[u_offset], ldu);
 148         F77_FUNC(dlaset,DLASET)("A", n, n, &zero, &one, &vt[vt_offset], ldvt);
 149     }
 150
 151     orgnrm = F77_FUNC(dlanst,DLANST)("M", n, &d__[1], &e[1]);
 152     if ( fabs(orgnrm)<GMX_DOUBLE_MIN) {
 153         return;
 154     }
 155     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &orgnrm, &one, n, &c_1, &d__[1], n, &ierr);
 156     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &orgnrm, &one, &nm1, &c_1, &e[1], &nm1, &ierr);
 157
 158     eps = GMX_DOUBLE_EPS;
 159
 160     mlvl = (int) (log((double) (*n) / (double) (smlsiz + 1)) /
 161             log(2.)) + 1;
 162     smlszp = smlsiz + 1;
 163
 164     if (icompq == 1) {
 165         iu = 1;
 166         ivt = smlsiz + 1;
 167         difl = ivt + smlszp;
 168         difr = difl + mlvl;
 169         z__ = difr + (mlvl << 1);
 170         ic = z__ + mlvl;
 171         is = ic + 1;
 172         poles = is + 1;
 173         givnum = poles + (mlvl << 1);
 174
 175         k = 1;
 176         givptr = 2;
 177         perm = 3;
 178         givcol = perm + mlvl;
 179     }
 180
 181     i__1 = *n;
 182     for (i__ = 1; i__ <= i__1; ++i__) {
 183         if (fabs(d__[i__]) < eps)
 184             d__[i__] = (d__[i__]>0) ? eps : -eps;
 185     }
 186
 187     start = 1;
 188     sqre = 0;
 189
 190     i__1 = nm1;
 191     for (i__ = 1; i__ <= i__1; ++i__) {
 192         if (fabs(e[i__]) < eps || i__ == nm1) {
 193             if (i__ < nm1) {
 194                 nsize = i__ - start + 1;
 195             } else if (fabs(e[i__]) >= eps) {
 196                 nsize = *n - start + 1;
 197             } else {
 198                 nsize = i__ - start + 1;
 199                 if (icompq == 2) {
 200                     u[*n + *n * u_dim1] = (d__[*n]>0) ? 1.0 : -1.0;
 201                     vt[*n + *n * vt_dim1] = 1.;
 202                 } else if (icompq == 1) {
 203                     q[*n + (qstart - 1) * *n] = (d__[*n]>0) ? 1.0 : -1.0;
 204                     q[*n + (smlsiz + qstart - 1) * *n] = 1.;
 205                 }
 206                 d__[*n] = fabs(d__[*n]);
 207             }
 208             if (icompq == 2) {
 209                 F77_FUNC(dlasd0,DLASD0)(&nsize, &sqre, &d__[start], &e[start],
 210                         &u[start + start * u_dim1], ldu,
 211                         &vt[start + start * vt_dim1],
 212                         ldvt, &smlsiz, &iwork[1], &work[wstart], info);
 213             } else {
 214                 F77_FUNC(dlasda,DLASDA)(&icompq, &smlsiz, &nsize, &sqre, &d__[start],
 215                         &e[start], &q[start + (iu + qstart - 2) * *n], n,
 216                         &q[start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
 217                         &q[start + (difl + qstart - 2) * *n],
 218                         &q[start + (difr + qstart - 2) * *n],
 219                         &q[start + (z__ + qstart - 2) * *n],
 220                         &q[start + (poles + qstart - 2) * *n],
 221                         &iq[start + givptr * *n], &iq[start + givcol * *n], n,
 222                         &iq[start + perm * *n],
 223                         &q[start + (givnum + qstart - 2) * *n],
 224                         &q[start + (ic + qstart - 2) * *n],
 225                         &q[start + (is + qstart - 2) * *n], &work[wstart],
 226                         &iwork[1], info);
 227                 if (*info != 0) {
 228                     return;
 229                 }
 230             }
 231             start = i__ + 1;
 232         }
 233     }
 234     F77_FUNC(dlascl,DLASCL)("G", &c_0, &c_0, &one, &orgnrm, n, &c_1, &d__[1], n, &ierr);
 235 L40:
 236     i__1 = *n;
 237     for (ii = 2; ii <= i__1; ++ii) {
 238         i__ = ii - 1;
 239         kk = i__;
 240         p = d__[i__];
 241         i__2 = *n;
 242         for (j = ii; j <= i__2; ++j) {
 243             if (d__[j] > p) {
 244                 kk = j;
 245                 p = d__[j];
 246             }
 247         }
 248         if (kk != i__) {
 249             d__[kk] = d__[i__];
 250             d__[i__] = p;
 251             if (icompq == 1) {
 252                 iq[i__] = kk;
 253             } else if (icompq == 2) {
 254                 F77_FUNC(dswap,DSWAP)(n, &u[i__ * u_dim1 + 1],&c_1,&u[kk*u_dim1+1],&c_1);
 255                 F77_FUNC(dswap,DSWAP)(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
 256             }
 257         } else if (icompq == 1) {
 258             iq[i__] = i__;
 259         }
 260     }
 261     if (icompq == 1) {
 262         if (iuplo == 1) {
 263             iq[*n] = 1;
 264         } else {
 265             iq[*n] = 0;
 266         }
 267     }
 268     if (iuplo == 2 && icompq == 2) {
 269         F77_FUNC(dlasr,DLASR)("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu);
 270     }
 271
 272     return;
 273 }