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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dsytri_(char *uplo, integer *n, doublereal *a, integer *
   4         lda, integer *ipiv, doublereal *work, integer *info)
   5 {
   6 /*  -- LAPACK routine (version 3.0) --
   7        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
   8        Courant Institute, Argonne National Lab, and Rice University
   9        March 31, 1993
  10
  11
  12     Purpose
  13     =======
  14
  15     DSYTRI computes the inverse of a real symmetric indefinite matrix
  16     A using the factorization A = U*D*U**T or A = L*D*L**T computed by
  17     DSYTRF.
  18
  19     Arguments
  20     =========
  21
  22     UPLO    (input) CHARACTER*1
  23             Specifies whether the details of the factorization are stored
  24             as an upper or lower triangular matrix.
  25             = 'U':  Upper triangular, form is A = U*D*U**T;
  26             = 'L':  Lower triangular, form is A = L*D*L**T.
  27
  28     N       (input) INTEGER
  29             The order of the matrix A.  N >= 0.
  30
  31     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  32             On entry, the block diagonal matrix D and the multipliers
  33             used to obtain the factor U or L as computed by DSYTRF.
  34
  35             On exit, if INFO = 0, the (symmetric) inverse of the original
  36             matrix.  If UPLO = 'U', the upper triangular part of the
  37             inverse is formed and the part of A below the diagonal is not
  38             referenced; if UPLO = 'L' the lower triangular part of the
  39             inverse is formed and the part of A above the diagonal is
  40             not referenced.
  41
  42     LDA     (input) INTEGER
  43             The leading dimension of the array A.  LDA >= max(1,N).
  44
  45     IPIV    (input) INTEGER array, dimension (N)
  46             Details of the interchanges and the block structure of D
  47             as determined by DSYTRF.
  48
  49     WORK    (workspace) DOUBLE PRECISION array, dimension (N)
  50
  51     INFO    (output) INTEGER
  52             = 0: successful exit
  53             < 0: if INFO = -i, the i-th argument had an illegal value
  54             > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
  55                  inverse could not be computed.
  56
  57     =====================================================================
  58
  59
  60        Test the input parameters.
  61
  62        Parameter adjustments */
  63     /* Table of constant values */
  64     static integer c__1 = 1;
  65     static doublereal c_b11 = -1.;
  66     static doublereal c_b13 = 0.;
  67
  68     /* System generated locals */
  69     integer a_dim1, a_offset, i__1;
  70     doublereal d__1;
  71     /* Local variables */
  72     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  73             integer *);
  74     static doublereal temp, akkp1, d__;
  75     static integer k;
  76     static doublereal t;
  77     extern logical lsame_(char *, char *);
  78     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  79             doublereal *, integer *), dswap_(integer *, doublereal *, integer
  80             *, doublereal *, integer *);
  81     static integer kstep;
  82     static logical upper;
  83     extern /* Subroutine */ int dsymv_(char *, integer *, doublereal *,
  84             doublereal *, integer *, doublereal *, integer *, doublereal *,
  85             doublereal *, integer *);
  86     static doublereal ak;
  87     static integer kp;
  88     extern /* Subroutine */ int xerbla_(char *, integer *);
  89     static doublereal akp1;
  90 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
  91
  92
  93     a_dim1 = *lda;
  94     a_offset = 1 + a_dim1 * 1;
  95     a -= a_offset;
  96     --ipiv;
  97     --work;
  98
  99     /* Function Body */
 100     *info = 0;
 101     upper = lsame_(uplo, "U");
 102     if (! upper && ! lsame_(uplo, "L")) {
 103         *info = -1;
 104     } else if (*n < 0) {
 105         *info = -2;
 106     } else if (*lda < max(1,*n)) {
 107         *info = -4;
 108     }
 109     if (*info != 0) {
 110         i__1 = -(*info);
 111         xerbla_("DSYTRI", &i__1);
 112         return 0;
 113     }
 114
 115 /*     Quick return if possible */
 116
 117     if (*n == 0) {
 118         return 0;
 119     }
 120
 121 /*     Check that the diagonal matrix D is nonsingular. */
 122
 123     if (upper) {
 124
 125 /*        Upper triangular storage: examine D from bottom to top */
 126
 127         for (*info = *n; *info >= 1; --(*info)) {
 128             if (ipiv[*info] > 0 && a_ref(*info, *info) == 0.) {
 129                 return 0;
 130             }
 131 /* L10: */
 132         }
 133     } else {
 134
 135 /*        Lower triangular storage: examine D from top to bottom. */
 136
 137         i__1 = *n;
 138         for (*info = 1; *info <= i__1; ++(*info)) {
 139             if (ipiv[*info] > 0 && a_ref(*info, *info) == 0.) {
 140                 return 0;
 141             }
 142 /* L20: */
 143         }
 144     }
 145     *info = 0;
 146
 147     if (upper) {
 148
 149 /*        Compute inv(A) from the factorization A = U*D*U'.
 150
 151           K is the main loop index, increasing from 1 to N in steps of
 152           1 or 2, depending on the size of the diagonal blocks. */
 153
 154         k = 1;
 155 L30:
 156
 157 /*        If K > N, exit from loop. */
 158
 159         if (k > *n) {
 160             goto L40;
 161         }
 162
 163         if (ipiv[k] > 0) {
 164
 165 /*           1 x 1 diagonal block
 166
 167              Invert the diagonal block. */
 168
 169             a_ref(k, k) = 1. / a_ref(k, k);
 170
 171 /*           Compute column K of the inverse. */
 172
 173             if (k > 1) {
 174                 i__1 = k - 1;
 175                 dcopy_(&i__1, &a_ref(1, k), &c__1, &work[1], &c__1);
 176                 i__1 = k - 1;
 177                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
 178                         c__1, &c_b13, &a_ref(1, k), &c__1);
 179                 i__1 = k - 1;
 180                 a_ref(k, k) = a_ref(k, k) - ddot_(&i__1, &work[1], &c__1, &
 181                         a_ref(1, k), &c__1);
 182             }
 183             kstep = 1;
 184         } else {
 185
 186 /*           2 x 2 diagonal block
 187
 188              Invert the diagonal block. */
 189
 190             t = (d__1 = a_ref(k, k + 1), abs(d__1));
 191             ak = a_ref(k, k) / t;
 192             akp1 = a_ref(k + 1, k + 1) / t;
 193             akkp1 = a_ref(k, k + 1) / t;
 194             d__ = t * (ak * akp1 - 1.);
 195             a_ref(k, k) = akp1 / d__;
 196             a_ref(k + 1, k + 1) = ak / d__;
 197             a_ref(k, k + 1) = -akkp1 / d__;
 198
 199 /*           Compute columns K and K+1 of the inverse. */
 200
 201             if (k > 1) {
 202                 i__1 = k - 1;
 203                 dcopy_(&i__1, &a_ref(1, k), &c__1, &work[1], &c__1);
 204                 i__1 = k - 1;
 205                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
 206                         c__1, &c_b13, &a_ref(1, k), &c__1);
 207                 i__1 = k - 1;
 208                 a_ref(k, k) = a_ref(k, k) - ddot_(&i__1, &work[1], &c__1, &
 209                         a_ref(1, k), &c__1);
 210                 i__1 = k - 1;
 211                 a_ref(k, k + 1) = a_ref(k, k + 1) - ddot_(&i__1, &a_ref(1, k),
 212                          &c__1, &a_ref(1, k + 1), &c__1);
 213                 i__1 = k - 1;
 214                 dcopy_(&i__1, &a_ref(1, k + 1), &c__1, &work[1], &c__1);
 215                 i__1 = k - 1;
 216                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
 217                         c__1, &c_b13, &a_ref(1, k + 1), &c__1);
 218                 i__1 = k - 1;
 219                 a_ref(k + 1, k + 1) = a_ref(k + 1, k + 1) - ddot_(&i__1, &
 220                         work[1], &c__1, &a_ref(1, k + 1), &c__1);
 221             }
 222             kstep = 2;
 223         }
 224
 225         kp = (i__1 = ipiv[k], abs(i__1));
 226         if (kp != k) {
 227
 228 /*           Interchange rows and columns K and KP in the leading
 229              submatrix A(1:k+1,1:k+1) */
 230
 231             i__1 = kp - 1;
 232             dswap_(&i__1, &a_ref(1, k), &c__1, &a_ref(1, kp), &c__1);
 233             i__1 = k - kp - 1;
 234             dswap_(&i__1, &a_ref(kp + 1, k), &c__1, &a_ref(kp, kp + 1), lda);
 235             temp = a_ref(k, k);
 236             a_ref(k, k) = a_ref(kp, kp);
 237             a_ref(kp, kp) = temp;
 238             if (kstep == 2) {
 239                 temp = a_ref(k, k + 1);
 240                 a_ref(k, k + 1) = a_ref(kp, k + 1);
 241                 a_ref(kp, k + 1) = temp;
 242             }
 243         }
 244
 245         k += kstep;
 246         goto L30;
 247 L40:
 248
 249         ;
 250     } else {
 251
 252 /*        Compute inv(A) from the factorization A = L*D*L'.
 253
 254           K is the main loop index, increasing from 1 to N in steps of
 255           1 or 2, depending on the size of the diagonal blocks. */
 256
 257         k = *n;
 258 L50:
 259
 260 /*        If K < 1, exit from loop. */
 261
 262         if (k < 1) {
 263             goto L60;
 264         }
 265
 266         if (ipiv[k] > 0) {
 267
 268 /*           1 x 1 diagonal block
 269
 270              Invert the diagonal block. */
 271
 272             a_ref(k, k) = 1. / a_ref(k, k);
 273
 274 /*           Compute column K of the inverse. */
 275
 276             if (k < *n) {
 277                 i__1 = *n - k;
 278                 dcopy_(&i__1, &a_ref(k + 1, k), &c__1, &work[1], &c__1);
 279                 i__1 = *n - k;
 280                 dsymv_(uplo, &i__1, &c_b11, &a_ref(k + 1, k + 1), lda, &work[
 281                         1], &c__1, &c_b13, &a_ref(k + 1, k), &c__1)
 282                         ;
 283                 i__1 = *n - k;
 284                 a_ref(k, k) = a_ref(k, k) - ddot_(&i__1, &work[1], &c__1, &
 285                         a_ref(k + 1, k), &c__1);
 286             }
 287             kstep = 1;
 288         } else {
 289
 290 /*           2 x 2 diagonal block
 291
 292              Invert the diagonal block. */
 293
 294             t = (d__1 = a_ref(k, k - 1), abs(d__1));
 295             ak = a_ref(k - 1, k - 1) / t;
 296             akp1 = a_ref(k, k) / t;
 297             akkp1 = a_ref(k, k - 1) / t;
 298             d__ = t * (ak * akp1 - 1.);
 299             a_ref(k - 1, k - 1) = akp1 / d__;
 300             a_ref(k, k) = ak / d__;
 301             a_ref(k, k - 1) = -akkp1 / d__;
 302
 303 /*           Compute columns K-1 and K of the inverse. */
 304
 305             if (k < *n) {
 306                 i__1 = *n - k;
 307                 dcopy_(&i__1, &a_ref(k + 1, k), &c__1, &work[1], &c__1);
 308                 i__1 = *n - k;
 309                 dsymv_(uplo, &i__1, &c_b11, &a_ref(k + 1, k + 1), lda, &work[
 310                         1], &c__1, &c_b13, &a_ref(k + 1, k), &c__1)
 311                         ;
 312                 i__1 = *n - k;
 313                 a_ref(k, k) = a_ref(k, k) - ddot_(&i__1, &work[1], &c__1, &
 314                         a_ref(k + 1, k), &c__1);
 315                 i__1 = *n - k;
 316                 a_ref(k, k - 1) = a_ref(k, k - 1) - ddot_(&i__1, &a_ref(k + 1,
 317                          k), &c__1, &a_ref(k + 1, k - 1), &c__1);
 318                 i__1 = *n - k;
 319                 dcopy_(&i__1, &a_ref(k + 1, k - 1), &c__1, &work[1], &c__1);
 320                 i__1 = *n - k;
 321                 dsymv_(uplo, &i__1, &c_b11, &a_ref(k + 1, k + 1), lda, &work[
 322                         1], &c__1, &c_b13, &a_ref(k + 1, k - 1), &c__1);
 323                 i__1 = *n - k;
 324                 a_ref(k - 1, k - 1) = a_ref(k - 1, k - 1) - ddot_(&i__1, &
 325                         work[1], &c__1, &a_ref(k + 1, k - 1), &c__1);
 326             }
 327             kstep = 2;
 328         }
 329
 330         kp = (i__1 = ipiv[k], abs(i__1));
 331         if (kp != k) {
 332
 333 /*           Interchange rows and columns K and KP in the trailing
 334              submatrix A(k-1:n,k-1:n) */
 335
 336             if (kp < *n) {
 337                 i__1 = *n - kp;
 338                 dswap_(&i__1, &a_ref(kp + 1, k), &c__1, &a_ref(kp + 1, kp), &
 339                         c__1);
 340             }
 341             i__1 = kp - k - 1;
 342             dswap_(&i__1, &a_ref(k + 1, k), &c__1, &a_ref(kp, k + 1), lda);
 343             temp = a_ref(k, k);
 344             a_ref(k, k) = a_ref(kp, kp);
 345             a_ref(kp, kp) = temp;
 346             if (kstep == 2) {
 347                 temp = a_ref(k, k - 1);
 348                 a_ref(k, k - 1) = a_ref(kp, k - 1);
 349                 a_ref(kp, k - 1) = temp;
 350             }
 351         }
 352
 353         k -= kstep;
 354         goto L50;
 355 L60:
 356         ;
 357     }
 358
 359     return 0;
 360
 361 /*     End of DSYTRI */
 362
 363 } /* dsytri_ */
 364
 365 #undef a_ref
 366
 367