zgbtrs.c
Go to the documentation of this file.
00001 /* zgbtrs.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 doublecomplex c_b1 = {1.,0.};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int zgbtrs_(char *trans, integer *n, integer *kl, integer *
00022         ku, integer *nrhs, doublecomplex *ab, integer *ldab, integer *ipiv, 
00023         doublecomplex *b, integer *ldb, integer *info)
00024 {
00025     /* System generated locals */
00026     integer ab_dim1, ab_offset, b_dim1, b_offset, i__1, i__2, i__3;
00027     doublecomplex z__1;
00028 
00029     /* Local variables */
00030     integer i__, j, l, kd, lm;
00031     extern logical lsame_(char *, char *);
00032     logical lnoti;
00033     extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
00034             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00035             integer *, doublecomplex *, doublecomplex *, integer *), 
00036             zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, 
00037             integer *, doublecomplex *, integer *, doublecomplex *, integer *)
00038             , zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 
00039             integer *), ztbsv_(char *, char *, char *, integer *, integer *, 
00040             doublecomplex *, integer *, doublecomplex *, integer *), xerbla_(char *, integer *), zlacgv_(
00041             integer *, doublecomplex *, integer *);
00042     logical notran;
00043 
00044 
00045 /*  -- LAPACK routine (version 3.2) -- */
00046 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00047 /*     November 2006 */
00048 
00049 /*     .. Scalar Arguments .. */
00050 /*     .. */
00051 /*     .. Array Arguments .. */
00052 /*     .. */
00053 
00054 /*  Purpose */
00055 /*  ======= */
00056 
00057 /*  ZGBTRS solves a system of linear equations */
00058 /*     A * X = B,  A**T * X = B,  or  A**H * X = B */
00059 /*  with a general band matrix A using the LU factorization computed */
00060 /*  by ZGBTRF. */
00061 
00062 /*  Arguments */
00063 /*  ========= */
00064 
00065 /*  TRANS   (input) CHARACTER*1 */
00066 /*          Specifies the form of the system of equations. */
00067 /*          = 'N':  A * X = B     (No transpose) */
00068 /*          = 'T':  A**T * X = B  (Transpose) */
00069 /*          = 'C':  A**H * X = B  (Conjugate transpose) */
00070 
00071 /*  N       (input) INTEGER */
00072 /*          The order of the matrix A.  N >= 0. */
00073 
00074 /*  KL      (input) INTEGER */
00075 /*          The number of subdiagonals within the band of A.  KL >= 0. */
00076 
00077 /*  KU      (input) INTEGER */
00078 /*          The number of superdiagonals within the band of A.  KU >= 0. */
00079 
00080 /*  NRHS    (input) INTEGER */
00081 /*          The number of right hand sides, i.e., the number of columns */
00082 /*          of the matrix B.  NRHS >= 0. */
00083 
00084 /*  AB      (input) COMPLEX*16 array, dimension (LDAB,N) */
00085 /*          Details of the LU factorization of the band matrix A, as */
00086 /*          computed by ZGBTRF.  U is stored as an upper triangular band */
00087 /*          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and */
00088 /*          the multipliers used during the factorization are stored in */
00089 /*          rows KL+KU+2 to 2*KL+KU+1. */
00090 
00091 /*  LDAB    (input) INTEGER */
00092 /*          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1. */
00093 
00094 /*  IPIV    (input) INTEGER array, dimension (N) */
00095 /*          The pivot indices; for 1 <= i <= N, row i of the matrix was */
00096 /*          interchanged with row IPIV(i). */
00097 
00098 /*  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS) */
00099 /*          On entry, the right hand side matrix B. */
00100 /*          On exit, the solution matrix X. */
00101 
00102 /*  LDB     (input) INTEGER */
00103 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00104 
00105 /*  INFO    (output) INTEGER */
00106 /*          = 0:  successful exit */
00107 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00108 
00109 /*  ===================================================================== */
00110 
00111 /*     .. Parameters .. */
00112 /*     .. */
00113 /*     .. Local Scalars .. */
00114 /*     .. */
00115 /*     .. External Functions .. */
00116 /*     .. */
00117 /*     .. External Subroutines .. */
00118 /*     .. */
00119 /*     .. Intrinsic Functions .. */
00120 /*     .. */
00121 /*     .. Executable Statements .. */
00122 
00123 /*     Test the input parameters. */
00124 
00125     /* Parameter adjustments */
00126     ab_dim1 = *ldab;
00127     ab_offset = 1 + ab_dim1;
00128     ab -= ab_offset;
00129     --ipiv;
00130     b_dim1 = *ldb;
00131     b_offset = 1 + b_dim1;
00132     b -= b_offset;
00133 
00134     /* Function Body */
00135     *info = 0;
00136     notran = lsame_(trans, "N");
00137     if (! notran && ! lsame_(trans, "T") && ! lsame_(
00138             trans, "C")) {
00139         *info = -1;
00140     } else if (*n < 0) {
00141         *info = -2;
00142     } else if (*kl < 0) {
00143         *info = -3;
00144     } else if (*ku < 0) {
00145         *info = -4;
00146     } else if (*nrhs < 0) {
00147         *info = -5;
00148     } else if (*ldab < (*kl << 1) + *ku + 1) {
00149         *info = -7;
00150     } else if (*ldb < max(1,*n)) {
00151         *info = -10;
00152     }
00153     if (*info != 0) {
00154         i__1 = -(*info);
00155         xerbla_("ZGBTRS", &i__1);
00156         return 0;
00157     }
00158 
00159 /*     Quick return if possible */
00160 
00161     if (*n == 0 || *nrhs == 0) {
00162         return 0;
00163     }
00164 
00165     kd = *ku + *kl + 1;
00166     lnoti = *kl > 0;
00167 
00168     if (notran) {
00169 
00170 /*        Solve  A*X = B. */
00171 
00172 /*        Solve L*X = B, overwriting B with X. */
00173 
00174 /*        L is represented as a product of permutations and unit lower */
00175 /*        triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1), */
00176 /*        where each transformation L(i) is a rank-one modification of */
00177 /*        the identity matrix. */
00178 
00179         if (lnoti) {
00180             i__1 = *n - 1;
00181             for (j = 1; j <= i__1; ++j) {
00182 /* Computing MIN */
00183                 i__2 = *kl, i__3 = *n - j;
00184                 lm = min(i__2,i__3);
00185                 l = ipiv[j];
00186                 if (l != j) {
00187                     zswap_(nrhs, &b[l + b_dim1], ldb, &b[j + b_dim1], ldb);
00188                 }
00189                 z__1.r = -1., z__1.i = -0.;
00190                 zgeru_(&lm, nrhs, &z__1, &ab[kd + 1 + j * ab_dim1], &c__1, &b[
00191                         j + b_dim1], ldb, &b[j + 1 + b_dim1], ldb);
00192 /* L10: */
00193             }
00194         }
00195 
00196         i__1 = *nrhs;
00197         for (i__ = 1; i__ <= i__1; ++i__) {
00198 
00199 /*           Solve U*X = B, overwriting B with X. */
00200 
00201             i__2 = *kl + *ku;
00202             ztbsv_("Upper", "No transpose", "Non-unit", n, &i__2, &ab[
00203                     ab_offset], ldab, &b[i__ * b_dim1 + 1], &c__1);
00204 /* L20: */
00205         }
00206 
00207     } else if (lsame_(trans, "T")) {
00208 
00209 /*        Solve A**T * X = B. */
00210 
00211         i__1 = *nrhs;
00212         for (i__ = 1; i__ <= i__1; ++i__) {
00213 
00214 /*           Solve U**T * X = B, overwriting B with X. */
00215 
00216             i__2 = *kl + *ku;
00217             ztbsv_("Upper", "Transpose", "Non-unit", n, &i__2, &ab[ab_offset], 
00218                      ldab, &b[i__ * b_dim1 + 1], &c__1);
00219 /* L30: */
00220         }
00221 
00222 /*        Solve L**T * X = B, overwriting B with X. */
00223 
00224         if (lnoti) {
00225             for (j = *n - 1; j >= 1; --j) {
00226 /* Computing MIN */
00227                 i__1 = *kl, i__2 = *n - j;
00228                 lm = min(i__1,i__2);
00229                 z__1.r = -1., z__1.i = -0.;
00230                 zgemv_("Transpose", &lm, nrhs, &z__1, &b[j + 1 + b_dim1], ldb, 
00231                          &ab[kd + 1 + j * ab_dim1], &c__1, &c_b1, &b[j + 
00232                         b_dim1], ldb);
00233                 l = ipiv[j];
00234                 if (l != j) {
00235                     zswap_(nrhs, &b[l + b_dim1], ldb, &b[j + b_dim1], ldb);
00236                 }
00237 /* L40: */
00238             }
00239         }
00240 
00241     } else {
00242 
00243 /*        Solve A**H * X = B. */
00244 
00245         i__1 = *nrhs;
00246         for (i__ = 1; i__ <= i__1; ++i__) {
00247 
00248 /*           Solve U**H * X = B, overwriting B with X. */
00249 
00250             i__2 = *kl + *ku;
00251             ztbsv_("Upper", "Conjugate transpose", "Non-unit", n, &i__2, &ab[
00252                     ab_offset], ldab, &b[i__ * b_dim1 + 1], &c__1);
00253 /* L50: */
00254         }
00255 
00256 /*        Solve L**H * X = B, overwriting B with X. */
00257 
00258         if (lnoti) {
00259             for (j = *n - 1; j >= 1; --j) {
00260 /* Computing MIN */
00261                 i__1 = *kl, i__2 = *n - j;
00262                 lm = min(i__1,i__2);
00263                 zlacgv_(nrhs, &b[j + b_dim1], ldb);
00264                 z__1.r = -1., z__1.i = -0.;
00265                 zgemv_("Conjugate transpose", &lm, nrhs, &z__1, &b[j + 1 + 
00266                         b_dim1], ldb, &ab[kd + 1 + j * ab_dim1], &c__1, &c_b1, 
00267                          &b[j + b_dim1], ldb);
00268                 zlacgv_(nrhs, &b[j + b_dim1], ldb);
00269                 l = ipiv[j];
00270                 if (l != j) {
00271                     zswap_(nrhs, &b[l + b_dim1], ldb, &b[j + b_dim1], ldb);
00272                 }
00273 /* L60: */
00274             }
00275         }
00276     }
00277     return 0;
00278 
00279 /*     End of ZGBTRS */
00280 
00281 } /* zgbtrs_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:30