git.maemo.org Git - opencv/blob - 3rdparty/lapack/dgetri.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   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.;
  10
  11 /* Subroutine */ int dgetri_(integer *n, doublereal *a, integer *lda, integer
  12         *ipiv, doublereal *work, integer *lwork, integer *info)
  13 {
  14     /* System generated locals */
  15     integer a_dim1, a_offset, i__1, i__2, i__3;
  16
  17     /* Local variables */
  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 *,
  24             integer *);
  25     integer nbmin;
  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_(
  30             char *, integer *);
  31     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  32             integer *, integer *);
  33     integer ldwork;
  34     extern /* Subroutine */ int dtrtri_(char *, char *, integer *, doublereal
  35             *, integer *, integer *);
  36     integer lwkopt;
  37     logical lquery;
  38
  39
  40 /*  -- LAPACK routine (version 3.1) -- */
  41 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  42 /*     November 2006 */
  43
  44 /*     .. Scalar Arguments .. */
  45 /*     .. */
  46 /*     .. Array Arguments .. */
  47 /*     .. */
  48
  49 /*  Purpose */
  50 /*  ======= */
  51
  52 /*  DGETRI computes the inverse of a matrix using the LU factorization */
  53 /*  computed by DGETRF. */
  54
  55 /*  This method inverts U and then computes inv(A) by solving the system */
  56 /*  inv(A)*L = inv(U) for inv(A). */
  57
  58 /*  Arguments */
  59 /*  ========= */
  60
  61 /*  N       (input) INTEGER */
  62 /*          The order of the matrix A.  N >= 0. */
  63
  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. */
  68
  69 /*  LDA     (input) INTEGER */
  70 /*          The leading dimension of the array A.  LDA >= max(1,N). */
  71
  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). */
  75
  76 /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
  77 /*          On exit, if INFO=0, then WORK(1) returns the optimal LWORK. */
  78
  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. */
  83
  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. */
  88
  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. */
  94
  95 /*  ===================================================================== */
  96
  97 /*     .. Parameters .. */
  98 /*     .. */
  99 /*     .. Local Scalars .. */
 100 /*     .. */
 101 /*     .. External Functions .. */
 102 /*     .. */
 103 /*     .. External Subroutines .. */
 104 /*     .. */
 105 /*     .. Intrinsic Functions .. */
 106 /*     .. */
 107 /*     .. Executable Statements .. */
 108
 109 /*     Test the input parameters. */
 110
 111     /* Parameter adjustments */
 112     a_dim1 = *lda;
 113     a_offset = 1 + a_dim1;
 114     a -= a_offset;
 115     --ipiv;
 116     --work;
 117
 118     /* Function Body */
 119     *info = 0;
 120     nb = ilaenv_(&c__1, "DGETRI", " ", n, &c_n1, &c_n1, &c_n1);
 121     lwkopt = *n * nb;
 122     work[1] = (doublereal) lwkopt;
 123     lquery = *lwork == -1;
 124     if (*n < 0) {
 125         *info = -1;
 126     } else if (*lda < max(1,*n)) {
 127         *info = -3;
 128     } else if (*lwork < max(1,*n) && ! lquery) {
 129         *info = -6;
 130     }
 131     if (*info != 0) {
 132         i__1 = -(*info);
 133         xerbla_("DGETRI", &i__1);
 134         return 0;
 135     } else if (lquery) {
 136         return 0;
 137     }
 138
 139 /*     Quick return if possible */
 140
 141     if (*n == 0) {
 142         return 0;
 143     }
 144
 145 /*     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular, */
 146 /*     and the inverse is not computed. */
 147
 148     dtrtri_("Upper", "Non-unit", n, &a[a_offset], lda, info);
 149     if (*info > 0) {
 150         return 0;
 151     }
 152
 153     nbmin = 2;
 154     ldwork = *n;
 155     if (nb > 1 && nb < *n) {
 156 /* Computing MAX */
 157         i__1 = ldwork * nb;
 158         iws = max(i__1,1);
 159         if (*lwork < iws) {
 160             nb = *lwork / ldwork;
 161 /* Computing MAX */
 162             i__1 = 2, i__2 = ilaenv_(&c__2, "DGETRI", " ", n, &c_n1, &c_n1, &
 163                     c_n1);
 164             nbmin = max(i__1,i__2);
 165         }
 166     } else {
 167         iws = *n;
 168     }
 169
 170 /*     Solve the equation inv(A)*L = inv(U) for inv(A). */
 171
 172     if (nb < nbmin || nb >= *n) {
 173
 174 /*        Use unblocked code. */
 175
 176         for (j = *n; j >= 1; --j) {
 177
 178 /*           Copy current column of L to WORK and replace with zeros. */
 179
 180             i__1 = *n;
 181             for (i__ = j + 1; i__ <= i__1; ++i__) {
 182                 work[i__] = a[i__ + j * a_dim1];
 183                 a[i__ + j * a_dim1] = 0.;
 184 /* L10: */
 185             }
 186
 187 /*           Compute current column of inv(A). */
 188
 189             if (j < *n) {
 190                 i__1 = *n - j;
 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
 193                         + 1], &c__1);
 194             }
 195 /* L20: */
 196         }
 197     } else {
 198
 199 /*        Use blocked code. */
 200
 201         nn = (*n - 1) / nb * nb + 1;
 202         i__1 = -nb;
 203         for (j = nn; i__1 < 0 ? j >= 1 : j <= 1; j += i__1) {
 204 /* Computing MIN */
 205             i__2 = nb, i__3 = *n - j + 1;
 206             jb = min(i__2,i__3);
 207
 208 /*           Copy current block column of L to WORK and replace with */
 209 /*           zeros. */
 210
 211             i__2 = j + jb - 1;
 212             for (jj = j; jj <= i__2; ++jj) {
 213                 i__3 = *n;
 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.;
 217 /* L30: */
 218                 }
 219 /* L40: */
 220             }
 221
 222 /*           Compute current block column of inv(A). */
 223
 224             if (j + jb <= *n) {
 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);
 229             }
 230             dtrsm_("Right", "Lower", "No transpose", "Unit", n, &jb, &c_b22, &
 231                     work[j], &ldwork, &a[j * a_dim1 + 1], lda);
 232 /* L50: */
 233         }
 234     }
 235
 236 /*     Apply column interchanges. */
 237
 238     for (j = *n - 1; j >= 1; --j) {
 239         jp = ipiv[j];
 240         if (jp != j) {
 241             dswap_(n, &a[j * a_dim1 + 1], &c__1, &a[jp * a_dim1 + 1], &c__1);
 242         }
 243 /* L60: */
 244     }
 245
 246     work[1] = (doublereal) iws;
 247     return 0;
 248
 249 /*     End of DGETRI */
 250
 251 } /* dgetri_ */