00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
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 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
00028 integer a_dim1, a_offset, i__1, i__2, i__3;
00029 doublereal d__1;
00030 doublecomplex z__1;
00031
00032
00033 double z_abs(doublecomplex *);
00034 void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00035
00036
00037 integer i__, j;
00038 doublecomplex wa, wb;
00039 doublereal wn;
00040 doublecomplex tau;
00041 extern 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 int xerbla_(char *, integer *), zlacgv_(
00049 integer *, doublecomplex *, integer *), zlarnv_(integer *,
00050 integer *, integer *, doublecomplex *);
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 --d__;
00126 a_dim1 = *lda;
00127 a_offset = 1 + a_dim1;
00128 a -= a_offset;
00129 --iseed;
00130 --work;
00131
00132
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
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
00160 }
00161
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
00169 }
00170
00171
00172
00173 for (i__ = min(*m,*n); i__ >= 1; --i__) {
00174 if (i__ < *m) {
00175
00176
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
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
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
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
00250 }
00251
00252
00253
00254
00255
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
00262
00263
00264 i__2 = *m - 1 - *kl;
00265 if (i__ <= min(i__2,*n)) {
00266
00267
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
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
00311 i__2 = *n - 1 - *ku;
00312 if (i__ <= min(i__2,*m)) {
00313
00314
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
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
00361
00362
00363
00364 i__2 = *n - 1 - *ku;
00365 if (i__ <= min(i__2,*m)) {
00366
00367
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
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
00413 i__2 = *m - 1 - *kl;
00414 if (i__ <= min(i__2,*n)) {
00415
00416
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
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
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
00472 }
00473
00474 }
00475 return 0;
00476
00477
00478
00479 }