src/gromacs/linearalgebra/gmx_lapack/dgetf2.c

   1 #include <math.h>
   2 #include "types/simple.h"
   3
   4 #include "../gmx_blas.h"
   5 #include "../gmx_lapack.h"
   6
   7
   8 void
   9 F77_FUNC(dgetf2,DGETF2)(int *m,
  10         int *n,
  11         double *a,
  12         int *lda,
  13         int *ipiv,
  14         int *info)
  15 {
  16   int j,jp,k,t1,t2,t3;
  17   double one,minusone;
  18   double tmp;
  19
  20   one = 1.0;
  21   minusone = -1.0;
  22
  23   if(*m<=0 || *n<=0)
  24     return;
  25
  26   k = (*m < *n) ? *m : *n;
  27   for(j=1;j<=k;j++) {
  28     t1 = *m-j+1;
  29     t2 = 1;
  30     jp = j - 1 + F77_FUNC(idamax,IDAMAX)(&t1,&(a[(j-1)*(*lda)+(j-1)]),&t2);
  31     ipiv[j-1] = jp;
  32     if( fabs(a[(j-1)*(*lda)+(jp-1)])>GMX_DOUBLE_MIN ) {
  33       if(jp != j)
  34         F77_FUNC(dswap,DSWAP)(n,&(a[ j-1 ]),lda,&(a[ jp-1 ]),lda);
  35
  36       if(j<*m) {
  37         t1 = *m-j;
  38         t2 = 1;
  39         tmp = 1.0/a[(j-1)*(*lda)+(j-1)];
  40         F77_FUNC(dscal,DSCAL)(&t1,&tmp,&(a[(j-1)*(*lda)+(j)]),&t2);
  41       }
  42     } else {
  43       *info = j;
  44     }
  45
  46     if(j<k) {
  47       t1 = *m-j;
  48       t2 = *n-j;
  49       t3 = 1;
  50       F77_FUNC(dger,DGER)(&t1,&t2,&minusone,&(a[(j-1)*(*lda)+(j)]),&t3,
  51             &(a[(j)*(*lda)+(j-1)]),lda, &(a[(j)*(*lda)+(j)]),lda);
  52     }
  53   }
  54   return;
  55 }