src/gmxlib/gmx_lapack/dlasd6.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38
  39 void
  40 F77_FUNC(dlasd6,DLASD6)(int *icompq,
  41         int *nl,
  42         int *nr,
  43         int *sqre,
  44         double *d__,
  45         double *vf,
  46         double *vl,
  47         double *alpha,
  48         double *beta,
  49         int *idxq,
  50         int *perm,
  51         int *givptr,
  52         int *givcol,
  53         int *ldgcol,
  54         double *givnum,
  55         int *ldgnum,
  56         double *poles,
  57         double *difl,
  58         double *difr,
  59         double *z__,
  60         int *k,
  61         double *c__,
  62         double *s,
  63         double *work,
  64         int *iwork,
  65         int *info)
  66 {
  67     int givcol_dim1, givcol_offset, givnum_dim1, givnum_offset,
  68             poles_dim1, poles_offset, i__1;
  69     double d__1, d__2;
  70
  71     int i__, m, n, n1, n2, iw, idx, idxc, idxp, ivfw, ivlw;
  72     int isigma;
  73     double orgnrm;
  74     int c__0 = 0;
  75     double one = 1.0;
  76     int c__1 = 1;
  77     int c_n1 = -1;
  78
  79     --d__;
  80     --vf;
  81     --vl;
  82     --idxq;
  83     --perm;
  84     givcol_dim1 = *ldgcol;
  85     givcol_offset = 1 + givcol_dim1;
  86     givcol -= givcol_offset;
  87     poles_dim1 = *ldgnum;
  88     poles_offset = 1 + poles_dim1;
  89     poles -= poles_offset;
  90     givnum_dim1 = *ldgnum;
  91     givnum_offset = 1 + givnum_dim1;
  92     givnum -= givnum_offset;
  93     --difl;
  94     --difr;
  95     --z__;
  96     --work;
  97     --iwork;
  98
  99     *info = 0;
 100     n = *nl + *nr + 1;
 101     m = n + *sqre;
 102
 103     isigma = 1;
 104     iw = isigma + n;
 105     ivfw = iw + m;
 106     ivlw = ivfw + m;
 107
 108     idx = 1;
 109     idxc = idx + n;
 110     idxp = idxc + n;
 111
 112     d__1 = fabs(*alpha);
 113     d__2 = fabs(*beta);
 114     orgnrm = (d__1 > d__2) ? d__1 : d__2;
 115     d__[*nl + 1] = 0.;
 116     i__1 = n;
 117     for (i__ = 1; i__ <= i__1; ++i__) {
 118       d__1 = fabs(d__[i__]);
 119         if (d__1 > orgnrm)
 120             orgnrm = d__1;
 121     }
 122     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &orgnrm, &one, &n, &c__1, &d__[1], &n, info);
 123     *alpha /= orgnrm;
 124     *beta /= orgnrm;
 125
 126     F77_FUNC(dlasd7,DLASD7)(icompq, nl, nr, sqre, k, &d__[1], &z__[1], &work[iw], &vf[1], &
 127             work[ivfw], &vl[1], &work[ivlw], alpha, beta, &work[isigma], &
 128             iwork[idx], &iwork[idxp], &idxq[1], &perm[1], givptr, &givcol[
 129             givcol_offset], ldgcol, &givnum[givnum_offset], ldgnum, c__, s,
 130             info);
 131
 132     F77_FUNC(dlasd8,DLASD8)(icompq, k, &d__[1], &z__[1], &vf[1], &vl[1], &difl[1], &difr[1],
 133             ldgnum, &work[isigma], &work[iw], info);
 134
 135     if (*icompq == 1) {
 136         F77_FUNC(dcopy,DCOPY)(k, &d__[1], &c__1, &poles[poles_dim1 + 1], &c__1);
 137         F77_FUNC(dcopy,DCOPY)(k, &work[isigma], &c__1, &poles[(poles_dim1 << 1) + 1], &c__1);
 138     }
 139
 140     F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &one, &orgnrm, &n, &c__1, &d__[1], &n, info);
 141
 142     n1 = *k;
 143     n2 = n - *k;
 144     F77_FUNC(dlamrg,DLAMRG)(&n1, &n2, &d__[1], &c__1, &c_n1, &idxq[1]);
 145
 146     return;
 147
 148 }
 149
 150