src/gromacs/linearalgebra/gmx_lapack/dlasd6.cpp

   1 #include <cmath>
   2 #include "../gmx_blas.h"
   3 #include "../gmx_lapack.h"
   4
   5 void
   6 F77_FUNC(dlasd6,DLASD6)(int *icompq,
   7         int *nl,
   8         int *nr,
   9         int *sqre,
  10         double *d__,
  11         double *vf,
  12         double *vl,
  13         double *alpha,
  14         double *beta,
  15         int *idxq,
  16         int *perm,
  17         int *givptr,
  18         int *givcol,
  19         int *ldgcol,
  20         double *givnum,
  21         int *ldgnum,
  22         double *poles,
  23         double *difl,
  24         double *difr,
  25         double *z__,
  26         int *k,
  27         double *c__,
  28         double *s,
  29         double *work,
  30         int *iwork,
  31         int *info)
  32 {
  33     int givcol_dim1, givcol_offset, givnum_dim1, givnum_offset,
  34             poles_dim1, poles_offset, i__1;
  35     double d__1, d__2;
  36
  37     int i__, m, n, n1, n2, iw, idx, idxc, idxp, ivfw, ivlw;
  38     int isigma;
  39     double orgnrm;
  40     int c__0 = 0;
  41     double one = 1.0;
  42     int c__1 = 1;
  43     int c_n1 = -1;
  44
  45     --d__;
  46     --vf;
  47     --vl;
  48     --idxq;
  49     --perm;
  50     givcol_dim1 = *ldgcol;
  51     givcol_offset = 1 + givcol_dim1;
  52     givcol -= givcol_offset;
  53     poles_dim1 = *ldgnum;
  54     poles_offset = 1 + poles_dim1;
  55     poles -= poles_offset;
  56     givnum_dim1 = *ldgnum;
  57     givnum_offset = 1 + givnum_dim1;
  58     givnum -= givnum_offset;
  59     --difl;
  60     --difr;
  61     --z__;
  62     --work;
  63     --iwork;
  64
  65     *info = 0;
  66     n = *nl + *nr + 1;
  67     m = n + *sqre;
  68
  69     isigma = 1;
  70     iw = isigma + n;
  71     ivfw = iw + m;
  72     ivlw = ivfw + m;
  73
  74     idx = 1;
  75     idxc = idx + n;
  76     idxp = idxc + n;
  77
  78     d__1 = std::abs(*alpha);
  79     d__2 = std::abs(*beta);
  80     orgnrm = (d__1 > d__2) ? d__1 : d__2;
  81     d__[*nl + 1] = 0.;
  82     i__1 = n;
  83     for (i__ = 1; i__ <= i__1; ++i__) {
  84       d__1 = std::abs(d__[i__]);
  85         if (d__1 > orgnrm)
  86             orgnrm = d__1;
  87     }
  88     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &orgnrm, &one, &n, &c__1, &d__[1], &n, info);
  89     *alpha /= orgnrm;
  90     *beta /= orgnrm;
  91
  92     F77_FUNC(dlasd7,DLASD7)(icompq, nl, nr, sqre, k, &d__[1], &z__[1], &work[iw], &vf[1], &
  93             work[ivfw], &vl[1], &work[ivlw], alpha, beta, &work[isigma], &
  94             iwork[idx], &iwork[idxp], &idxq[1], &perm[1], givptr, &givcol[
  95             givcol_offset], ldgcol, &givnum[givnum_offset], ldgnum, c__, s,
  96             info);
  97
  98     F77_FUNC(dlasd8,DLASD8)(icompq, k, &d__[1], &z__[1], &vf[1], &vl[1], &difl[1], &difr[1],
  99             ldgnum, &work[isigma], &work[iw], info);
 100
 101     if (*icompq == 1) {
 102         F77_FUNC(dcopy,DCOPY)(k, &d__[1], &c__1, &poles[poles_dim1 + 1], &c__1);
 103         F77_FUNC(dcopy,DCOPY)(k, &work[isigma], &c__1, &poles[(poles_dim1 << 1) + 1], &c__1);
 104     }
 105
 106     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &one, &orgnrm, &n, &c__1, &d__[1], &n, info);
 107
 108     n1 = *k;
 109     n2 = n - *k;
 110     F77_FUNC(dlamrg,DLAMRG)(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);
 111
 112     return;
 113
 114 }
 115
 116