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

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:34