src/gmxlib/gmx_lapack/sgetrf.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 void
   6 F77_FUNC(sgetrf,SGETRF)(int *m,
   7         int *n,
   8         float *a,
   9         int *lda,
  10         int *ipiv,
  11         int *info)
  12 {
  13   int mindim,jb;
  14   int i,j,k,l;
  15   int iinfo;
  16   float minusone = -1.0;
  17   float one = 1.0;
  18
  19   if(*m<=0 || *n<=0)
  20     return;
  21
  22   *info = 0;
  23
  24   mindim = (*m < *n) ? *m : *n;
  25
  26   if(DGETRF_BLOCKSIZE>=mindim) {
  27
  28     /* unblocked code */
  29     F77_FUNC(sgetf2,SGETF2)(m,n,a,lda,ipiv,info);
  30
  31   } else {
  32
  33     /* blocked case */
  34
  35     for(j=1;j<=mindim;j+=DGETRF_BLOCKSIZE) {
  36       jb = ( DGETRF_BLOCKSIZE < (mindim-j+1)) ? DGETRF_BLOCKSIZE : (mindim-j+1);
  37       /* factor diag. and subdiag blocks and test for singularity */
  38       k = *m-j+1;
  39       F77_FUNC(sgetf2,SGETF2)(&k,&jb,&(a[(j-1)*(*lda)+(j-1)]),lda,&(ipiv[j-1]),&iinfo);
  40
  41       if(*info==0 && iinfo>0)
  42         *info = iinfo + j - 1;
  43
  44       /* adjust pivot indices */
  45       k = (*m < (j+jb-1)) ? *m : (j+jb-1);
  46       for(i=j;i<=k;i++)
  47         ipiv[i-1] += j - 1;
  48
  49       /* Apply to columns 1 throughj j-1 */
  50       k = j - 1;
  51       i = j + jb - 1;
  52       l = 1;
  53       F77_FUNC(slaswp,SLASWP)(&k,a,lda,&j,&i,ipiv,&l);
  54       if((j+jb)<=*n) {
  55         /* Apply to cols. j+jb through n */
  56         k = *n-j-jb+1;
  57         i = j+jb-1;
  58         l = 1;
  59         F77_FUNC(slaswp,SLASWP)(&k,&(a[(j+jb-1)*(*lda)+0]),lda,&j,&i,ipiv,&l);
  60         /* Compute block row of U */
  61         k = *n-j-jb+1;
  62         F77_FUNC(strsm,STRSM)("Left","Lower","No transpose","Unit",&jb,&k,&one,
  63                &(a[(j-1)*(*lda)+(j-1)]),lda,&(a[(j+jb-1)*(*lda)+(j-1)]),lda);
  64
  65         if((j+jb)<=*m) {
  66           /* Update trailing submatrix */
  67           k = *m-j-jb+1;
  68           i = *n-j-jb+1;
  69           F77_FUNC(sgemm,SGEMM)("No transpose","No transpose",&k,&i,&jb,&minusone,
  70                  &(a[(j-1)*(*lda)+(j+jb-1)]),lda,
  71                  &(a[(j+jb-1)*(*lda)+(j-1)]),lda,&one,
  72                  &(a[(j+jb-1)*(*lda)+(j+jb-1)]),lda);
  73         }
  74
  75       }
  76     }
  77   }
  78 }