swiftnav: ssptri.c Source File

Go to the documentation of this file.
00001 /* ssptri.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 ssptri_(char *uplo, integer *n, real *ap, integer *ipiv, 
00023         real *work, integer *info)
00024 {
00025     /* System generated locals */
00026     integer i__1;
00027     real r__1;
00028 
00029     /* Local variables */
00030     real d__;
00031     integer j, k;
00032     real t, ak;
00033     integer kc, kp, kx, kpc, npp;
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 ), sspmv_(char *, integer *, real *, real *, real *, integer *, 
00043             real *, real *, integer *), xerbla_(char *, integer *);
00044     integer kcnext;
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 /*  SSPTRI computes the inverse of a real symmetric indefinite matrix */
00060 /*  A in packed storage using the factorization A = U*D*U**T or */
00061 /*  A = L*D*L**T computed by SSPTRF. */
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 /*  AP      (input/output) REAL array, dimension (N*(N+1)/2) */
00076 /*          On entry, the block diagonal matrix D and the multipliers */
00077 /*          used to obtain the factor U or L as computed by SSPTRF, */
00078 /*          stored as a packed triangular matrix. */
00079 
00080 /*          On exit, if INFO = 0, the (symmetric) inverse of the original */
00081 /*          matrix, stored as a packed triangular matrix. The j-th column */
00082 /*          of inv(A) is stored in the array AP as follows: */
00083 /*          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j; */
00084 /*          if UPLO = 'L', */
00085 /*             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n. */
00086 
00087 /*  IPIV    (input) INTEGER array, dimension (N) */
00088 /*          Details of the interchanges and the block structure of D */
00089 /*          as determined by SSPTRF. */
00090 
00091 /*  WORK    (workspace) REAL array, dimension (N) */
00092 
00093 /*  INFO    (output) INTEGER */
00094 /*          = 0: successful exit */
00095 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00096 /*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its */
00097 /*               inverse could not be computed. */
00098 
00099 /*  ===================================================================== */
00100 
00101 /*     .. Parameters .. */
00102 /*     .. */
00103 /*     .. Local Scalars .. */
00104 /*     .. */
00105 /*     .. External Functions .. */
00106 /*     .. */
00107 /*     .. External Subroutines .. */
00108 /*     .. */
00109 /*     .. Intrinsic Functions .. */
00110 /*     .. */
00111 /*     .. Executable Statements .. */
00112 
00113 /*     Test the input parameters. */
00114 
00115     /* Parameter adjustments */
00116     --work;
00117     --ipiv;
00118     --ap;
00119 
00120     /* Function Body */
00121     *info = 0;
00122     upper = lsame_(uplo, "U");
00123     if (! upper && ! lsame_(uplo, "L")) {
00124         *info = -1;
00125     } else if (*n < 0) {
00126         *info = -2;
00127     }
00128     if (*info != 0) {
00129         i__1 = -(*info);
00130         xerbla_("SSPTRI", &i__1);
00131         return 0;
00132     }
00133 
00134 /*     Quick return if possible */
00135 
00136     if (*n == 0) {
00137         return 0;
00138     }
00139 
00140 /*     Check that the diagonal matrix D is nonsingular. */
00141 
00142     if (upper) {
00143 
00144 /*        Upper triangular storage: examine D from bottom to top */
00145 
00146         kp = *n * (*n + 1) / 2;
00147         for (*info = *n; *info >= 1; --(*info)) {
00148             if (ipiv[*info] > 0 && ap[kp] == 0.f) {
00149                 return 0;
00150             }
00151             kp -= *info;
00152 /* L10: */
00153         }
00154     } else {
00155 
00156 /*        Lower triangular storage: examine D from top to bottom. */
00157 
00158         kp = 1;
00159         i__1 = *n;
00160         for (*info = 1; *info <= i__1; ++(*info)) {
00161             if (ipiv[*info] > 0 && ap[kp] == 0.f) {
00162                 return 0;
00163             }
00164             kp = kp + *n - *info + 1;
00165 /* L20: */
00166         }
00167     }
00168     *info = 0;
00169 
00170     if (upper) {
00171 
00172 /*        Compute inv(A) from the factorization A = U*D*U'. */
00173 
00174 /*        K is the main loop index, increasing from 1 to N in steps of */
00175 /*        1 or 2, depending on the size of the diagonal blocks. */
00176 
00177         k = 1;
00178         kc = 1;
00179 L30:
00180 
00181 /*        If K > N, exit from loop. */
00182 
00183         if (k > *n) {
00184             goto L50;
00185         }
00186 
00187         kcnext = kc + k;
00188         if (ipiv[k] > 0) {
00189 
00190 /*           1 x 1 diagonal block */
00191 
00192 /*           Invert the diagonal block. */
00193 
00194             ap[kc + k - 1] = 1.f / ap[kc + k - 1];
00195 
00196 /*           Compute column K of the inverse. */
00197 
00198             if (k > 1) {
00199                 i__1 = k - 1;
00200                 scopy_(&i__1, &ap[kc], &c__1, &work[1], &c__1);
00201                 i__1 = k - 1;
00202                 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, &
00203                         ap[kc], &c__1);
00204                 i__1 = k - 1;
00205                 ap[kc + k - 1] -= sdot_(&i__1, &work[1], &c__1, &ap[kc], &
00206                         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 = ap[kcnext + k - 1], dabs(r__1));
00216             ak = ap[kc + k - 1] / t;
00217             akp1 = ap[kcnext + k] / t;
00218             akkp1 = ap[kcnext + k - 1] / t;
00219             d__ = t * (ak * akp1 - 1.f);
00220             ap[kc + k - 1] = akp1 / d__;
00221             ap[kcnext + k] = ak / d__;
00222             ap[kcnext + k - 1] = -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, &ap[kc], &c__1, &work[1], &c__1);
00229                 i__1 = k - 1;
00230                 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, &
00231                         ap[kc], &c__1);
00232                 i__1 = k - 1;
00233                 ap[kc + k - 1] -= sdot_(&i__1, &work[1], &c__1, &ap[kc], &
00234                         c__1);
00235                 i__1 = k - 1;
00236                 ap[kcnext + k - 1] -= sdot_(&i__1, &ap[kc], &c__1, &ap[kcnext]
00237 , &c__1);
00238                 i__1 = k - 1;
00239                 scopy_(&i__1, &ap[kcnext], &c__1, &work[1], &c__1);
00240                 i__1 = k - 1;
00241                 sspmv_(uplo, &i__1, &c_b11, &ap[1], &work[1], &c__1, &c_b13, &
00242                         ap[kcnext], &c__1);
00243                 i__1 = k - 1;
00244                 ap[kcnext + k] -= sdot_(&i__1, &work[1], &c__1, &ap[kcnext], &
00245                         c__1);
00246             }
00247             kstep = 2;
00248             kcnext = kcnext + k + 1;
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             kpc = (kp - 1) * kp / 2 + 1;
00258             i__1 = kp - 1;
00259             sswap_(&i__1, &ap[kc], &c__1, &ap[kpc], &c__1);
00260             kx = kpc + kp - 1;
00261             i__1 = k - 1;
00262             for (j = kp + 1; j <= i__1; ++j) {
00263                 kx = kx + j - 1;
00264                 temp = ap[kc + j - 1];
00265                 ap[kc + j - 1] = ap[kx];
00266                 ap[kx] = temp;
00267 /* L40: */
00268             }
00269             temp = ap[kc + k - 1];
00270             ap[kc + k - 1] = ap[kpc + kp - 1];
00271             ap[kpc + kp - 1] = temp;
00272             if (kstep == 2) {
00273                 temp = ap[kc + k + k - 1];
00274                 ap[kc + k + k - 1] = ap[kc + k + kp - 1];
00275                 ap[kc + k + kp - 1] = temp;
00276             }
00277         }
00278 
00279         k += kstep;
00280         kc = kcnext;
00281         goto L30;
00282 L50:
00283 
00284         ;
00285     } else {
00286 
00287 /*        Compute inv(A) from the factorization A = L*D*L'. */
00288 
00289 /*        K is the main loop index, increasing from 1 to N in steps of */
00290 /*        1 or 2, depending on the size of the diagonal blocks. */
00291 
00292         npp = *n * (*n + 1) / 2;
00293         k = *n;
00294         kc = npp;
00295 L60:
00296 
00297 /*        If K < 1, exit from loop. */
00298 
00299         if (k < 1) {
00300             goto L80;
00301         }
00302 
00303         kcnext = kc - (*n - k + 2);
00304         if (ipiv[k] > 0) {
00305 
00306 /*           1 x 1 diagonal block */
00307 
00308 /*           Invert the diagonal block. */
00309 
00310             ap[kc] = 1.f / ap[kc];
00311 
00312 /*           Compute column K of the inverse. */
00313 
00314             if (k < *n) {
00315                 i__1 = *n - k;
00316                 scopy_(&i__1, &ap[kc + 1], &c__1, &work[1], &c__1);
00317                 i__1 = *n - k;
00318                 sspmv_(uplo, &i__1, &c_b11, &ap[kc + *n - k + 1], &work[1], &
00319                         c__1, &c_b13, &ap[kc + 1], &c__1);
00320                 i__1 = *n - k;
00321                 ap[kc] -= sdot_(&i__1, &work[1], &c__1, &ap[kc + 1], &c__1);
00322             }
00323             kstep = 1;
00324         } else {
00325 
00326 /*           2 x 2 diagonal block */
00327 
00328 /*           Invert the diagonal block. */
00329 
00330             t = (r__1 = ap[kcnext + 1], dabs(r__1));
00331             ak = ap[kcnext] / t;
00332             akp1 = ap[kc] / t;
00333             akkp1 = ap[kcnext + 1] / t;
00334             d__ = t * (ak * akp1 - 1.f);
00335             ap[kcnext] = akp1 / d__;
00336             ap[kc] = ak / d__;
00337             ap[kcnext + 1] = -akkp1 / d__;
00338 
00339 /*           Compute columns K-1 and K of the inverse. */
00340 
00341             if (k < *n) {
00342                 i__1 = *n - k;
00343                 scopy_(&i__1, &ap[kc + 1], &c__1, &work[1], &c__1);
00344                 i__1 = *n - k;
00345                 sspmv_(uplo, &i__1, &c_b11, &ap[kc + (*n - k + 1)], &work[1], 
00346                         &c__1, &c_b13, &ap[kc + 1], &c__1);
00347                 i__1 = *n - k;
00348                 ap[kc] -= sdot_(&i__1, &work[1], &c__1, &ap[kc + 1], &c__1);
00349                 i__1 = *n - k;
00350                 ap[kcnext + 1] -= sdot_(&i__1, &ap[kc + 1], &c__1, &ap[kcnext 
00351                         + 2], &c__1);
00352                 i__1 = *n - k;
00353                 scopy_(&i__1, &ap[kcnext + 2], &c__1, &work[1], &c__1);
00354                 i__1 = *n - k;
00355                 sspmv_(uplo, &i__1, &c_b11, &ap[kc + (*n - k + 1)], &work[1], 
00356                         &c__1, &c_b13, &ap[kcnext + 2], &c__1);
00357                 i__1 = *n - k;
00358                 ap[kcnext] -= sdot_(&i__1, &work[1], &c__1, &ap[kcnext + 2], &
00359                         c__1);
00360             }
00361             kstep = 2;
00362             kcnext -= *n - k + 3;
00363         }
00364 
00365         kp = (i__1 = ipiv[k], abs(i__1));
00366         if (kp != k) {
00367 
00368 /*           Interchange rows and columns K and KP in the trailing */
00369 /*           submatrix A(k-1:n,k-1:n) */
00370 
00371             kpc = npp - (*n - kp + 1) * (*n - kp + 2) / 2 + 1;
00372             if (kp < *n) {
00373                 i__1 = *n - kp;
00374                 sswap_(&i__1, &ap[kc + kp - k + 1], &c__1, &ap[kpc + 1], &
00375                         c__1);
00376             }
00377             kx = kc + kp - k;
00378             i__1 = kp - 1;
00379             for (j = k + 1; j <= i__1; ++j) {
00380                 kx = kx + *n - j + 1;
00381                 temp = ap[kc + j - k];
00382                 ap[kc + j - k] = ap[kx];
00383                 ap[kx] = temp;
00384 /* L70: */
00385             }
00386             temp = ap[kc];
00387             ap[kc] = ap[kpc];
00388             ap[kpc] = temp;
00389             if (kstep == 2) {
00390                 temp = ap[kc - *n + k - 1];
00391                 ap[kc - *n + k - 1] = ap[kc - *n + kp - 1];
00392                 ap[kc - *n + kp - 1] = temp;
00393             }
00394         }
00395 
00396         k -= kstep;
00397         kc = kcnext;
00398         goto L60;
00399 L80:
00400         ;
00401     }
00402 
00403     return 0;
00404 
00405 /*     End of SSPTRI */
00406 
00407 } /* ssptri_ */