src/gmxlib/gmx_lapack/dsytd2.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 <ctype.h>
  36 #include <math.h>
  37
  38 #include <types/simple.h>
  39
  40 #include "gmx_blas.h"
  41 #include "gmx_lapack.h"
  42
  43 void
  44 F77_FUNC(dsytd2,DSYTD2)(const char *    uplo,
  45         int *     n,
  46         double *  a,
  47         int *     lda,
  48         double *  d,
  49         double *  e,
  50         double *  tau,
  51         int *     info)
  52 {
  53   double minusone,zero;
  54   double taui,alpha,tmp;
  55   int ti1,ti2,ti3,i;
  56   const char ch=toupper(*uplo);
  57
  58   zero = 0.0;
  59   minusone = -1.0;
  60
  61   if(*n<=0)
  62     return;
  63
  64   if(ch=='U') {
  65     for(i=*n-1;i>=1;i--) {
  66
  67       ti1 = 1;
  68       F77_FUNC(dlarfg,DLARFG)(&i,&(a[i*(*lda)+(i-1)]),&(a[i*(*lda)+0]),&ti1,&taui);
  69       e[i-1] = a[i*(*lda) + (i-1)];
  70       if(fabs(taui)>GMX_DOUBLE_MIN) {
  71         a[i*(*lda)+(i-1)] = 1.0;
  72
  73         ti1 = 1;
  74         F77_FUNC(dsymv,DSYMV)("U",&i,&taui,a,lda,&(a[i*(*lda)+0]),&ti1,&zero,tau,&ti1);
  75
  76         tmp = F77_FUNC(ddot,DDOT)(&i,tau,&ti1,&(a[i*(*lda)+0]),&ti1);
  77
  78         alpha = -0.5*taui*tmp;
  79
  80         F77_FUNC(daxpy,DAXPY)(&i,&alpha,&(a[i*(*lda)+0]),&ti1,tau,&ti1);
  81
  82         F77_FUNC(dsyr2,DSYR2)("U",&i,&minusone,&(a[i*(*lda)+0]),&ti1,tau,&ti1,a,lda);
  83
  84         a[i*(*lda)+(i-1)] = e[i-1];
  85
  86       }
  87       d[i] = a[i*(*lda)+i];
  88       tau[i-1] = taui;
  89     }
  90     d[0] = a[0];
  91
  92   } else {
  93     /* lower */
  94
  95     for(i=1;i<*n;i++) {
  96
  97       ti1 = *n - i;
  98       ti2 = ( *n < i+2) ? *n : i+2;
  99       ti3 = 1;
 100       F77_FUNC(dlarfg,DLARFG)(&ti1,&(a[(i-1)*(*lda)+(i)]),&(a[(i-1)*(*lda)+ti2-1]),&ti3,&taui);
 101
 102       e[i-1] = a[(i-1)*(*lda) + (i)];
 103
 104       if(fabs(taui)>GMX_DOUBLE_MIN) {
 105         a[(i-1)*(*lda)+(i)] = 1.0;
 106
 107         ti1 = *n - i;
 108         ti2 = 1;
 109         F77_FUNC(dsymv,DSYMV)(uplo,&ti1,&taui,&(a[i*(*lda)+i]),lda,&(a[(i-1)*(*lda)+i]),
 110                &ti2,&zero,&(tau[i-1]),&ti2);
 111
 112         tmp = F77_FUNC(ddot,DDOT)(&ti1,&(tau[i-1]),&ti2,&(a[(i-1)*(*lda)+i]),&ti2);
 113
 114         alpha = -0.5*taui*tmp;
 115
 116         F77_FUNC(daxpy,DAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+i]),&ti2,&(tau[i-1]),&ti2);
 117
 118         F77_FUNC(dsyr2,DSYR2)(uplo,&ti1,&minusone,&(a[(i-1)*(*lda)+i]),&ti2,&(tau[i-1]),&ti2,
 119                &(a[(i)*(*lda)+i]),lda);
 120
 121         a[(i-1)*(*lda)+(i)] = e[i-1];
 122
 123       }
 124       d[i-1] = a[(i-1)*(*lda)+i-1];
 125       tau[i-1] = taui;
 126     }
 127     d[*n-1] = a[(*n-1)*(*lda)+(*n-1)];
 128
 129   }
 130   return;
 131 }