swiftnav: zlagge.c Source File

Go to the documentation of this file.
00001 /* zlagge.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 = {0.,0.};
00019 static doublecomplex c_b2 = {1.,0.};
00020 static integer c__3 = 3;
00021 static integer c__1 = 1;
00022 
00023 /* Subroutine */ int zlagge_(integer *m, integer *n, integer *kl, integer *ku, 
00024          doublereal *d__, doublecomplex *a, integer *lda, integer *iseed, 
00025         doublecomplex *work, integer *info)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, i__1, i__2, i__3;
00029     doublereal d__1;
00030     doublecomplex z__1;
00031 
00032     /* Builtin functions */
00033     double z_abs(doublecomplex *);
00034     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00035 
00036     /* Local variables */
00037     integer i__, j;
00038     doublecomplex wa, wb;
00039     doublereal wn;
00040     doublecomplex tau;
00041     extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 
00042             doublecomplex *, integer *, doublecomplex *, integer *, 
00043             doublecomplex *, integer *), zscal_(integer *, doublecomplex *, 
00044             doublecomplex *, integer *), zgemv_(char *, integer *, integer *, 
00045             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00046             integer *, doublecomplex *, doublecomplex *, integer *);
00047     extern doublereal dznrm2_(integer *, doublecomplex *, integer *);
00048     extern /* Subroutine */ int xerbla_(char *, integer *), zlacgv_(
00049             integer *, doublecomplex *, integer *), zlarnv_(integer *, 
00050             integer *, integer *, doublecomplex *);
00051 
00052 
00053 /*  -- LAPACK auxiliary test routine (version 3.1) -- */
00054 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00055 /*     November 2006 */
00056 
00057 /*     .. Scalar Arguments .. */
00058 /*     .. */
00059 /*     .. Array Arguments .. */
00060 /*     .. */
00061 
00062 /*  Purpose */
00063 /*  ======= */
00064 
00065 /*  ZLAGGE generates a complex general m by n matrix A, by pre- and post- */
00066 /*  multiplying a real diagonal matrix D with random unitary matrices: */
00067 /*  A = U*D*V. The lower and upper bandwidths may then be reduced to */
00068 /*  kl and ku by additional unitary transformations. */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*  M       (input) INTEGER */
00074 /*          The number of rows of the matrix A.  M >= 0. */
00075 
00076 /*  N       (input) INTEGER */
00077 /*          The number of columns of the matrix A.  N >= 0. */
00078 
00079 /*  KL      (input) INTEGER */
00080 /*          The number of nonzero subdiagonals within the band of A. */
00081 /*          0 <= KL <= M-1. */
00082 
00083 /*  KU      (input) INTEGER */
00084 /*          The number of nonzero superdiagonals within the band of A. */
00085 /*          0 <= KU <= N-1. */
00086 
00087 /*  D       (input) DOUBLE PRECISION array, dimension (min(M,N)) */
00088 /*          The diagonal elements of the diagonal matrix D. */
00089 
00090 /*  A       (output) COMPLEX*16 array, dimension (LDA,N) */
00091 /*          The generated m by n matrix A. */
00092 
00093 /*  LDA     (input) INTEGER */
00094 /*          The leading dimension of the array A.  LDA >= M. */
00095 
00096 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00097 /*          On entry, the seed of the random number generator; the array */
00098 /*          elements must be between 0 and 4095, and ISEED(4) must be */
00099 /*          odd. */
00100 /*          On exit, the seed is updated. */
00101 
00102 /*  WORK    (workspace) COMPLEX*16 array, dimension (M+N) */
00103 
00104 /*  INFO    (output) INTEGER */
00105 /*          = 0: successful exit */
00106 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00107 
00108 /*  ===================================================================== */
00109 
00110 /*     .. Parameters .. */
00111 /*     .. */
00112 /*     .. Local Scalars .. */
00113 /*     .. */
00114 /*     .. External Subroutines .. */
00115 /*     .. */
00116 /*     .. Intrinsic Functions .. */
00117 /*     .. */
00118 /*     .. External Functions .. */
00119 /*     .. */
00120 /*     .. Executable Statements .. */
00121 
00122 /*     Test the input arguments */
00123 
00124     /* Parameter adjustments */
00125     --d__;
00126     a_dim1 = *lda;
00127     a_offset = 1 + a_dim1;
00128     a -= a_offset;
00129     --iseed;
00130     --work;
00131 
00132     /* Function Body */
00133     *info = 0;
00134     if (*m < 0) {
00135         *info = -1;
00136     } else if (*n < 0) {
00137         *info = -2;
00138     } else if (*kl < 0 || *kl > *m - 1) {
00139         *info = -3;
00140     } else if (*ku < 0 || *ku > *n - 1) {
00141         *info = -4;
00142     } else if (*lda < max(1,*m)) {
00143         *info = -7;
00144     }
00145     if (*info < 0) {
00146         i__1 = -(*info);
00147         xerbla_("ZLAGGE", &i__1);
00148         return 0;
00149     }
00150 
00151 /*     initialize A to diagonal matrix */
00152 
00153     i__1 = *n;
00154     for (j = 1; j <= i__1; ++j) {
00155         i__2 = *m;
00156         for (i__ = 1; i__ <= i__2; ++i__) {
00157             i__3 = i__ + j * a_dim1;
00158             a[i__3].r = 0., a[i__3].i = 0.;
00159 /* L10: */
00160         }
00161 /* L20: */
00162     }
00163     i__1 = min(*m,*n);
00164     for (i__ = 1; i__ <= i__1; ++i__) {
00165         i__2 = i__ + i__ * a_dim1;
00166         i__3 = i__;
00167         a[i__2].r = d__[i__3], a[i__2].i = 0.;
00168 /* L30: */
00169     }
00170 
00171 /*     pre- and post-multiply A by random unitary matrices */
00172 
00173     for (i__ = min(*m,*n); i__ >= 1; --i__) {
00174         if (i__ < *m) {
00175 
00176 /*           generate random reflection */
00177 
00178             i__1 = *m - i__ + 1;
00179             zlarnv_(&c__3, &iseed[1], &i__1, &work[1]);
00180             i__1 = *m - i__ + 1;
00181             wn = dznrm2_(&i__1, &work[1], &c__1);
00182             d__1 = wn / z_abs(&work[1]);
00183             z__1.r = d__1 * work[1].r, z__1.i = d__1 * work[1].i;
00184             wa.r = z__1.r, wa.i = z__1.i;
00185             if (wn == 0.) {
00186                 tau.r = 0., tau.i = 0.;
00187             } else {
00188                 z__1.r = work[1].r + wa.r, z__1.i = work[1].i + wa.i;
00189                 wb.r = z__1.r, wb.i = z__1.i;
00190                 i__1 = *m - i__;
00191                 z_div(&z__1, &c_b2, &wb);
00192                 zscal_(&i__1, &z__1, &work[2], &c__1);
00193                 work[1].r = 1., work[1].i = 0.;
00194                 z_div(&z__1, &wb, &wa);
00195                 d__1 = z__1.r;
00196                 tau.r = d__1, tau.i = 0.;
00197             }
00198 
00199 /*           multiply A(i:m,i:n) by random reflection from the left */
00200 
00201             i__1 = *m - i__ + 1;
00202             i__2 = *n - i__ + 1;
00203             zgemv_("Conjugate transpose", &i__1, &i__2, &c_b2, &a[i__ + i__ * 
00204                     a_dim1], lda, &work[1], &c__1, &c_b1, &work[*m + 1], &
00205                     c__1);
00206             i__1 = *m - i__ + 1;
00207             i__2 = *n - i__ + 1;
00208             z__1.r = -tau.r, z__1.i = -tau.i;
00209             zgerc_(&i__1, &i__2, &z__1, &work[1], &c__1, &work[*m + 1], &c__1, 
00210                      &a[i__ + i__ * a_dim1], lda);
00211         }
00212         if (i__ < *n) {
00213 
00214 /*           generate random reflection */
00215 
00216             i__1 = *n - i__ + 1;
00217             zlarnv_(&c__3, &iseed[1], &i__1, &work[1]);
00218             i__1 = *n - i__ + 1;
00219             wn = dznrm2_(&i__1, &work[1], &c__1);
00220             d__1 = wn / z_abs(&work[1]);
00221             z__1.r = d__1 * work[1].r, z__1.i = d__1 * work[1].i;
00222             wa.r = z__1.r, wa.i = z__1.i;
00223             if (wn == 0.) {
00224                 tau.r = 0., tau.i = 0.;
00225             } else {
00226                 z__1.r = work[1].r + wa.r, z__1.i = work[1].i + wa.i;
00227                 wb.r = z__1.r, wb.i = z__1.i;
00228                 i__1 = *n - i__;
00229                 z_div(&z__1, &c_b2, &wb);
00230                 zscal_(&i__1, &z__1, &work[2], &c__1);
00231                 work[1].r = 1., work[1].i = 0.;
00232                 z_div(&z__1, &wb, &wa);
00233                 d__1 = z__1.r;
00234                 tau.r = d__1, tau.i = 0.;
00235             }
00236 
00237 /*           multiply A(i:m,i:n) by random reflection from the right */
00238 
00239             i__1 = *m - i__ + 1;
00240             i__2 = *n - i__ + 1;
00241             zgemv_("No transpose", &i__1, &i__2, &c_b2, &a[i__ + i__ * a_dim1]
00242 , lda, &work[1], &c__1, &c_b1, &work[*n + 1], &c__1);
00243             i__1 = *m - i__ + 1;
00244             i__2 = *n - i__ + 1;
00245             z__1.r = -tau.r, z__1.i = -tau.i;
00246             zgerc_(&i__1, &i__2, &z__1, &work[*n + 1], &c__1, &work[1], &c__1, 
00247                      &a[i__ + i__ * a_dim1], lda);
00248         }
00249 /* L40: */
00250     }
00251 
00252 /*     Reduce number of subdiagonals to KL and number of superdiagonals */
00253 /*     to KU */
00254 
00255 /* Computing MAX */
00256     i__2 = *m - 1 - *kl, i__3 = *n - 1 - *ku;
00257     i__1 = max(i__2,i__3);
00258     for (i__ = 1; i__ <= i__1; ++i__) {
00259         if (*kl <= *ku) {
00260 
00261 /*           annihilate subdiagonal elements first (necessary if KL = 0) */
00262 
00263 /* Computing MIN */
00264             i__2 = *m - 1 - *kl;
00265             if (i__ <= min(i__2,*n)) {
00266 
00267 /*              generate reflection to annihilate A(kl+i+1:m,i) */
00268 
00269                 i__2 = *m - *kl - i__ + 1;
00270                 wn = dznrm2_(&i__2, &a[*kl + i__ + i__ * a_dim1], &c__1);
00271                 d__1 = wn / z_abs(&a[*kl + i__ + i__ * a_dim1]);
00272                 i__2 = *kl + i__ + i__ * a_dim1;
00273                 z__1.r = d__1 * a[i__2].r, z__1.i = d__1 * a[i__2].i;
00274                 wa.r = z__1.r, wa.i = z__1.i;
00275                 if (wn == 0.) {
00276                     tau.r = 0., tau.i = 0.;
00277                 } else {
00278                     i__2 = *kl + i__ + i__ * a_dim1;
00279                     z__1.r = a[i__2].r + wa.r, z__1.i = a[i__2].i + wa.i;
00280                     wb.r = z__1.r, wb.i = z__1.i;
00281                     i__2 = *m - *kl - i__;
00282                     z_div(&z__1, &c_b2, &wb);
00283                     zscal_(&i__2, &z__1, &a[*kl + i__ + 1 + i__ * a_dim1], &
00284                             c__1);
00285                     i__2 = *kl + i__ + i__ * a_dim1;
00286                     a[i__2].r = 1., a[i__2].i = 0.;
00287                     z_div(&z__1, &wb, &wa);
00288                     d__1 = z__1.r;
00289                     tau.r = d__1, tau.i = 0.;
00290                 }
00291 
00292 /*              apply reflection to A(kl+i:m,i+1:n) from the left */
00293 
00294                 i__2 = *m - *kl - i__ + 1;
00295                 i__3 = *n - i__;
00296                 zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*kl + 
00297                         i__ + (i__ + 1) * a_dim1], lda, &a[*kl + i__ + i__ * 
00298                         a_dim1], &c__1, &c_b1, &work[1], &c__1);
00299                 i__2 = *m - *kl - i__ + 1;
00300                 i__3 = *n - i__;
00301                 z__1.r = -tau.r, z__1.i = -tau.i;
00302                 zgerc_(&i__2, &i__3, &z__1, &a[*kl + i__ + i__ * a_dim1], &
00303                         c__1, &work[1], &c__1, &a[*kl + i__ + (i__ + 1) * 
00304                         a_dim1], lda);
00305                 i__2 = *kl + i__ + i__ * a_dim1;
00306                 z__1.r = -wa.r, z__1.i = -wa.i;
00307                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00308             }
00309 
00310 /* Computing MIN */
00311             i__2 = *n - 1 - *ku;
00312             if (i__ <= min(i__2,*m)) {
00313 
00314 /*              generate reflection to annihilate A(i,ku+i+1:n) */
00315 
00316                 i__2 = *n - *ku - i__ + 1;
00317                 wn = dznrm2_(&i__2, &a[i__ + (*ku + i__) * a_dim1], lda);
00318                 d__1 = wn / z_abs(&a[i__ + (*ku + i__) * a_dim1]);
00319                 i__2 = i__ + (*ku + i__) * a_dim1;
00320                 z__1.r = d__1 * a[i__2].r, z__1.i = d__1 * a[i__2].i;
00321                 wa.r = z__1.r, wa.i = z__1.i;
00322                 if (wn == 0.) {
00323                     tau.r = 0., tau.i = 0.;
00324                 } else {
00325                     i__2 = i__ + (*ku + i__) * a_dim1;
00326                     z__1.r = a[i__2].r + wa.r, z__1.i = a[i__2].i + wa.i;
00327                     wb.r = z__1.r, wb.i = z__1.i;
00328                     i__2 = *n - *ku - i__;
00329                     z_div(&z__1, &c_b2, &wb);
00330                     zscal_(&i__2, &z__1, &a[i__ + (*ku + i__ + 1) * a_dim1], 
00331                             lda);
00332                     i__2 = i__ + (*ku + i__) * a_dim1;
00333                     a[i__2].r = 1., a[i__2].i = 0.;
00334                     z_div(&z__1, &wb, &wa);
00335                     d__1 = z__1.r;
00336                     tau.r = d__1, tau.i = 0.;
00337                 }
00338 
00339 /*              apply reflection to A(i+1:m,ku+i:n) from the right */
00340 
00341                 i__2 = *n - *ku - i__ + 1;
00342                 zlacgv_(&i__2, &a[i__ + (*ku + i__) * a_dim1], lda);
00343                 i__2 = *m - i__;
00344                 i__3 = *n - *ku - i__ + 1;
00345                 zgemv_("No transpose", &i__2, &i__3, &c_b2, &a[i__ + 1 + (*ku 
00346                         + i__) * a_dim1], lda, &a[i__ + (*ku + i__) * a_dim1], 
00347                          lda, &c_b1, &work[1], &c__1);
00348                 i__2 = *m - i__;
00349                 i__3 = *n - *ku - i__ + 1;
00350                 z__1.r = -tau.r, z__1.i = -tau.i;
00351                 zgerc_(&i__2, &i__3, &z__1, &work[1], &c__1, &a[i__ + (*ku + 
00352                         i__) * a_dim1], lda, &a[i__ + 1 + (*ku + i__) * 
00353                         a_dim1], lda);
00354                 i__2 = i__ + (*ku + i__) * a_dim1;
00355                 z__1.r = -wa.r, z__1.i = -wa.i;
00356                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00357             }
00358         } else {
00359 
00360 /*           annihilate superdiagonal elements first (necessary if */
00361 /*           KU = 0) */
00362 
00363 /* Computing MIN */
00364             i__2 = *n - 1 - *ku;
00365             if (i__ <= min(i__2,*m)) {
00366 
00367 /*              generate reflection to annihilate A(i,ku+i+1:n) */
00368 
00369                 i__2 = *n - *ku - i__ + 1;
00370                 wn = dznrm2_(&i__2, &a[i__ + (*ku + i__) * a_dim1], lda);
00371                 d__1 = wn / z_abs(&a[i__ + (*ku + i__) * a_dim1]);
00372                 i__2 = i__ + (*ku + i__) * a_dim1;
00373                 z__1.r = d__1 * a[i__2].r, z__1.i = d__1 * a[i__2].i;
00374                 wa.r = z__1.r, wa.i = z__1.i;
00375                 if (wn == 0.) {
00376                     tau.r = 0., tau.i = 0.;
00377                 } else {
00378                     i__2 = i__ + (*ku + i__) * a_dim1;
00379                     z__1.r = a[i__2].r + wa.r, z__1.i = a[i__2].i + wa.i;
00380                     wb.r = z__1.r, wb.i = z__1.i;
00381                     i__2 = *n - *ku - i__;
00382                     z_div(&z__1, &c_b2, &wb);
00383                     zscal_(&i__2, &z__1, &a[i__ + (*ku + i__ + 1) * a_dim1], 
00384                             lda);
00385                     i__2 = i__ + (*ku + i__) * a_dim1;
00386                     a[i__2].r = 1., a[i__2].i = 0.;
00387                     z_div(&z__1, &wb, &wa);
00388                     d__1 = z__1.r;
00389                     tau.r = d__1, tau.i = 0.;
00390                 }
00391 
00392 /*              apply reflection to A(i+1:m,ku+i:n) from the right */
00393 
00394                 i__2 = *n - *ku - i__ + 1;
00395                 zlacgv_(&i__2, &a[i__ + (*ku + i__) * a_dim1], lda);
00396                 i__2 = *m - i__;
00397                 i__3 = *n - *ku - i__ + 1;
00398                 zgemv_("No transpose", &i__2, &i__3, &c_b2, &a[i__ + 1 + (*ku 
00399                         + i__) * a_dim1], lda, &a[i__ + (*ku + i__) * a_dim1], 
00400                          lda, &c_b1, &work[1], &c__1);
00401                 i__2 = *m - i__;
00402                 i__3 = *n - *ku - i__ + 1;
00403                 z__1.r = -tau.r, z__1.i = -tau.i;
00404                 zgerc_(&i__2, &i__3, &z__1, &work[1], &c__1, &a[i__ + (*ku + 
00405                         i__) * a_dim1], lda, &a[i__ + 1 + (*ku + i__) * 
00406                         a_dim1], lda);
00407                 i__2 = i__ + (*ku + i__) * a_dim1;
00408                 z__1.r = -wa.r, z__1.i = -wa.i;
00409                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00410             }
00411 
00412 /* Computing MIN */
00413             i__2 = *m - 1 - *kl;
00414             if (i__ <= min(i__2,*n)) {
00415 
00416 /*              generate reflection to annihilate A(kl+i+1:m,i) */
00417 
00418                 i__2 = *m - *kl - i__ + 1;
00419                 wn = dznrm2_(&i__2, &a[*kl + i__ + i__ * a_dim1], &c__1);
00420                 d__1 = wn / z_abs(&a[*kl + i__ + i__ * a_dim1]);
00421                 i__2 = *kl + i__ + i__ * a_dim1;
00422                 z__1.r = d__1 * a[i__2].r, z__1.i = d__1 * a[i__2].i;
00423                 wa.r = z__1.r, wa.i = z__1.i;
00424                 if (wn == 0.) {
00425                     tau.r = 0., tau.i = 0.;
00426                 } else {
00427                     i__2 = *kl + i__ + i__ * a_dim1;
00428                     z__1.r = a[i__2].r + wa.r, z__1.i = a[i__2].i + wa.i;
00429                     wb.r = z__1.r, wb.i = z__1.i;
00430                     i__2 = *m - *kl - i__;
00431                     z_div(&z__1, &c_b2, &wb);
00432                     zscal_(&i__2, &z__1, &a[*kl + i__ + 1 + i__ * a_dim1], &
00433                             c__1);
00434                     i__2 = *kl + i__ + i__ * a_dim1;
00435                     a[i__2].r = 1., a[i__2].i = 0.;
00436                     z_div(&z__1, &wb, &wa);
00437                     d__1 = z__1.r;
00438                     tau.r = d__1, tau.i = 0.;
00439                 }
00440 
00441 /*              apply reflection to A(kl+i:m,i+1:n) from the left */
00442 
00443                 i__2 = *m - *kl - i__ + 1;
00444                 i__3 = *n - i__;
00445                 zgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*kl + 
00446                         i__ + (i__ + 1) * a_dim1], lda, &a[*kl + i__ + i__ * 
00447                         a_dim1], &c__1, &c_b1, &work[1], &c__1);
00448                 i__2 = *m - *kl - i__ + 1;
00449                 i__3 = *n - i__;
00450                 z__1.r = -tau.r, z__1.i = -tau.i;
00451                 zgerc_(&i__2, &i__3, &z__1, &a[*kl + i__ + i__ * a_dim1], &
00452                         c__1, &work[1], &c__1, &a[*kl + i__ + (i__ + 1) * 
00453                         a_dim1], lda);
00454                 i__2 = *kl + i__ + i__ * a_dim1;
00455                 z__1.r = -wa.r, z__1.i = -wa.i;
00456                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00457             }
00458         }
00459 
00460         i__2 = *m;
00461         for (j = *kl + i__ + 1; j <= i__2; ++j) {
00462             i__3 = j + i__ * a_dim1;
00463             a[i__3].r = 0., a[i__3].i = 0.;
00464 /* L50: */
00465         }
00466 
00467         i__2 = *n;
00468         for (j = *ku + i__ + 1; j <= i__2; ++j) {
00469             i__3 = i__ + j * a_dim1;
00470             a[i__3].r = 0., a[i__3].i = 0.;
00471 /* L60: */
00472         }
00473 /* L70: */
00474     }
00475     return 0;
00476 
00477 /*     End of ZLAGGE */
00478 
00479 } /* zlagge_ */