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 real c_b20 = -1.f;
9 static real c_b22 = 1.f;
11 /* Subroutine */ int sgetri_(integer *n, real *a, integer *lda, integer *ipiv,
12 real *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, nbmin;
19 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
20 integer *, real *, real *, integer *, real *, integer *, real *,
21 real *, integer *), sgemv_(char *, integer *,
22 integer *, real *, real *, integer *, real *, integer *, real *,
23 real *, integer *), sswap_(integer *, real *, integer *,
24 real *, integer *), strsm_(char *, char *, char *, char *,
25 integer *, integer *, real *, real *, integer *, real *, integer *
26 ), xerbla_(char *, integer *);
27 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
28 integer *, integer *);
29 integer ldwork, lwkopt;
31 extern /* Subroutine */ int strtri_(char *, char *, integer *, real *,
32 integer *, integer *);
35 /* -- LAPACK routine (version 3.1) -- */
36 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
39 /* .. Scalar Arguments .. */
41 /* .. Array Arguments .. */
47 /* SGETRI computes the inverse of a matrix using the LU factorization */
48 /* computed by SGETRF. */
50 /* This method inverts U and then computes inv(A) by solving the system */
51 /* inv(A)*L = inv(U) for inv(A). */
56 /* N (input) INTEGER */
57 /* The order of the matrix A. N >= 0. */
59 /* A (input/output) REAL array, dimension (LDA,N) */
60 /* On entry, the factors L and U from the factorization */
61 /* A = P*L*U as computed by SGETRF. */
62 /* On exit, if INFO = 0, the inverse of the original matrix A. */
64 /* LDA (input) INTEGER */
65 /* The leading dimension of the array A. LDA >= max(1,N). */
67 /* IPIV (input) INTEGER array, dimension (N) */
68 /* The pivot indices from SGETRF; for 1<=i<=N, row i of the */
69 /* matrix was interchanged with row IPIV(i). */
71 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
72 /* On exit, if INFO=0, then WORK(1) returns the optimal LWORK. */
74 /* LWORK (input) INTEGER */
75 /* The dimension of the array WORK. LWORK >= max(1,N). */
76 /* For optimal performance LWORK >= N*NB, where NB is */
77 /* the optimal blocksize returned by ILAENV. */
79 /* If LWORK = -1, then a workspace query is assumed; the routine */
80 /* only calculates the optimal size of the WORK array, returns */
81 /* this value as the first entry of the WORK array, and no error */
82 /* message related to LWORK is issued by XERBLA. */
84 /* INFO (output) INTEGER */
85 /* = 0: successful exit */
86 /* < 0: if INFO = -i, the i-th argument had an illegal value */
87 /* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is */
88 /* singular and its inverse could not be computed. */
90 /* ===================================================================== */
92 /* .. Parameters .. */
94 /* .. Local Scalars .. */
96 /* .. External Functions .. */
98 /* .. External Subroutines .. */
100 /* .. Intrinsic Functions .. */
102 /* .. Executable Statements .. */
104 /* Test the input parameters. */
106 /* Parameter adjustments */
108 a_offset = 1 + a_dim1;
115 nb = ilaenv_(&c__1, "SGETRI", " ", n, &c_n1, &c_n1, &c_n1);
117 work[1] = (real) lwkopt;
118 lquery = *lwork == -1;
121 } else if (*lda < max(1,*n)) {
123 } else if (*lwork < max(1,*n) && ! lquery) {
128 xerbla_("SGETRI", &i__1);
134 /* Quick return if possible */
140 /* Form inv(U). If INFO > 0 from STRTRI, then U is singular, */
141 /* and the inverse is not computed. */
143 strtri_("Upper", "Non-unit", n, &a[a_offset], lda, info);
150 if (nb > 1 && nb < *n) {
155 nb = *lwork / ldwork;
157 i__1 = 2, i__2 = ilaenv_(&c__2, "SGETRI", " ", n, &c_n1, &c_n1, &
159 nbmin = max(i__1,i__2);
165 /* Solve the equation inv(A)*L = inv(U) for inv(A). */
167 if (nb < nbmin || nb >= *n) {
169 /* Use unblocked code. */
171 for (j = *n; j >= 1; --j) {
173 /* Copy current column of L to WORK and replace with zeros. */
176 for (i__ = j + 1; i__ <= i__1; ++i__) {
177 work[i__] = a[i__ + j * a_dim1];
178 a[i__ + j * a_dim1] = 0.f;
182 /* Compute current column of inv(A). */
186 sgemv_("No transpose", n, &i__1, &c_b20, &a[(j + 1) * a_dim1
187 + 1], lda, &work[j + 1], &c__1, &c_b22, &a[j * a_dim1
194 /* Use blocked code. */
196 nn = (*n - 1) / nb * nb + 1;
198 for (j = nn; i__1 < 0 ? j >= 1 : j <= 1; j += i__1) {
200 i__2 = nb, i__3 = *n - j + 1;
203 /* Copy current block column of L to WORK and replace with */
207 for (jj = j; jj <= i__2; ++jj) {
209 for (i__ = jj + 1; i__ <= i__3; ++i__) {
210 work[i__ + (jj - j) * ldwork] = a[i__ + jj * a_dim1];
211 a[i__ + jj * a_dim1] = 0.f;
217 /* Compute current block column of inv(A). */
220 i__2 = *n - j - jb + 1;
221 sgemm_("No transpose", "No transpose", n, &jb, &i__2, &c_b20,
222 &a[(j + jb) * a_dim1 + 1], lda, &work[j + jb], &
223 ldwork, &c_b22, &a[j * a_dim1 + 1], lda);
225 strsm_("Right", "Lower", "No transpose", "Unit", n, &jb, &c_b22, &
226 work[j], &ldwork, &a[j * a_dim1 + 1], lda);
231 /* Apply column interchanges. */
233 for (j = *n - 1; j >= 1; --j) {
236 sswap_(n, &a[j * a_dim1 + 1], &c__1, &a[jp * a_dim1 + 1], &c__1);
241 work[1] = (real) iws;