src/gmxlib/gmx_lapack/dlasd0.c

   1 #include "gmx_lapack.h"
   2
   3 void
   4 F77_FUNC(dlasd0,DLASD0)(int *n,
   5         int *sqre,
   6         double *d__,
   7         double *e,
   8         double *u,
   9         int *ldu,
  10         double *vt,
  11         int *ldvt,
  12         int *smlsiz,
  13         int *iwork,
  14         double *work,
  15         int *info)
  16 {
  17     int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  18
  19     int i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf,
  20             iwk, lvl, ndb1, nlp1, nrp1;
  21     double beta;
  22     int idxq, nlvl;
  23     double alpha;
  24     int inode, ndiml, idxqc, ndimr, itemp, sqrei;
  25     int c__0 = 0;
  26
  27
  28     --d__;
  29     --e;
  30     u_dim1 = *ldu;
  31     u_offset = 1 + u_dim1;
  32     u -= u_offset;
  33     vt_dim1 = *ldvt;
  34     vt_offset = 1 + vt_dim1;
  35     vt -= vt_offset;
  36     --iwork;
  37     --work;
  38
  39     *info = 0;
  40
  41     if (*n < 0) {
  42         *info = -1;
  43     } else if (*sqre < 0 || *sqre > 1) {
  44         *info = -2;
  45     }
  46
  47     m = *n + *sqre;
  48
  49     if (*ldu < *n) {
  50         *info = -6;
  51     } else if (*ldvt < m) {
  52         *info = -8;
  53     } else if (*smlsiz < 3) {
  54         *info = -9;
  55     }
  56     if (*info != 0) {
  57         i__1 = -(*info);
  58         return;
  59     }
  60
  61     if (*n <= *smlsiz) {
  62         F77_FUNC(dlasdq,DLASDQ)("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset],
  63                 ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[1], info);
  64         return;
  65     }
  66
  67     inode = 1;
  68     ndiml = inode + *n;
  69     ndimr = ndiml + *n;
  70     idxq = ndimr + *n;
  71     iwk = idxq + *n;
  72     F77_FUNC(dlasdt,DLASDT)(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
  73             smlsiz);
  74
  75     ndb1 = (nd + 1) / 2;
  76     ncc = 0;
  77     i__1 = nd;
  78     for (i__ = ndb1; i__ <= i__1; ++i__) {
  79
  80         i1 = i__ - 1;
  81         ic = iwork[inode + i1];
  82         nl = iwork[ndiml + i1];
  83         nlp1 = nl + 1;
  84         nr = iwork[ndimr + i1];
  85         nrp1 = nr + 1;
  86         nlf = ic - nl;
  87         nrf = ic + 1;
  88         sqrei = 1;
  89         F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &vt[
  90                 nlf + nlf * vt_dim1], ldvt, &u[nlf + nlf * u_dim1], ldu, &u[
  91                 nlf + nlf * u_dim1], ldu, &work[1], info);
  92         if (*info != 0) {
  93             return;
  94         }
  95         itemp = idxq + nlf - 2;
  96         i__2 = nl;
  97         for (j = 1; j <= i__2; ++j) {
  98             iwork[itemp + j] = j;
  99         }
 100         if (i__ == nd) {
 101             sqrei = *sqre;
 102         } else {
 103             sqrei = 1;
 104         }
 105         nrp1 = nr + sqrei;
 106         F77_FUNC(dlasdq,DLASDQ)("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &vt[
 107                 nrf + nrf * vt_dim1], ldvt, &u[nrf + nrf * u_dim1], ldu, &u[
 108                 nrf + nrf * u_dim1], ldu, &work[1], info);
 109         if (*info != 0) {
 110             return;
 111         }
 112         itemp = idxq + ic;
 113         i__2 = nr;
 114         for (j = 1; j <= i__2; ++j) {
 115             iwork[itemp + j - 1] = j;
 116         }
 117     }
 118
 119     for (lvl = nlvl; lvl >= 1; --lvl) {
 120
 121         if (lvl == 1) {
 122             lf = 1;
 123             ll = 1;
 124         } else {
 125             i__1 = lvl - 1;
 126             lf = (1 << i__1);
 127             ll = (lf << 1) - 1;
 128         }
 129         i__1 = ll;
 130         for (i__ = lf; i__ <= i__1; ++i__) {
 131             im1 = i__ - 1;
 132             ic = iwork[inode + im1];
 133             nl = iwork[ndiml + im1];
 134             nr = iwork[ndimr + im1];
 135             nlf = ic - nl;
 136             if (*sqre == 0 && i__ == ll) {
 137                 sqrei = *sqre;
 138             } else {
 139                 sqrei = 1;
 140             }
 141             idxqc = idxq + nlf - 1;
 142             alpha = d__[ic];
 143             beta = e[ic];
 144             F77_FUNC(dlasd1,DLASD1)(&nl, &nr, &sqrei, &d__[nlf], &alpha, &beta, &u[nlf + nlf *
 145                      u_dim1], ldu, &vt[nlf + nlf * vt_dim1], ldvt, &iwork[
 146                     idxqc], &iwork[iwk], &work[1], info);
 147             if (*info != 0) {
 148                 return;
 149             }
 150         }
 151     }
 152
 153     return;
 154
 155 }