2 #include "gmx_lapack.h"
3 #include "lapack_limits.h"
6 F77_FUNC(dgetri,DGETRI)(int *n,
14 int a_dim1, a_offset, i__1, i__2, i__3;
16 int i__, j, jb, nb, jj, jp, nn, iws;
25 a_offset = 1 + a_dim1;
31 nb = DGETRI_BLOCKSIZE;
33 work[1] = (double) lwkopt;
37 } else if (*lda < (*n)) {
39 } else if (*lwork < (*n) && *lwork!=-1) {
45 } else if (*lwork == -1) {
53 F77_FUNC(dtrtri,DTRTRI)("Upper", "Non-unit", n, &a[a_offset], lda, info);
60 if (nb > 1 && nb < *n) {
62 iws = (i__1>1) ? i__1 : 1;
65 nbmin = DGETRI_MINBLOCKSIZE;
71 if (nb < nbmin || nb >= *n) {
73 for (j = *n; j >= 1; --j) {
76 for (i__ = j + 1; i__ <= i__1; ++i__) {
77 work[i__] = a[i__ + j * a_dim1];
78 a[i__ + j * a_dim1] = 0.;
83 F77_FUNC(dgemv,DGEMV)("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
90 nn = (*n - 1) / nb * nb + 1;
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;
97 for (jj = j; jj <= i__2; ++jj) {
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.;
106 i__2 = *n - j - jb + 1;
107 F77_FUNC(dgemm,DGEMM)("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);
111 F77_FUNC(dtrsm,DTRSM)("Right", "Lower", "No transpose", "Unit", n, &jb, &c_b22, &
112 work[j], &ldwork, &a[j * a_dim1 + 1], lda);
116 for (j = *n - 1; j >= 1; --j) {
119 F77_FUNC(dswap,DSWAP)(n, &a[j * a_dim1 + 1], &c__1, &a[jp * a_dim1 + 1], &c__1);
123 work[1] = (double) iws;