src/gromacs/linearalgebra/gmx_lapack/dorml2.c

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