src/gromacs/linearalgebra/gmx_lapack/dtrtri.c

   1 #include <math.h>
   2 #include "../gmx_blas.h"
   3 #include "../gmx_lapack.h"
   4 #include "lapack_limits.h"
   5
   6 #include "gromacs/utility/real.h"
   7
   8 void
   9 F77_FUNC(dtrtri,DTRTRI)(const char *uplo,
  10         const char *diag,
  11         int *n,
  12         double *a,
  13         int *lda,
  14         int *info)
  15 {
  16     int a_dim1, a_offset, i__1, i__3, i__4, i__5;
  17     int j, jb, nb, nn;
  18     double c_b18 = 1.;
  19     double c_b22 = -1.;
  20
  21     int upper;
  22     int nounit;
  23
  24     a_dim1 = *lda;
  25     a_offset = 1 + a_dim1;
  26     a -= a_offset;
  27
  28     *info = 0;
  29     upper = (*uplo=='U' || *uplo=='u');
  30     nounit = (*diag=='N' || *diag=='n');
  31
  32     if (*info != 0) {
  33         i__1 = -(*info);
  34         return;
  35     }
  36
  37     if (*n == 0) {
  38         return;
  39     }
  40
  41     if (nounit) {
  42         i__1 = *n;
  43         for (*info = 1; *info <= i__1; ++(*info)) {
  44             if (fabs(a[*info + *info * a_dim1])<GMX_DOUBLE_MIN) {
  45                 return;
  46             }
  47         }
  48         *info = 0;
  49     }
  50
  51     nb = DTRTRI_BLOCKSIZE;
  52     if (nb <= 1 || nb >= *n) {
  53
  54         F77_FUNC(dtrti2,DTRTI2)(uplo, diag, n, &a[a_offset], lda, info);
  55     } else {
  56
  57         if (upper) {
  58
  59             i__1 = *n;
  60             i__3 = nb;
  61             for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
  62                 i__4 = nb, i__5 = *n - j + 1;
  63                 jb = (i__4<i__5) ? i__4 : i__5;
  64
  65                 i__4 = j - 1;
  66                 F77_FUNC(dtrmm,DTRMM)("Left", "Upper", "No transpose", diag, &i__4, &jb, &
  67                         c_b18, &a[a_offset], lda, &a[j * a_dim1 + 1], lda);
  68                 i__4 = j - 1;
  69                 F77_FUNC(dtrsm,DTRSM)("Right", "Upper", "No transpose", diag, &i__4, &jb, &
  70                         c_b22, &a[j + j * a_dim1], lda, &a[j * a_dim1 + 1],
  71                         lda);
  72
  73                 F77_FUNC(dtrti2,DTRTI2)("Upper", diag, &jb, &a[j + j * a_dim1], lda, info);
  74             }
  75         } else {
  76
  77             nn = (*n - 1) / nb * nb + 1;
  78             i__3 = -nb;
  79             for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) {
  80                 i__1 = nb, i__4 = *n - j + 1;
  81                 jb = (i__1<i__4) ? i__1 : i__4;
  82                 if (j + jb <= *n) {
  83
  84                     i__1 = *n - j - jb + 1;
  85                     F77_FUNC(dtrmm,DTRMM)("Left", "Lower", "No transpose", diag, &i__1, &jb,
  86                             &c_b18, &a[j + jb + (j + jb) * a_dim1], lda, &a[j
  87                             + jb + j * a_dim1], lda);
  88                     i__1 = *n - j - jb + 1;
  89                     F77_FUNC(dtrsm,DTRSM)("Right", "Lower", "No transpose", diag, &i__1, &jb,
  90                              &c_b22, &a[j + j * a_dim1], lda, &a[j + jb + j *
  91                             a_dim1], lda);
  92                 }
  93
  94                 F77_FUNC(dtrti2,DTRTI2)("Lower", diag, &jb, &a[j + j * a_dim1], lda, info);
  95             }
  96         }
  97     }
  98     return;
  99 }
 100
 101