swiftnav: ssytri.c Source File

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