src/gmxlib/gmx_lapack/slatrd.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 "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38 #include "lapack_limits.h"
  39
  40
  41 void
  42 F77_FUNC(slatrd,SLATRD)(const char *  uplo,
  43        int  *   n,
  44        int  *   nb,
  45        float * a,
  46        int *    lda,
  47        float * e,
  48        float * tau,
  49        float * w,
  50        int *    ldw)
  51 {
  52   int i,iw;
  53   int ti1,ti2,ti3;
  54   float one,zero,minusone,alpha;
  55   const char ch=toupper(*uplo);
  56
  57   one=1.0;
  58   minusone=-1.0;
  59   zero=0.0;
  60
  61   if(*n<=0)
  62     return;
  63
  64   if(ch=='U') {
  65     for(i=*n;i>=(*n-*nb+1);i--) {
  66       iw = i -*n + *nb;
  67
  68       if(i<*n) {
  69         ti1 = *n-i;
  70         ti2 = 1;
  71         /* BLAS */
  72         F77_FUNC(sgemv,SGEMV)("N",&i,&ti1,&minusone, &(a[ i*(*lda) + 0]),lda,&(w[iw*(*ldw)+(i-1)]),
  73                ldw,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
  74         /* BLAS */
  75         F77_FUNC(sgemv,SGEMV)("N",&i,&ti1,&minusone, &(w[ iw*(*ldw) + 0]),ldw,&(a[i*(*lda)+(i-1)]),
  76                lda,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
  77       }
  78
  79       if(i>1) {
  80         /*  Generate elementary reflector H(i) to annihilate
  81          *              A(1:i-2,i)
  82          */
  83         ti1 = i-1;
  84         ti2 = 1;
  85
  86         /* LAPACK */
  87         F77_FUNC(slarfg,SLARFG)(&ti1,&(a[(i-1)*(*lda)+(i-2)]),&(a[(i-1)*(*lda)+0]),&ti2,&(tau[i-2]));
  88
  89         e[i-2] = a[(i-1)*(*lda)+(i-2)];
  90         a[(i-1)*(*lda)+(i-2)] = 1.0;
  91
  92         /* Compute W(1:i-1,i) */
  93         ti1 = i-1;
  94         ti2 = 1;
  95
  96         /* BLAS */
  97         F77_FUNC(ssymv,SSYMV)("U",&ti1,&one,a,lda,&(a[(i-1)*(*lda)+0]),&ti2,&zero,
  98                &(w[(iw-1)*(*ldw)+0]),&ti2);
  99         if(i<*n) {
 100           ti1 = i-1;
 101           ti2 = *n-i;
 102           ti3 = 1;
 103           /* BLAS */
 104           F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(w[iw*(*ldw)+0]),ldw,&(a[(i-1)*(*lda)+0]),&ti3,
 105                  &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
 106
 107           /* BLAS */
 108           F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(a[i*(*lda)+0]),lda,&(w[(iw-1)*(*ldw)+i]),&ti3,
 109                  &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
 110
 111           /* BLAS */
 112           F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(a[i*(*lda)+0]),lda,&(a[(i-1)*(*lda)+0]),&ti3,
 113                  &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
 114
 115           /* BLAS */
 116           F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(w[iw*(*ldw)+0]),ldw,&(w[(iw-1)*(*ldw)+i]),&ti3,
 117                  &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
 118         }
 119
 120         ti1 = i-1;
 121         ti2 = 1;
 122         /* BLAS */
 123         F77_FUNC(sscal,SSCAL)(&ti1,&(tau[i-2]),&(w[(iw-1)*(*ldw)+0]),&ti2);
 124
 125         alpha = -0.5*tau[i-2]*F77_FUNC(sdot,SDOT)(&ti1,&(w[(iw-1)*(*ldw)+0]),&ti2,
 126                                     &(a[(i-1)*(*lda)+0]),&ti2);
 127
 128         /* BLAS */
 129         F77_FUNC(saxpy,SAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+0]),&ti2,&(w[(iw-1)*(*ldw)+0]),&ti2);
 130
 131       }
 132     }
 133   } else {
 134     /* lower */
 135     for(i=1;i<=*nb;i++) {
 136
 137       ti1 = *n-i+1;
 138       ti2 = i-1;
 139       ti3 = 1;
 140       /* BLAS */
 141       F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone, &(a[ i-1 ]),lda,&(w[ i-1 ]),
 142                ldw,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
 143       /* BLAS */
 144       F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone, &(w[ i-1 ]),ldw,&(a[ i-1 ]),
 145                lda,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
 146
 147       if(i<*n) {
 148         ti1 = *n - i;
 149         ti2 = (*n < i+2 ) ? *n : (i+2);
 150         ti3 = 1;
 151         /* LAPACK */
 152         F77_FUNC(slarfg,SLARFG)(&ti1,&(a[(i-1)*(*lda)+(i)]),&(a[(i-1)*(*lda)+(ti2-1)]),&ti3,&(tau[i-1]));
 153         e[i-1] = a[(i-1)*(*lda)+(i)];
 154         a[(i-1)*(*lda)+(i)] = 1.0;
 155
 156         ti1 = *n - i;
 157         ti2 = 1;
 158         F77_FUNC(ssymv,SSYMV)("L",&ti1,&one,&(a[i*(*lda)+i]),lda,&(a[(i-1)*(*lda)+i]),&ti2,
 159                &zero,&(w[(i-1)*(*ldw)+i]),&ti2);
 160         ti1 = *n - i;
 161         ti2 = i-1;
 162         ti3 = 1;
 163         /* BLAS */
 164         F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(w[ i ]),ldw,&(a[(i-1)*(*lda)+i]),&ti3,
 165                &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
 166
 167         /* BLAS */
 168         F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(a[ i ]),lda,&(w[(i-1)*(*ldw)+0]),&ti3,
 169                &one,&(w[(i-1)*(*ldw)+i]),&ti3);
 170
 171         /* BLAS */
 172         F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(a[ i ]),lda,&(a[(i-1)*(*lda)+i]),&ti3,
 173                &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
 174
 175         /* BLAS */
 176         F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(w[ i ]),ldw,&(w[(i-1)*(*ldw)+0]),&ti3,
 177                &one,&(w[(i-1)*(*ldw)+i]),&ti3);
 178
 179         F77_FUNC(sscal,SSCAL)(&ti1,&(tau[i-1]),&(w[(i-1)*(*ldw)+i]),&ti3);
 180         alpha = -0.5*tau[i-1]*F77_FUNC(sdot,SDOT)(&ti1,&(w[(i-1)*(*ldw)+i]),&ti3,
 181                                    &(a[(i-1)*(*lda)+i]),&ti3);
 182
 183         F77_FUNC(saxpy,SAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+i]),&ti3,&(w[(i-1)*(*ldw)+i]),&ti3);
 184       }
 185     }
 186   }
 187   return;
 188 }
 189
 190
 191
 192