src/gmxlib/gmx_lapack/sorgl2.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(sorgl2,SORGL2)(int *m,
   6                         int *n,
   7                         int *k,
   8                         float *a,
   9                         int *lda,
  10                         float *tau,
  11                         float *work,
  12                         int *info)
  13 {
  14     int a_dim1, a_offset, i__1, i__2;
  15     float r__1;
  16
  17     int i__, j, l;
  18
  19     a_dim1 = *lda;
  20     a_offset = 1 + a_dim1;
  21     a -= a_offset;
  22     --tau;
  23     --work;
  24
  25     i__ = (*m > 1) ? *m : 1;
  26
  27     *info = 0;
  28     if (*m < 0) {
  29         *info = -1;
  30     } else if (*n < *m) {
  31         *info = -2;
  32     } else if (*k < 0 || *k > *m) {
  33         *info = -3;
  34     } else if (*lda < i__) {
  35         *info = -5;
  36     }
  37     if (*info != 0) {
  38         return;
  39     }
  40     if (*m <= 0) {
  41         return;
  42     }
  43
  44     if (*k < *m) {
  45         i__1 = *n;
  46         for (j = 1; j <= i__1; ++j) {
  47             i__2 = *m;
  48             for (l = *k + 1; l <= i__2; ++l) {
  49                 a[l + j * a_dim1] = 0.0;
  50             }
  51             if (j > *k && j <= *m) {
  52                 a[j + j * a_dim1] = 1.0;
  53             }
  54         }
  55     }
  56
  57     for (i__ = *k; i__ >= 1; --i__) {
  58         if (i__ < *n) {
  59             if (i__ < *m) {
  60                 a[i__ + i__ * a_dim1] = 1.0;
  61                 i__1 = *m - i__;
  62                 i__2 = *n - i__ + 1;
  63                 F77_FUNC(slarf,SLARF)("R", &i__1, &i__2, &a[i__ + i__ * a_dim1], lda,
  64                &tau[i__], &a[i__ + 1 + i__ * a_dim1], lda, &work[1]);
  65             }
  66             i__1 = *n - i__;
  67             r__1 = -tau[i__];
  68             F77_FUNC(sscal,SSCAL)(&i__1, &r__1, &a[i__ + (i__ + 1) * a_dim1], lda);
  69         }
  70         a[i__ + i__ * a_dim1] = 1.0 - tau[i__];
  71         i__1 = i__ - 1;
  72         for (l = 1; l <= i__1; ++l) {
  73             a[i__ + l * a_dim1] = 0.0;
  74         }
  75     }
  76     return;
  77
  78 }
  79
  80
  81