swiftnav: chegs2.c Source File

Go to the documentation of this file.
00001 /* chegs2.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 integer c__1 = 1;
00020 
00021 /* Subroutine */ int chegs2_(integer *itype, char *uplo, integer *n, complex *
00022         a, integer *lda, complex *b, integer *ldb, integer *info)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
00026     real r__1, r__2;
00027     complex q__1;
00028 
00029     /* Local variables */
00030     integer k;
00031     complex ct;
00032     real akk, bkk;
00033     extern /* Subroutine */ int cher2_(char *, integer *, complex *, complex *
00034 , integer *, complex *, integer *, complex *, integer *);
00035     extern logical lsame_(char *, char *);
00036     extern /* Subroutine */ int caxpy_(integer *, complex *, complex *, 
00037             integer *, complex *, integer *);
00038     logical upper;
00039     extern /* Subroutine */ int ctrmv_(char *, char *, char *, integer *, 
00040             complex *, integer *, complex *, integer *), ctrsv_(char *, char *, char *, integer *, complex *, 
00041             integer *, complex *, integer *), clacgv_(
00042             integer *, complex *, integer *), csscal_(integer *, real *, 
00043             complex *, integer *), xerbla_(char *, integer *);
00044 
00045 
00046 /*  -- LAPACK routine (version 3.2) -- */
00047 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00048 /*     November 2006 */
00049 
00050 /*     .. Scalar Arguments .. */
00051 /*     .. */
00052 /*     .. Array Arguments .. */
00053 /*     .. */
00054 
00055 /*  Purpose */
00056 /*  ======= */
00057 
00058 /*  CHEGS2 reduces a complex Hermitian-definite generalized */
00059 /*  eigenproblem to standard form. */
00060 
00061 /*  If ITYPE = 1, the problem is A*x = lambda*B*x, */
00062 /*  and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L') */
00063 
00064 /*  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or */
00065 /*  B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L. */
00066 
00067 /*  B must have been previously factorized as U'*U or L*L' by CPOTRF. */
00068 
00069 /*  Arguments */
00070 /*  ========= */
00071 
00072 /*  ITYPE   (input) INTEGER */
00073 /*          = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L'); */
00074 /*          = 2 or 3: compute U*A*U' or L'*A*L. */
00075 
00076 /*  UPLO    (input) CHARACTER*1 */
00077 /*          Specifies whether the upper or lower triangular part of the */
00078 /*          Hermitian matrix A is stored, and how B has been factorized. */
00079 /*          = 'U':  Upper triangular */
00080 /*          = 'L':  Lower triangular */
00081 
00082 /*  N       (input) INTEGER */
00083 /*          The order of the matrices A and B.  N >= 0. */
00084 
00085 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00086 /*          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading */
00087 /*          n by n upper triangular part of A contains the upper */
00088 /*          triangular part of the matrix A, and the strictly lower */
00089 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
00090 /*          leading n by n lower triangular part of A contains the lower */
00091 /*          triangular part of the matrix A, and the strictly upper */
00092 /*          triangular part of A is not referenced. */
00093 
00094 /*          On exit, if INFO = 0, the transformed matrix, stored in the */
00095 /*          same format as A. */
00096 
00097 /*  LDA     (input) INTEGER */
00098 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00099 
00100 /*  B       (input) COMPLEX array, dimension (LDB,N) */
00101 /*          The triangular factor from the Cholesky factorization of B, */
00102 /*          as returned by CPOTRF. */
00103 
00104 /*  LDB     (input) INTEGER */
00105 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00106 
00107 /*  INFO    (output) INTEGER */
00108 /*          = 0:  successful exit. */
00109 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00110 
00111 /*  ===================================================================== */
00112 
00113 /*     .. Parameters .. */
00114 /*     .. */
00115 /*     .. Local Scalars .. */
00116 /*     .. */
00117 /*     .. External Subroutines .. */
00118 /*     .. */
00119 /*     .. Intrinsic Functions .. */
00120 /*     .. */
00121 /*     .. External Functions .. */
00122 /*     .. */
00123 /*     .. Executable Statements .. */
00124 
00125 /*     Test the input parameters. */
00126 
00127     /* Parameter adjustments */
00128     a_dim1 = *lda;
00129     a_offset = 1 + a_dim1;
00130     a -= a_offset;
00131     b_dim1 = *ldb;
00132     b_offset = 1 + b_dim1;
00133     b -= b_offset;
00134 
00135     /* Function Body */
00136     *info = 0;
00137     upper = lsame_(uplo, "U");
00138     if (*itype < 1 || *itype > 3) {
00139         *info = -1;
00140     } else if (! upper && ! lsame_(uplo, "L")) {
00141         *info = -2;
00142     } else if (*n < 0) {
00143         *info = -3;
00144     } else if (*lda < max(1,*n)) {
00145         *info = -5;
00146     } else if (*ldb < max(1,*n)) {
00147         *info = -7;
00148     }
00149     if (*info != 0) {
00150         i__1 = -(*info);
00151         xerbla_("CHEGS2", &i__1);
00152         return 0;
00153     }
00154 
00155     if (*itype == 1) {
00156         if (upper) {
00157 
00158 /*           Compute inv(U')*A*inv(U) */
00159 
00160             i__1 = *n;
00161             for (k = 1; k <= i__1; ++k) {
00162 
00163 /*              Update the upper triangle of A(k:n,k:n) */
00164 
00165                 i__2 = k + k * a_dim1;
00166                 akk = a[i__2].r;
00167                 i__2 = k + k * b_dim1;
00168                 bkk = b[i__2].r;
00169 /* Computing 2nd power */
00170                 r__1 = bkk;
00171                 akk /= r__1 * r__1;
00172                 i__2 = k + k * a_dim1;
00173                 a[i__2].r = akk, a[i__2].i = 0.f;
00174                 if (k < *n) {
00175                     i__2 = *n - k;
00176                     r__1 = 1.f / bkk;
00177                     csscal_(&i__2, &r__1, &a[k + (k + 1) * a_dim1], lda);
00178                     r__1 = akk * -.5f;
00179                     ct.r = r__1, ct.i = 0.f;
00180                     i__2 = *n - k;
00181                     clacgv_(&i__2, &a[k + (k + 1) * a_dim1], lda);
00182                     i__2 = *n - k;
00183                     clacgv_(&i__2, &b[k + (k + 1) * b_dim1], ldb);
00184                     i__2 = *n - k;
00185                     caxpy_(&i__2, &ct, &b[k + (k + 1) * b_dim1], ldb, &a[k + (
00186                             k + 1) * a_dim1], lda);
00187                     i__2 = *n - k;
00188                     q__1.r = -1.f, q__1.i = -0.f;
00189                     cher2_(uplo, &i__2, &q__1, &a[k + (k + 1) * a_dim1], lda, 
00190                             &b[k + (k + 1) * b_dim1], ldb, &a[k + 1 + (k + 1) 
00191                             * a_dim1], lda);
00192                     i__2 = *n - k;
00193                     caxpy_(&i__2, &ct, &b[k + (k + 1) * b_dim1], ldb, &a[k + (
00194                             k + 1) * a_dim1], lda);
00195                     i__2 = *n - k;
00196                     clacgv_(&i__2, &b[k + (k + 1) * b_dim1], ldb);
00197                     i__2 = *n - k;
00198                     ctrsv_(uplo, "Conjugate transpose", "Non-unit", &i__2, &b[
00199                             k + 1 + (k + 1) * b_dim1], ldb, &a[k + (k + 1) * 
00200                             a_dim1], lda);
00201                     i__2 = *n - k;
00202                     clacgv_(&i__2, &a[k + (k + 1) * a_dim1], lda);
00203                 }
00204 /* L10: */
00205             }
00206         } else {
00207 
00208 /*           Compute inv(L)*A*inv(L') */
00209 
00210             i__1 = *n;
00211             for (k = 1; k <= i__1; ++k) {
00212 
00213 /*              Update the lower triangle of A(k:n,k:n) */
00214 
00215                 i__2 = k + k * a_dim1;
00216                 akk = a[i__2].r;
00217                 i__2 = k + k * b_dim1;
00218                 bkk = b[i__2].r;
00219 /* Computing 2nd power */
00220                 r__1 = bkk;
00221                 akk /= r__1 * r__1;
00222                 i__2 = k + k * a_dim1;
00223                 a[i__2].r = akk, a[i__2].i = 0.f;
00224                 if (k < *n) {
00225                     i__2 = *n - k;
00226                     r__1 = 1.f / bkk;
00227                     csscal_(&i__2, &r__1, &a[k + 1 + k * a_dim1], &c__1);
00228                     r__1 = akk * -.5f;
00229                     ct.r = r__1, ct.i = 0.f;
00230                     i__2 = *n - k;
00231                     caxpy_(&i__2, &ct, &b[k + 1 + k * b_dim1], &c__1, &a[k + 
00232                             1 + k * a_dim1], &c__1);
00233                     i__2 = *n - k;
00234                     q__1.r = -1.f, q__1.i = -0.f;
00235                     cher2_(uplo, &i__2, &q__1, &a[k + 1 + k * a_dim1], &c__1, 
00236                             &b[k + 1 + k * b_dim1], &c__1, &a[k + 1 + (k + 1) 
00237                             * a_dim1], lda);
00238                     i__2 = *n - k;
00239                     caxpy_(&i__2, &ct, &b[k + 1 + k * b_dim1], &c__1, &a[k + 
00240                             1 + k * a_dim1], &c__1);
00241                     i__2 = *n - k;
00242                     ctrsv_(uplo, "No transpose", "Non-unit", &i__2, &b[k + 1 
00243                             + (k + 1) * b_dim1], ldb, &a[k + 1 + k * a_dim1], 
00244                             &c__1);
00245                 }
00246 /* L20: */
00247             }
00248         }
00249     } else {
00250         if (upper) {
00251 
00252 /*           Compute U*A*U' */
00253 
00254             i__1 = *n;
00255             for (k = 1; k <= i__1; ++k) {
00256 
00257 /*              Update the upper triangle of A(1:k,1:k) */
00258 
00259                 i__2 = k + k * a_dim1;
00260                 akk = a[i__2].r;
00261                 i__2 = k + k * b_dim1;
00262                 bkk = b[i__2].r;
00263                 i__2 = k - 1;
00264                 ctrmv_(uplo, "No transpose", "Non-unit", &i__2, &b[b_offset], 
00265                         ldb, &a[k * a_dim1 + 1], &c__1);
00266                 r__1 = akk * .5f;
00267                 ct.r = r__1, ct.i = 0.f;
00268                 i__2 = k - 1;
00269                 caxpy_(&i__2, &ct, &b[k * b_dim1 + 1], &c__1, &a[k * a_dim1 + 
00270                         1], &c__1);
00271                 i__2 = k - 1;
00272                 cher2_(uplo, &i__2, &c_b1, &a[k * a_dim1 + 1], &c__1, &b[k * 
00273                         b_dim1 + 1], &c__1, &a[a_offset], lda);
00274                 i__2 = k - 1;
00275                 caxpy_(&i__2, &ct, &b[k * b_dim1 + 1], &c__1, &a[k * a_dim1 + 
00276                         1], &c__1);
00277                 i__2 = k - 1;
00278                 csscal_(&i__2, &bkk, &a[k * a_dim1 + 1], &c__1);
00279                 i__2 = k + k * a_dim1;
00280 /* Computing 2nd power */
00281                 r__2 = bkk;
00282                 r__1 = akk * (r__2 * r__2);
00283                 a[i__2].r = r__1, a[i__2].i = 0.f;
00284 /* L30: */
00285             }
00286         } else {
00287 
00288 /*           Compute L'*A*L */
00289 
00290             i__1 = *n;
00291             for (k = 1; k <= i__1; ++k) {
00292 
00293 /*              Update the lower triangle of A(1:k,1:k) */
00294 
00295                 i__2 = k + k * a_dim1;
00296                 akk = a[i__2].r;
00297                 i__2 = k + k * b_dim1;
00298                 bkk = b[i__2].r;
00299                 i__2 = k - 1;
00300                 clacgv_(&i__2, &a[k + a_dim1], lda);
00301                 i__2 = k - 1;
00302                 ctrmv_(uplo, "Conjugate transpose", "Non-unit", &i__2, &b[
00303                         b_offset], ldb, &a[k + a_dim1], lda);
00304                 r__1 = akk * .5f;
00305                 ct.r = r__1, ct.i = 0.f;
00306                 i__2 = k - 1;
00307                 clacgv_(&i__2, &b[k + b_dim1], ldb);
00308                 i__2 = k - 1;
00309                 caxpy_(&i__2, &ct, &b[k + b_dim1], ldb, &a[k + a_dim1], lda);
00310                 i__2 = k - 1;
00311                 cher2_(uplo, &i__2, &c_b1, &a[k + a_dim1], lda, &b[k + b_dim1]
00312 , ldb, &a[a_offset], lda);
00313                 i__2 = k - 1;
00314                 caxpy_(&i__2, &ct, &b[k + b_dim1], ldb, &a[k + a_dim1], lda);
00315                 i__2 = k - 1;
00316                 clacgv_(&i__2, &b[k + b_dim1], ldb);
00317                 i__2 = k - 1;
00318                 csscal_(&i__2, &bkk, &a[k + a_dim1], lda);
00319                 i__2 = k - 1;
00320                 clacgv_(&i__2, &a[k + a_dim1], lda);
00321                 i__2 = k + k * a_dim1;
00322 /* Computing 2nd power */
00323                 r__2 = bkk;
00324                 r__1 = akk * (r__2 * r__2);
00325                 a[i__2].r = r__1, a[i__2].i = 0.f;
00326 /* L40: */
00327             }
00328         }
00329     }
00330     return 0;
00331 
00332 /*     End of CHEGS2 */
00333 
00334 } /* chegs2_ */