src/gmxlib/gmx_lapack/sorml2.c

   1 #include <ctype.h>
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(sorml2,SORML2)(const char *side,
   6                         const char *trans,
   7                         int *m,
   8                         int *n,
   9                         int *k,
  10                         float *a,
  11                         int *lda,
  12                         float *tau,
  13                         float *c,
  14                         int *ldc,
  15                         float *work,
  16                         int *info)
  17 {
  18   const char xside=toupper(*side);
  19   const char xtrans=toupper(*trans);
  20   int i,i1,i2,i3,ni,mi,ic,jc;
  21   float aii;
  22
  23   if(*m<=0 || *n<=0 || *k<=0)
  24     return;
  25
  26   ic = jc = 0;
  27
  28   if((xside=='L' && xtrans=='N') || (xside!='L' && xtrans!='N')) {
  29     i1 = 0;
  30     i2 = *k;
  31     i3 = 1;
  32   } else {
  33     i1 = *k-1;
  34     i2 = -1;
  35     i3 = -1;
  36   }
  37
  38   if(xside=='L') {
  39     ni = *n;
  40     jc = 0;
  41   } else {
  42     mi = *m;
  43     ic = 0;
  44   }
  45
  46   for(i=i1;i!=i2;i+=i3) {
  47     if(xside=='L') {
  48       mi = *m - i;
  49       ic = i;
  50     } else {
  51       ni = *n - i;
  52       jc = i;
  53     }
  54     aii = a[i*(*lda)+i];
  55     a[i*(*lda)+i] = 1.0;
  56     F77_FUNC(slarf,SLARF)(side,&mi,&ni,&(a[i*(*lda)+i]),lda,tau+i,
  57            &(c[jc*(*ldc)+ic]),ldc,work);
  58     a[i*(*lda)+i] = aii;
  59   }
  60   return;
  61 }
  62