3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static integer c__2 = 2;
8 static doublereal c_b20 = -1.;
9 static doublereal c_b22 = 1.;
11 /* Subroutine */ int dgetri_(integer *n, doublereal *a, integer *lda, integer
12 *ipiv, doublereal *work, integer *lwork, integer *info)
14 /* System generated locals */
15 integer a_dim1, a_offset, i__1, i__2, i__3;
18 integer i__, j, jb, nb, jj, jp, nn, iws;
19 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
20 integer *, doublereal *, doublereal *, integer *, doublereal *,
21 integer *, doublereal *, doublereal *, integer *),
22 dgemv_(char *, integer *, integer *, doublereal *, doublereal *,
23 integer *, doublereal *, integer *, doublereal *, doublereal *,
26 extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
27 doublereal *, integer *), dtrsm_(char *, char *, char *, char *,
28 integer *, integer *, doublereal *, doublereal *, integer *,
29 doublereal *, integer *), xerbla_(
31 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
32 integer *, integer *);
34 extern /* Subroutine */ int dtrtri_(char *, char *, integer *, doublereal
35 *, integer *, integer *);
40 /* -- LAPACK routine (version 3.1) -- */
41 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
44 /* .. Scalar Arguments .. */
46 /* .. Array Arguments .. */
52 /* DGETRI computes the inverse of a matrix using the LU factorization */
53 /* computed by DGETRF. */
55 /* This method inverts U and then computes inv(A) by solving the system */
56 /* inv(A)*L = inv(U) for inv(A). */
61 /* N (input) INTEGER */
62 /* The order of the matrix A. N >= 0. */
64 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
65 /* On entry, the factors L and U from the factorization */
66 /* A = P*L*U as computed by DGETRF. */
67 /* On exit, if INFO = 0, the inverse of the original matrix A. */
69 /* LDA (input) INTEGER */
70 /* The leading dimension of the array A. LDA >= max(1,N). */
72 /* IPIV (input) INTEGER array, dimension (N) */
73 /* The pivot indices from DGETRF; for 1<=i<=N, row i of the */
74 /* matrix was interchanged with row IPIV(i). */
76 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
77 /* On exit, if INFO=0, then WORK(1) returns the optimal LWORK. */
79 /* LWORK (input) INTEGER */
80 /* The dimension of the array WORK. LWORK >= max(1,N). */
81 /* For optimal performance LWORK >= N*NB, where NB is */
82 /* the optimal blocksize returned by ILAENV. */
84 /* If LWORK = -1, then a workspace query is assumed; the routine */
85 /* only calculates the optimal size of the WORK array, returns */
86 /* this value as the first entry of the WORK array, and no error */
87 /* message related to LWORK is issued by XERBLA. */
89 /* INFO (output) INTEGER */
90 /* = 0: successful exit */
91 /* < 0: if INFO = -i, the i-th argument had an illegal value */
92 /* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is */
93 /* singular and its inverse could not be computed. */
95 /* ===================================================================== */
97 /* .. Parameters .. */
99 /* .. Local Scalars .. */
101 /* .. External Functions .. */
103 /* .. External Subroutines .. */
105 /* .. Intrinsic Functions .. */
107 /* .. Executable Statements .. */
109 /* Test the input parameters. */
111 /* Parameter adjustments */
113 a_offset = 1 + a_dim1;
120 nb = ilaenv_(&c__1, "DGETRI", " ", n, &c_n1, &c_n1, &c_n1);
122 work[1] = (doublereal) lwkopt;
123 lquery = *lwork == -1;
126 } else if (*lda < max(1,*n)) {
128 } else if (*lwork < max(1,*n) && ! lquery) {
133 xerbla_("DGETRI", &i__1);
139 /* Quick return if possible */
145 /* Form inv(U). If INFO > 0 from DTRTRI, then U is singular, */
146 /* and the inverse is not computed. */
148 dtrtri_("Upper", "Non-unit", n, &a[a_offset], lda, info);
155 if (nb > 1 && nb < *n) {
160 nb = *lwork / ldwork;
162 i__1 = 2, i__2 = ilaenv_(&c__2, "DGETRI", " ", n, &c_n1, &c_n1, &
164 nbmin = max(i__1,i__2);
170 /* Solve the equation inv(A)*L = inv(U) for inv(A). */
172 if (nb < nbmin || nb >= *n) {
174 /* Use unblocked code. */
176 for (j = *n; j >= 1; --j) {
178 /* Copy current column of L to WORK and replace with zeros. */
181 for (i__ = j + 1; i__ <= i__1; ++i__) {
182 work[i__] = a[i__ + j * a_dim1];
183 a[i__ + j * a_dim1] = 0.;
187 /* Compute current column of inv(A). */
191 dgemv_("No transpose", n, &i__1, &c_b20, &a[(j + 1) * a_dim1
192 + 1], lda, &work[j + 1], &c__1, &c_b22, &a[j * a_dim1
199 /* Use blocked code. */
201 nn = (*n - 1) / nb * nb + 1;
203 for (j = nn; i__1 < 0 ? j >= 1 : j <= 1; j += i__1) {
205 i__2 = nb, i__3 = *n - j + 1;
208 /* Copy current block column of L to WORK and replace with */
212 for (jj = j; jj <= i__2; ++jj) {
214 for (i__ = jj + 1; i__ <= i__3; ++i__) {
215 work[i__ + (jj - j) * ldwork] = a[i__ + jj * a_dim1];
216 a[i__ + jj * a_dim1] = 0.;
222 /* Compute current block column of inv(A). */
225 i__2 = *n - j - jb + 1;
226 dgemm_("No transpose", "No transpose", n, &jb, &i__2, &c_b20,
227 &a[(j + jb) * a_dim1 + 1], lda, &work[j + jb], &
228 ldwork, &c_b22, &a[j * a_dim1 + 1], lda);
230 dtrsm_("Right", "Lower", "No transpose", "Unit", n, &jb, &c_b22, &
231 work[j], &ldwork, &a[j * a_dim1 + 1], lda);
236 /* Apply column interchanges. */
238 for (j = *n - 1; j >= 1; --j) {
241 dswap_(n, &a[j * a_dim1 + 1], &c__1, &a[jp * a_dim1 + 1], &c__1);
246 work[1] = (doublereal) iws;