src/gmxlib/gmx_lapack/dlasda.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(dlasda,DLASDA)(int *icompq,
   6         int *smlsiz,
   7         int *n,
   8         int *sqre,
   9         double *d__,
  10         double *e,
  11         double *u,
  12         int *ldu,
  13         double *vt,
  14         int *k,
  15         double *difl,
  16         double *difr,
  17         double *z__,
  18         double *poles,
  19         int *givptr,
  20         int *givcol,
  21         int *ldgcol,
  22         int *perm,
  23         double *givnum,
  24         double *c__,
  25         double *s,
  26         double *work,
  27         int *iwork,
  28         int *info)
  29 {
  30     int givcol_dim1, givcol_offset, perm_dim1, perm_offset, difl_dim1,
  31             difl_offset, difr_dim1, difr_offset, givnum_dim1, givnum_offset,
  32             poles_dim1, poles_offset, u_dim1, u_offset, vt_dim1, vt_offset,
  33             z_dim1, z_offset, i__1, i__2;
  34
  35     int i__, j, m, i1, ic, lf, nd, ll, nl, vf, nr, vl, im1, ncc,
  36             nlf, nrf, vfi, iwk, vli, lvl, nru, ndb1, nlp1, lvl2, nrp1;
  37     double beta;
  38     int idxq, nlvl;
  39     double alpha;
  40     int inode, ndiml, ndimr, idxqi, itemp;
  41     int sqrei;
  42     int nwork1, nwork2;
  43     int smlszp;
  44     int c__0 = 0;
  45     double zero = 0.0;
  46     double one = 1.;
  47     int c__1 = 1;
  48     --d__;
  49     --e;
  50     givnum_dim1 = *ldu;
  51     givnum_offset = 1 + givnum_dim1;
  52     givnum -= givnum_offset;
  53     poles_dim1 = *ldu;
  54     poles_offset = 1 + poles_dim1;
  55     poles -= poles_offset;
  56     z_dim1 = *ldu;
  57     z_offset = 1 + z_dim1;
  58     z__ -= z_offset;
  59     difr_dim1 = *ldu;
  60     difr_offset = 1 + difr_dim1;
  61     difr -= difr_offset;
  62     difl_dim1 = *ldu;
  63     difl_offset = 1 + difl_dim1;
  64     difl -= difl_offset;
  65     vt_dim1 = *ldu;
  66     vt_offset = 1 + vt_dim1;
  67     vt -= vt_offset;
  68     u_dim1 = *ldu;
  69     u_offset = 1 + u_dim1;
  70     u -= u_offset;
  71     --k;
  72     --givptr;
  73     perm_dim1 = *ldgcol;
  74     perm_offset = 1 + perm_dim1;
  75     perm -= perm_offset;
  76     givcol_dim1 = *ldgcol;
  77     givcol_offset = 1 + givcol_dim1;
  78     givcol -= givcol_offset;
  79     --c__;
  80     --s;
  81     --work;
  82     --iwork;
  83     *info = 0;
  84
  85     m = *n + *sqre;
  86
  87     if (*n <= *smlsiz) {
  88         if (*icompq == 0) {
  89             F77_FUNC(dlasdq,DLASDQ)("U", sqre, n, &c__0, &c__0, &c__0, &d__[1], &e[1], &vt[
  90                     vt_offset], ldu, &u[u_offset], ldu, &u[u_offset], ldu, &
  91                     work[1], info);
  92         } else {
  93             F77_FUNC(dlasdq,DLASDQ)("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset]
  94                     , ldu, &u[u_offset], ldu, &u[u_offset], ldu, &work[1],
  95                     info);
  96         }
  97         return;
  98     }
  99
 100     inode = 1;
 101     ndiml = inode + *n;
 102     ndimr = ndiml + *n;
 103     idxq = ndimr + *n;
 104     iwk = idxq + *n;
 105
 106     ncc = 0;
 107     nru = 0;
 108
 109     smlszp = *smlsiz + 1;
 110     vf = 1;
 111     vl = vf + m;
 112     nwork1 = vl + m;
 113     nwork2 = nwork1 + smlszp * smlszp;
 114
 115     F77_FUNC(dlasdt,DLASDT)(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
 116             smlsiz);
 117
 118     ndb1 = (nd + 1) / 2;
 119     i__1 = nd;
 120     for (i__ = ndb1; i__ <= i__1; ++i__) {
 121         i1 = i__ - 1;
 122         ic = iwork[inode + i1];
 123         nl = iwork[ndiml + i1];
 124         nlp1 = nl + 1;
 125         nr = iwork[ndimr + i1];
 126         nlf = ic - nl;
 127         nrf = ic + 1;
 128         idxqi = idxq + nlf - 2;
 129         vfi = vf + nlf - 1;
 130         vli = vl + nlf - 1;
 131         sqrei = 1;
 132         if (*icompq == 0) {
 133             F77_FUNC(dlaset,DLASET)("A", &nlp1, &nlp1, &zero, &one, &work[nwork1], &smlszp);
 134             F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nl, &nlp1, &nru, &ncc, &d__[nlf], &e[nlf], &
 135                     work[nwork1], &smlszp, &work[nwork2], &nl, &work[nwork2],
 136                     &nl, &work[nwork2], info);
 137             itemp = nwork1 + nl * smlszp;
 138             F77_FUNC(dcopy,DCOPY)(&nlp1, &work[nwork1], &c__1, &work[vfi], &c__1);
 139             F77_FUNC(dcopy,DCOPY)(&nlp1, &work[itemp], &c__1, &work[vli], &c__1);
 140         } else {
 141             F77_FUNC(dlaset,DLASET)("A", &nl, &nl, &zero, &one, &u[nlf + u_dim1], ldu);
 142             F77_FUNC(dlaset,DLASET)("A", &nlp1, &nlp1, &zero, &one, &vt[nlf + vt_dim1],
 143                     ldu);
 144             F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &
 145                     vt[nlf + vt_dim1], ldu, &u[nlf + u_dim1], ldu, &u[nlf +
 146                     u_dim1], ldu, &work[nwork1], info);
 147             F77_FUNC(dcopy,DCOPY)(&nlp1, &vt[nlf + vt_dim1], &c__1, &work[vfi], &c__1);
 148             F77_FUNC(dcopy,DCOPY)(&nlp1, &vt[nlf + nlp1 * vt_dim1], &c__1, &work[vli], &c__1)
 149                     ;
 150         }
 151         if (*info != 0) {
 152             return;
 153         }
 154         i__2 = nl;
 155         for (j = 1; j <= i__2; ++j) {
 156             iwork[idxqi + j] = j;
 157         }
 158         if (i__ == nd && *sqre == 0) {
 159             sqrei = 0;
 160         } else {
 161             sqrei = 1;
 162         }
 163         idxqi += nlp1;
 164         vfi += nlp1;
 165         vli += nlp1;
 166         nrp1 = nr + sqrei;
 167         if (*icompq == 0) {
 168             F77_FUNC(dlaset,DLASET)("A", &nrp1, &nrp1, &zero, &one, &work[nwork1], &smlszp);
 169             F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nr, &nrp1, &nru, &ncc, &d__[nrf], &e[nrf], &
 170                     work[nwork1], &smlszp, &work[nwork2], &nr, &work[nwork2],
 171                     &nr, &work[nwork2], info);
 172             itemp = nwork1 + (nrp1 - 1) * smlszp;
 173             F77_FUNC(dcopy,DCOPY)(&nrp1, &work[nwork1], &c__1, &work[vfi], &c__1);
 174             F77_FUNC(dcopy,DCOPY)(&nrp1, &work[itemp], &c__1, &work[vli], &c__1);
 175         } else {
 176             F77_FUNC(dlaset,DLASET)("A", &nr, &nr, &zero, &one, &u[nrf + u_dim1], ldu);
 177             F77_FUNC(dlaset,DLASET)("A", &nrp1, &nrp1, &zero, &one, &vt[nrf + vt_dim1],
 178                     ldu);
 179             F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &
 180                     vt[nrf + vt_dim1], ldu, &u[nrf + u_dim1], ldu, &u[nrf +
 181                     u_dim1], ldu, &work[nwork1], info);
 182             F77_FUNC(dcopy,DCOPY)(&nrp1, &vt[nrf + vt_dim1], &c__1, &work[vfi], &c__1);
 183             F77_FUNC(dcopy,DCOPY)(&nrp1, &vt[nrf + nrp1 * vt_dim1], &c__1, &work[vli], &c__1)
 184                     ;
 185         }
 186         if (*info != 0) {
 187             return;
 188         }
 189         i__2 = nr;
 190         for (j = 1; j <= i__2; ++j) {
 191             iwork[idxqi + j] = j;
 192         }
 193     }
 194
 195     j = (1 << nlvl);
 196
 197     for (lvl = nlvl; lvl >= 1; --lvl) {
 198         lvl2 = (lvl << 1) - 1;
 199
 200         if (lvl == 1) {
 201             lf = 1;
 202             ll = 1;
 203         } else {
 204             i__1 = lvl - 1;
 205             lf = (1 << (lvl-1));
 206             ll = (lf << 1) - 1;
 207         }
 208         i__1 = ll;
 209         for (i__ = lf; i__ <= i__1; ++i__) {
 210             im1 = i__ - 1;
 211             ic = iwork[inode + im1];
 212             nl = iwork[ndiml + im1];
 213             nr = iwork[ndimr + im1];
 214             nlf = ic - nl;
 215             nrf = ic + 1;
 216             if (i__ == ll) {
 217                 sqrei = *sqre;
 218             } else {
 219                 sqrei = 1;
 220             }
 221             vfi = vf + nlf - 1;
 222             vli = vl + nlf - 1;
 223             idxqi = idxq + nlf - 1;
 224             alpha = d__[ic];
 225             beta = e[ic];
 226             if (*icompq == 0) {
 227                 F77_FUNC(dlasd6,DLASD6)(icompq, &nl, &nr, &sqrei, &d__[nlf], &work[vfi], &
 228                         work[vli], &alpha, &beta, &iwork[idxqi], &perm[
 229                         perm_offset], &givptr[1], &givcol[givcol_offset],
 230                         ldgcol, &givnum[givnum_offset], ldu, &poles[
 231                         poles_offset], &difl[difl_offset], &difr[difr_offset],
 232                          &z__[z_offset], &k[1], &c__[1], &s[1], &work[nwork1],
 233                          &iwork[iwk], info);
 234             } else {
 235                 --j;
 236                 F77_FUNC(dlasd6,DLASD6)(icompq, &nl, &nr, &sqrei, &d__[nlf], &work[vfi], &
 237                         work[vli], &alpha, &beta, &iwork[idxqi], &perm[nlf +
 238                         lvl * perm_dim1], &givptr[j], &givcol[nlf + lvl2 *
 239                         givcol_dim1], ldgcol, &givnum[nlf + lvl2 *
 240                         givnum_dim1], ldu, &poles[nlf + lvl2 * poles_dim1], &
 241                         difl[nlf + lvl * difl_dim1], &difr[nlf + lvl2 *
 242                         difr_dim1], &z__[nlf + lvl * z_dim1], &k[j], &c__[j],
 243                         &s[j], &work[nwork1], &iwork[iwk], info);
 244             }
 245             if (*info != 0) {
 246                 return;
 247             }
 248         }
 249     }
 250
 251     return;
 252
 253 }
 254
 255