src/gmxlib/gmx_lapack/sgetri.c

   1 #include "gmx_blas.h"
   2 #include "gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 void
   6 F77_FUNC(sgetri,SGETRI)(int *n,
   7         float *a,
   8         int *lda,
   9         int *ipiv,
  10         float *work,
  11         int *lwork,
  12         int *info)
  13 {
  14     int a_dim1, a_offset, i__1, i__2, i__3;
  15
  16     int i__, j, jb, nb, jj, jp, nn, iws;
  17     int nbmin;
  18     int ldwork;
  19     int lwkopt;
  20     int c__1 = 1;
  21     float c_b20 = -1.;
  22     float c_b22 = 1.;
  23
  24     a_dim1 = *lda;
  25     a_offset = 1 + a_dim1;
  26     a -= a_offset;
  27     --ipiv;
  28     --work;
  29
  30     *info = 0;
  31     nb = DGETRI_BLOCKSIZE;
  32     lwkopt = *n * nb;
  33     work[1] = (float) lwkopt;
  34
  35     if (*n < 0) {
  36         *info = -1;
  37     } else if (*lda < (*n)) {
  38         *info = -3;
  39     } else if (*lwork < (*n) && *lwork!=-1) {
  40         *info = -6;
  41     }
  42     if (*info != 0) {
  43         i__1 = -(*info);
  44         return;
  45     } else if (*lwork == -1) {
  46         return;
  47     }
  48
  49     if (*n == 0) {
  50         return;
  51     }
  52
  53     F77_FUNC(strtri,STRTRI)("Upper", "Non-unit", n, &a[a_offset], lda, info);
  54     if (*info > 0) {
  55         return;
  56     }
  57
  58     nbmin = 2;
  59     ldwork = *n;
  60     if (nb > 1 && nb < *n) {
  61         i__1 = ldwork * nb;
  62         iws = (i__1>1) ? i__1 : 1;
  63         if (*lwork < iws) {
  64             nb = *lwork / ldwork;
  65             nbmin = DGETRI_MINBLOCKSIZE;
  66         }
  67     } else {
  68         iws = *n;
  69     }
  70
  71     if (nb < nbmin || nb >= *n) {
  72
  73         for (j = *n; j >= 1; --j) {
  74
  75             i__1 = *n;
  76             for (i__ = j + 1; i__ <= i__1; ++i__) {
  77                 work[i__] = a[i__ + j * a_dim1];
  78                 a[i__ + j * a_dim1] = 0.;
  79             }
  80
  81             if (j < *n) {
  82                 i__1 = *n - j;
  83                 F77_FUNC(sgemv,SGEMV)("No transpose", n, &i__1, &c_b20, &a[(j + 1) * a_dim1
  84                         + 1], lda, &work[j + 1], &c__1, &c_b22, &a[j * a_dim1
  85                         + 1], &c__1);
  86             }
  87         }
  88     } else {
  89
  90         nn = (*n - 1) / nb * nb + 1;
  91         i__1 = -nb;
  92         for (j = nn; i__1 < 0 ? j >= 1 : j <= 1; j += i__1) {
  93             i__2 = nb, i__3 = *n - j + 1;
  94             jb = (i__2<i__3) ? i__2 : i__3;
  95
  96             i__2 = j + jb - 1;
  97             for (jj = j; jj <= i__2; ++jj) {
  98                 i__3 = *n;
  99                 for (i__ = jj + 1; i__ <= i__3; ++i__) {
 100                     work[i__ + (jj - j) * ldwork] = a[i__ + jj * a_dim1];
 101                     a[i__ + jj * a_dim1] = 0.;
 102                 }
 103             }
 104
 105             if (j + jb <= *n) {
 106                 i__2 = *n - j - jb + 1;
 107                 F77_FUNC(sgemm,SGEMM)("No transpose", "No transpose", n, &jb, &i__2, &c_b20,
 108                         &a[(j + jb) * a_dim1 + 1], lda, &work[j + jb], &
 109                         ldwork, &c_b22, &a[j * a_dim1 + 1], lda);
 110             }
 111             F77_FUNC(strsm,STRSM)("Right", "Lower", "No transpose", "Unit", n, &jb, &c_b22, &
 112                     work[j], &ldwork, &a[j * a_dim1 + 1], lda);
 113         }
 114     }
 115
 116     for (j = *n - 1; j >= 1; --j) {
 117         jp = ipiv[j];
 118         if (jp != j) {
 119             F77_FUNC(sswap,SSWAP)(n, &a[j * a_dim1 + 1], &c__1, &a[jp * a_dim1 + 1], &c__1);
 120         }
 121     }
 122
 123     work[1] = (float) iws;
 124     return;
 125
 126 }
 127
 128