swiftnav: dsytri.c Source File

Go to the documentation of this file.
00001 /* dsytri.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static doublereal c_b11 = -1.;
00020 static doublereal c_b13 = 0.;
00021 
00022 /* Subroutine */ int dsytri_(char *uplo, integer *n, doublereal *a, integer *
00023         lda, integer *ipiv, doublereal *work, integer *info)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, i__1;
00027     doublereal d__1;
00028 
00029     /* Local variables */
00030     doublereal d__;
00031     integer k;
00032     doublereal t, ak;
00033     integer kp;
00034     doublereal akp1;
00035     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
00036             integer *);
00037     doublereal temp, akkp1;
00038     extern logical lsame_(char *, char *);
00039     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
00040             doublereal *, integer *), dswap_(integer *, doublereal *, integer 
00041             *, doublereal *, integer *);
00042     integer kstep;
00043     logical upper;
00044     extern /* Subroutine */ int dsymv_(char *, integer *, doublereal *, 
00045             doublereal *, integer *, doublereal *, integer *, doublereal *, 
00046             doublereal *, integer *), xerbla_(char *, integer *);
00047 
00048 
00049 /*  -- LAPACK routine (version 3.2) -- */
00050 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00051 /*     November 2006 */
00052 
00053 /*     .. Scalar Arguments .. */
00054 /*     .. */
00055 /*     .. Array Arguments .. */
00056 /*     .. */
00057 
00058 /*  Purpose */
00059 /*  ======= */
00060 
00061 /*  DSYTRI computes the inverse of a real symmetric indefinite matrix */
00062 /*  A using the factorization A = U*D*U**T or A = L*D*L**T computed by */
00063 /*  DSYTRF. */
00064 
00065 /*  Arguments */
00066 /*  ========= */
00067 
00068 /*  UPLO    (input) CHARACTER*1 */
00069 /*          Specifies whether the details of the factorization are stored */
00070 /*          as an upper or lower triangular matrix. */
00071 /*          = 'U':  Upper triangular, form is A = U*D*U**T; */
00072 /*          = 'L':  Lower triangular, form is A = L*D*L**T. */
00073 
00074 /*  N       (input) INTEGER */
00075 /*          The order of the matrix A.  N >= 0. */
00076 
00077 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
00078 /*          On entry, the block diagonal matrix D and the multipliers */
00079 /*          used to obtain the factor U or L as computed by DSYTRF. */
00080 
00081 /*          On exit, if INFO = 0, the (symmetric) inverse of the original */
00082 /*          matrix.  If UPLO = 'U', the upper triangular part of the */
00083 /*          inverse is formed and the part of A below the diagonal is not */
00084 /*          referenced; if UPLO = 'L' the lower triangular part of the */
00085 /*          inverse is formed and the part of A above the diagonal is */
00086 /*          not referenced. */
00087 
00088 /*  LDA     (input) INTEGER */
00089 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00090 
00091 /*  IPIV    (input) INTEGER array, dimension (N) */
00092 /*          Details of the interchanges and the block structure of D */
00093 /*          as determined by DSYTRF. */
00094 
00095 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (N) */
00096 
00097 /*  INFO    (output) INTEGER */
00098 /*          = 0: successful exit */
00099 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00100 /*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its */
00101 /*               inverse could not be computed. */
00102 
00103 /*  ===================================================================== */
00104 
00105 /*     .. Parameters .. */
00106 /*     .. */
00107 /*     .. Local Scalars .. */
00108 /*     .. */
00109 /*     .. External Functions .. */
00110 /*     .. */
00111 /*     .. External Subroutines .. */
00112 /*     .. */
00113 /*     .. Intrinsic Functions .. */
00114 /*     .. */
00115 /*     .. Executable Statements .. */
00116 
00117 /*     Test the input parameters. */
00118 
00119     /* Parameter adjustments */
00120     a_dim1 = *lda;
00121     a_offset = 1 + a_dim1;
00122     a -= a_offset;
00123     --ipiv;
00124     --work;
00125 
00126     /* Function Body */
00127     *info = 0;
00128     upper = lsame_(uplo, "U");
00129     if (! upper && ! lsame_(uplo, "L")) {
00130         *info = -1;
00131     } else if (*n < 0) {
00132         *info = -2;
00133     } else if (*lda < max(1,*n)) {
00134         *info = -4;
00135     }
00136     if (*info != 0) {
00137         i__1 = -(*info);
00138         xerbla_("DSYTRI", &i__1);
00139         return 0;
00140     }
00141 
00142 /*     Quick return if possible */
00143 
00144     if (*n == 0) {
00145         return 0;
00146     }
00147 
00148 /*     Check that the diagonal matrix D is nonsingular. */
00149 
00150     if (upper) {
00151 
00152 /*        Upper triangular storage: examine D from bottom to top */
00153 
00154         for (*info = *n; *info >= 1; --(*info)) {
00155             if (ipiv[*info] > 0 && a[*info + *info * a_dim1] == 0.) {
00156                 return 0;
00157             }
00158 /* L10: */
00159         }
00160     } else {
00161 
00162 /*        Lower triangular storage: examine D from top to bottom. */
00163 
00164         i__1 = *n;
00165         for (*info = 1; *info <= i__1; ++(*info)) {
00166             if (ipiv[*info] > 0 && a[*info + *info * a_dim1] == 0.) {
00167                 return 0;
00168             }
00169 /* L20: */
00170         }
00171     }
00172     *info = 0;
00173 
00174     if (upper) {
00175 
00176 /*        Compute inv(A) from the factorization A = U*D*U'. */
00177 
00178 /*        K is the main loop index, increasing from 1 to N in steps of */
00179 /*        1 or 2, depending on the size of the diagonal blocks. */
00180 
00181         k = 1;
00182 L30:
00183 
00184 /*        If K > N, exit from loop. */
00185 
00186         if (k > *n) {
00187             goto L40;
00188         }
00189 
00190         if (ipiv[k] > 0) {
00191 
00192 /*           1 x 1 diagonal block */
00193 
00194 /*           Invert the diagonal block. */
00195 
00196             a[k + k * a_dim1] = 1. / a[k + k * a_dim1];
00197 
00198 /*           Compute column K of the inverse. */
00199 
00200             if (k > 1) {
00201                 i__1 = k - 1;
00202                 dcopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &work[1], &c__1);
00203                 i__1 = k - 1;
00204                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
00205                         c__1, &c_b13, &a[k * a_dim1 + 1], &c__1);
00206                 i__1 = k - 1;
00207                 a[k + k * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &a[k * 
00208                         a_dim1 + 1], &c__1);
00209             }
00210             kstep = 1;
00211         } else {
00212 
00213 /*           2 x 2 diagonal block */
00214 
00215 /*           Invert the diagonal block. */
00216 
00217             t = (d__1 = a[k + (k + 1) * a_dim1], abs(d__1));
00218             ak = a[k + k * a_dim1] / t;
00219             akp1 = a[k + 1 + (k + 1) * a_dim1] / t;
00220             akkp1 = a[k + (k + 1) * a_dim1] / t;
00221             d__ = t * (ak * akp1 - 1.);
00222             a[k + k * a_dim1] = akp1 / d__;
00223             a[k + 1 + (k + 1) * a_dim1] = ak / d__;
00224             a[k + (k + 1) * a_dim1] = -akkp1 / d__;
00225 
00226 /*           Compute columns K and K+1 of the inverse. */
00227 
00228             if (k > 1) {
00229                 i__1 = k - 1;
00230                 dcopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &work[1], &c__1);
00231                 i__1 = k - 1;
00232                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
00233                         c__1, &c_b13, &a[k * a_dim1 + 1], &c__1);
00234                 i__1 = k - 1;
00235                 a[k + k * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &a[k * 
00236                         a_dim1 + 1], &c__1);
00237                 i__1 = k - 1;
00238                 a[k + (k + 1) * a_dim1] -= ddot_(&i__1, &a[k * a_dim1 + 1], &
00239                         c__1, &a[(k + 1) * a_dim1 + 1], &c__1);
00240                 i__1 = k - 1;
00241                 dcopy_(&i__1, &a[(k + 1) * a_dim1 + 1], &c__1, &work[1], &
00242                         c__1);
00243                 i__1 = k - 1;
00244                 dsymv_(uplo, &i__1, &c_b11, &a[a_offset], lda, &work[1], &
00245                         c__1, &c_b13, &a[(k + 1) * a_dim1 + 1], &c__1);
00246                 i__1 = k - 1;
00247                 a[k + 1 + (k + 1) * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &
00248                         a[(k + 1) * a_dim1 + 1], &c__1);
00249             }
00250             kstep = 2;
00251         }
00252 
00253         kp = (i__1 = ipiv[k], abs(i__1));
00254         if (kp != k) {
00255 
00256 /*           Interchange rows and columns K and KP in the leading */
00257 /*           submatrix A(1:k+1,1:k+1) */
00258 
00259             i__1 = kp - 1;
00260             dswap_(&i__1, &a[k * a_dim1 + 1], &c__1, &a[kp * a_dim1 + 1], &
00261                     c__1);
00262             i__1 = k - kp - 1;
00263             dswap_(&i__1, &a[kp + 1 + k * a_dim1], &c__1, &a[kp + (kp + 1) * 
00264                     a_dim1], lda);
00265             temp = a[k + k * a_dim1];
00266             a[k + k * a_dim1] = a[kp + kp * a_dim1];
00267             a[kp + kp * a_dim1] = temp;
00268             if (kstep == 2) {
00269                 temp = a[k + (k + 1) * a_dim1];
00270                 a[k + (k + 1) * a_dim1] = a[kp + (k + 1) * a_dim1];
00271                 a[kp + (k + 1) * a_dim1] = temp;
00272             }
00273         }
00274 
00275         k += kstep;
00276         goto L30;
00277 L40:
00278 
00279         ;
00280     } else {
00281 
00282 /*        Compute inv(A) from the factorization A = L*D*L'. */
00283 
00284 /*        K is the main loop index, increasing from 1 to N in steps of */
00285 /*        1 or 2, depending on the size of the diagonal blocks. */
00286 
00287         k = *n;
00288 L50:
00289 
00290 /*        If K < 1, exit from loop. */
00291 
00292         if (k < 1) {
00293             goto L60;
00294         }
00295 
00296         if (ipiv[k] > 0) {
00297 
00298 /*           1 x 1 diagonal block */
00299 
00300 /*           Invert the diagonal block. */
00301 
00302             a[k + k * a_dim1] = 1. / a[k + k * a_dim1];
00303 
00304 /*           Compute column K of the inverse. */
00305 
00306             if (k < *n) {
00307                 i__1 = *n - k;
00308                 dcopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &work[1], &c__1);
00309                 i__1 = *n - k;
00310                 dsymv_(uplo, &i__1, &c_b11, &a[k + 1 + (k + 1) * a_dim1], lda, 
00311                          &work[1], &c__1, &c_b13, &a[k + 1 + k * a_dim1], &
00312                         c__1);
00313                 i__1 = *n - k;
00314                 a[k + k * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &a[k + 1 + 
00315                         k * a_dim1], &c__1);
00316             }
00317             kstep = 1;
00318         } else {
00319 
00320 /*           2 x 2 diagonal block */
00321 
00322 /*           Invert the diagonal block. */
00323 
00324             t = (d__1 = a[k + (k - 1) * a_dim1], abs(d__1));
00325             ak = a[k - 1 + (k - 1) * a_dim1] / t;
00326             akp1 = a[k + k * a_dim1] / t;
00327             akkp1 = a[k + (k - 1) * a_dim1] / t;
00328             d__ = t * (ak * akp1 - 1.);
00329             a[k - 1 + (k - 1) * a_dim1] = akp1 / d__;
00330             a[k + k * a_dim1] = ak / d__;
00331             a[k + (k - 1) * a_dim1] = -akkp1 / d__;
00332 
00333 /*           Compute columns K-1 and K of the inverse. */
00334 
00335             if (k < *n) {
00336                 i__1 = *n - k;
00337                 dcopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &work[1], &c__1);
00338                 i__1 = *n - k;
00339                 dsymv_(uplo, &i__1, &c_b11, &a[k + 1 + (k + 1) * a_dim1], lda, 
00340                          &work[1], &c__1, &c_b13, &a[k + 1 + k * a_dim1], &
00341                         c__1);
00342                 i__1 = *n - k;
00343                 a[k + k * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &a[k + 1 + 
00344                         k * a_dim1], &c__1);
00345                 i__1 = *n - k;
00346                 a[k + (k - 1) * a_dim1] -= ddot_(&i__1, &a[k + 1 + k * a_dim1]
00347 , &c__1, &a[k + 1 + (k - 1) * a_dim1], &c__1);
00348                 i__1 = *n - k;
00349                 dcopy_(&i__1, &a[k + 1 + (k - 1) * a_dim1], &c__1, &work[1], &
00350                         c__1);
00351                 i__1 = *n - k;
00352                 dsymv_(uplo, &i__1, &c_b11, &a[k + 1 + (k + 1) * a_dim1], lda, 
00353                          &work[1], &c__1, &c_b13, &a[k + 1 + (k - 1) * a_dim1]
00354 , &c__1);
00355                 i__1 = *n - k;
00356                 a[k - 1 + (k - 1) * a_dim1] -= ddot_(&i__1, &work[1], &c__1, &
00357                         a[k + 1 + (k - 1) * a_dim1], &c__1);
00358             }
00359             kstep = 2;
00360         }
00361 
00362         kp = (i__1 = ipiv[k], abs(i__1));
00363         if (kp != k) {
00364 
00365 /*           Interchange rows and columns K and KP in the trailing */
00366 /*           submatrix A(k-1:n,k-1:n) */
00367 
00368             if (kp < *n) {
00369                 i__1 = *n - kp;
00370                 dswap_(&i__1, &a[kp + 1 + k * a_dim1], &c__1, &a[kp + 1 + kp *
00371                          a_dim1], &c__1);
00372             }
00373             i__1 = kp - k - 1;
00374             dswap_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &a[kp + (k + 1) * 
00375                     a_dim1], lda);
00376             temp = a[k + k * a_dim1];
00377             a[k + k * a_dim1] = a[kp + kp * a_dim1];
00378             a[kp + kp * a_dim1] = temp;
00379             if (kstep == 2) {
00380                 temp = a[k + (k - 1) * a_dim1];
00381                 a[k + (k - 1) * a_dim1] = a[kp + (k - 1) * a_dim1];
00382                 a[kp + (k - 1) * a_dim1] = temp;
00383             }
00384         }
00385 
00386         k -= kstep;
00387         goto L50;
00388 L60:
00389         ;
00390     }
00391 
00392     return 0;
00393 
00394 /*     End of DSYTRI */
00395 
00396 } /* dsytri_ */