src/gmxlib/gmx_lapack/dlasd1.c

   1 #include <math.h>
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(dlasd1,DLASD1)(int *nl,
   6         int *nr,
   7         int *sqre,
   8         double *d__,
   9         double *alpha,
  10         double *beta,
  11         double *u,
  12         int *ldu,
  13         double *vt,
  14         int *ldvt,
  15         int *idxq,
  16         int *iwork,
  17         double *work,
  18         int *info)
  19 {
  20     int u_dim1, u_offset, vt_dim1, vt_offset, i__1;
  21     double d__1, d__2;
  22
  23     int i__, k, m, n, n1, n2, iq, iz, iu2, ldq, idx, ldu2, ivt2,
  24             idxc, idxp, ldvt2;
  25     int isigma;
  26     double orgnrm;
  27     int coltyp;
  28     int c__0 = 0;
  29     double one = 1.0;
  30     int c__1 = 1;
  31     int c_n1 = -1;
  32
  33     --d__;
  34     u_dim1 = *ldu;
  35     u_offset = 1 + u_dim1;
  36     u -= u_offset;
  37     vt_dim1 = *ldvt;
  38     vt_offset = 1 + vt_dim1;
  39     vt -= vt_offset;
  40     --idxq;
  41     --iwork;
  42     --work;
  43
  44     *info = 0;
  45
  46     if (*nl < 1) {
  47         *info = -1;
  48     } else if (*nr < 1) {
  49         *info = -2;
  50     } else if (*sqre < 0 || *sqre > 1) {
  51         *info = -3;
  52     }
  53     if (*info != 0) {
  54         i__1 = -(*info);
  55         return;
  56     }
  57
  58     n = *nl + *nr + 1;
  59     m = n + *sqre;
  60
  61
  62     ldu2 = n;
  63     ldvt2 = m;
  64
  65     iz = 1;
  66     isigma = iz + m;
  67     iu2 = isigma + n;
  68     ivt2 = iu2 + ldu2 * n;
  69     iq = ivt2 + ldvt2 * m;
  70
  71     idx = 1;
  72     idxc = idx + n;
  73     coltyp = idxc + n;
  74     idxp = coltyp + n;
  75
  76     d__1 = fabs(*alpha);
  77     d__2 = fabs(*beta);
  78     orgnrm = (d__1>d__2) ? d__1 : d__2;
  79     d__[*nl + 1] = 0.;
  80     i__1 = n;
  81     for (i__ = 1; i__ <= i__1; ++i__) {
  82         if (fabs(d__[i__]) > orgnrm) {
  83             orgnrm = fabs(d__[i__]);
  84         }
  85     }
  86     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &orgnrm, &one, &n, &c__1, &d__[1], &n, info);
  87     *alpha /= orgnrm;
  88     *beta /= orgnrm;
  89
  90     F77_FUNC(dlasd2,DLASD2)(nl, nr, sqre, &k, &d__[1], &work[iz], alpha, beta, &u[u_offset],
  91             ldu, &vt[vt_offset], ldvt, &work[isigma], &work[iu2], &ldu2, &
  92             work[ivt2], &ldvt2, &iwork[idxp], &iwork[idx], &iwork[idxc], &
  93             idxq[1], &iwork[coltyp], info);
  94
  95     ldq = k;
  96     F77_FUNC(dlasd3,DLASD3)(nl, nr, sqre, &k, &d__[1], &work[iq], &ldq, &work[isigma], &u[
  97             u_offset], ldu, &work[iu2], &ldu2, &vt[vt_offset], ldvt, &work[
  98             ivt2], &ldvt2, &iwork[idxc], &iwork[coltyp], &work[iz], info);
  99     if (*info != 0) {
 100         return;
 101     }
 102     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &one, &orgnrm, &n, &c__1, &d__[1], &n, info);
 103
 104     n1 = k;
 105     n2 = n - k;
 106     F77_FUNC(dlamrg,DLAMRG)(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);
 107
 108     return;
 109
 110 }