src/gmxlib/gmx_lapack/dlatrd.c

   1 #include <ctype.h>
   2 #include "gmx_blas.h"
   3 #include "gmx_lapack.h"
   4 #include "lapack_limits.h"
   5
   6
   7 void
   8 F77_FUNC(dlatrd,DLATRD)(const char *  uplo,
   9        int  *   n,
  10        int  *   nb,
  11        double * a,
  12        int *    lda,
  13        double * e,
  14        double * tau,
  15        double * w,
  16        int *    ldw)
  17 {
  18   int i,iw;
  19   int ti1,ti2,ti3;
  20   double one,zero,minusone,alpha;
  21   const char ch=toupper(*uplo);
  22
  23   one=1.0;
  24   minusone=-1.0;
  25   zero=0.0;
  26
  27   if(*n<=0)
  28     return;
  29
  30   if(ch=='U') {
  31     for(i=*n;i>=(*n-*nb+1);i--) {
  32       iw = i -*n + *nb;
  33
  34       if(i<*n) {
  35         ti1 = *n-i;
  36         ti2 = 1;
  37         /* BLAS */
  38         F77_FUNC(dgemv,DGEMV)("N",&i,&ti1,&minusone, &(a[ i*(*lda) + 0]),lda,&(w[iw*(*ldw)+(i-1)]),
  39                ldw,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
  40         /* BLAS */
  41         F77_FUNC(dgemv,DGEMV)("N",&i,&ti1,&minusone, &(w[ iw*(*ldw) + 0]),ldw,&(a[i*(*lda)+(i-1)]),
  42                lda,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
  43       }
  44
  45       if(i>1) {
  46         /*  Generate elementary reflector H(i) to annihilate
  47          *              A(1:i-2,i)
  48          */
  49         ti1 = i-1;
  50         ti2 = 1;
  51
  52         /* LAPACK */
  53         F77_FUNC(dlarfg,DLARFG)(&ti1,&(a[(i-1)*(*lda)+(i-2)]),&(a[(i-1)*(*lda)+0]),&ti2,&(tau[i-2]));
  54
  55         e[i-2] = a[(i-1)*(*lda)+(i-2)];
  56         a[(i-1)*(*lda)+(i-2)] = 1.0;
  57
  58         /* Compute W(1:i-1,i) */
  59         ti1 = i-1;
  60         ti2 = 1;
  61
  62         /* BLAS */
  63         F77_FUNC(dsymv,DSYMV)("U",&ti1,&one,a,lda,&(a[(i-1)*(*lda)+0]),&ti2,&zero,
  64                &(w[(iw-1)*(*ldw)+0]),&ti2);
  65         if(i<*n) {
  66           ti1 = i-1;
  67           ti2 = *n-i;
  68           ti3 = 1;
  69           /* BLAS */
  70           F77_FUNC(dgemv,DGEMV)("T",&ti1,&ti2,&one,&(w[iw*(*ldw)+0]),ldw,&(a[(i-1)*(*lda)+0]),&ti3,
  71                  &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
  72
  73           /* BLAS */
  74           F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone,&(a[i*(*lda)+0]),lda,&(w[(iw-1)*(*ldw)+i]),&ti3,
  75                  &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
  76
  77           /* BLAS */
  78           F77_FUNC(dgemv,DGEMV)("T",&ti1,&ti2,&one,&(a[i*(*lda)+0]),lda,&(a[(i-1)*(*lda)+0]),&ti3,
  79                  &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
  80
  81           /* BLAS */
  82           F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone,&(w[iw*(*ldw)+0]),ldw,&(w[(iw-1)*(*ldw)+i]),&ti3,
  83                  &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
  84         }
  85
  86         ti1 = i-1;
  87         ti2 = 1;
  88         /* BLAS */
  89         F77_FUNC(dscal,DSCAL)(&ti1,&(tau[i-2]),&(w[(iw-1)*(*ldw)+0]),&ti2);
  90
  91         alpha = -0.5*tau[i-2]*F77_FUNC(ddot,DDOT)(&ti1,&(w[(iw-1)*(*ldw)+0]),&ti2,
  92                                     &(a[(i-1)*(*lda)+0]),&ti2);
  93
  94         /* BLAS */
  95         F77_FUNC(daxpy,DAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+0]),&ti2,&(w[(iw-1)*(*ldw)+0]),&ti2);
  96
  97       }
  98     }
  99   } else {
 100     /* lower */
 101     for(i=1;i<=*nb;i++) {
 102
 103       ti1 = *n-i+1;
 104       ti2 = i-1;
 105       ti3 = 1;
 106       /* BLAS */
 107       F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone, &(a[ i-1 ]),lda,&(w[ i-1 ]),
 108                ldw,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
 109       /* BLAS */
 110       F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone, &(w[ i-1 ]),ldw,&(a[ i-1 ]),
 111                lda,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
 112
 113       if(i<*n) {
 114         ti1 = *n - i;
 115         ti2 = (*n < i+2 ) ? *n : (i+2);
 116         ti3 = 1;
 117         /* LAPACK */
 118         F77_FUNC(dlarfg,DLARFG)(&ti1,&(a[(i-1)*(*lda)+(i)]),&(a[(i-1)*(*lda)+(ti2-1)]),&ti3,&(tau[i-1]));
 119         e[i-1] = a[(i-1)*(*lda)+(i)];
 120         a[(i-1)*(*lda)+(i)] = 1.0;
 121
 122         ti1 = *n - i;
 123         ti2 = 1;
 124         F77_FUNC(dsymv,DSYMV)("L",&ti1,&one,&(a[i*(*lda)+i]),lda,&(a[(i-1)*(*lda)+i]),&ti2,
 125                &zero,&(w[(i-1)*(*ldw)+i]),&ti2);
 126         ti1 = *n - i;
 127         ti2 = i-1;
 128         ti3 = 1;
 129         /* BLAS */
 130         F77_FUNC(dgemv,DGEMV)("T",&ti1,&ti2,&one,&(w[ i ]),ldw,&(a[(i-1)*(*lda)+i]),&ti3,
 131                &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
 132
 133         /* BLAS */
 134         F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone,&(a[ i ]),lda,&(w[(i-1)*(*ldw)+0]),&ti3,
 135                &one,&(w[(i-1)*(*ldw)+i]),&ti3);
 136
 137         /* BLAS */
 138         F77_FUNC(dgemv,DGEMV)("T",&ti1,&ti2,&one,&(a[ i ]),lda,&(a[(i-1)*(*lda)+i]),&ti3,
 139                &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
 140
 141         /* BLAS */
 142         F77_FUNC(dgemv,DGEMV)("N",&ti1,&ti2,&minusone,&(w[ i ]),ldw,&(w[(i-1)*(*ldw)+0]),&ti3,
 143                &one,&(w[(i-1)*(*ldw)+i]),&ti3);
 144
 145         F77_FUNC(dscal,DSCAL)(&ti1,&(tau[i-1]),&(w[(i-1)*(*ldw)+i]),&ti3);
 146         alpha = -0.5*tau[i-1]*F77_FUNC(ddot,DDOT)(&ti1,&(w[(i-1)*(*ldw)+i]),&ti3,
 147                                    &(a[(i-1)*(*lda)+i]),&ti3);
 148
 149         F77_FUNC(daxpy,DAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+i]),&ti3,&(w[(i-1)*(*ldw)+i]),&ti3);
 150       }
 151     }
 152   }
 153   return;
 154 }
 155
 156
 157
 158