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 complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 static integer c__1 = 1;
00021
00022 int clahr2_(integer *n, integer *k, integer *nb, complex *a,
00023 integer *lda, complex *tau, complex *t, integer *ldt, complex *y,
00024 integer *ldy)
00025 {
00026
00027 integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2,
00028 i__3;
00029 complex q__1;
00030
00031
00032 integer i__;
00033 complex ei;
00034 extern int cscal_(integer *, complex *, complex *,
00035 integer *), cgemm_(char *, char *, integer *, integer *, integer *
00036 , complex *, complex *, integer *, complex *, integer *, complex *
00037 , complex *, integer *), cgemv_(char *, integer *,
00038 integer *, complex *, complex *, integer *, complex *, integer *,
00039 complex *, complex *, integer *), ccopy_(integer *,
00040 complex *, integer *, complex *, integer *), ctrmm_(char *, char *
00041 , char *, char *, integer *, integer *, complex *, complex *,
00042 integer *, complex *, integer *),
00043 caxpy_(integer *, complex *, complex *, integer *, complex *,
00044 integer *), ctrmv_(char *, char *, char *, integer *, complex *,
00045 integer *, complex *, integer *), clarfg_(
00046 integer *, complex *, complex *, integer *, complex *), clacgv_(
00047 integer *, complex *, integer *), clacpy_(char *, integer *,
00048 integer *, complex *, integer *, complex *, integer *);
00049
00050
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
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 --tau;
00170 a_dim1 = *lda;
00171 a_offset = 1 + a_dim1;
00172 a -= a_offset;
00173 t_dim1 = *ldt;
00174 t_offset = 1 + t_dim1;
00175 t -= t_offset;
00176 y_dim1 = *ldy;
00177 y_offset = 1 + y_dim1;
00178 y -= y_offset;
00179
00180
00181 if (*n <= 1) {
00182 return 0;
00183 }
00184
00185 i__1 = *nb;
00186 for (i__ = 1; i__ <= i__1; ++i__) {
00187 if (i__ > 1) {
00188
00189
00190
00191
00192
00193 i__2 = i__ - 1;
00194 clacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
00195 i__2 = *n - *k;
00196 i__3 = i__ - 1;
00197 q__1.r = -1.f, q__1.i = -0.f;
00198 cgemv_("NO TRANSPOSE", &i__2, &i__3, &q__1, &y[*k + 1 + y_dim1],
00199 ldy, &a[*k + i__ - 1 + a_dim1], lda, &c_b2, &a[*k + 1 +
00200 i__ * a_dim1], &c__1);
00201 i__2 = i__ - 1;
00202 clacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 i__2 = i__ - 1;
00215 ccopy_(&i__2, &a[*k + 1 + i__ * a_dim1], &c__1, &t[*nb * t_dim1 +
00216 1], &c__1);
00217 i__2 = i__ - 1;
00218 ctrmv_("Lower", "Conjugate transpose", "UNIT", &i__2, &a[*k + 1 +
00219 a_dim1], lda, &t[*nb * t_dim1 + 1], &c__1);
00220
00221
00222
00223 i__2 = *n - *k - i__ + 1;
00224 i__3 = i__ - 1;
00225 cgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ +
00226 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b2, &
00227 t[*nb * t_dim1 + 1], &c__1);
00228
00229
00230
00231 i__2 = i__ - 1;
00232 ctrmv_("Upper", "Conjugate transpose", "NON-UNIT", &i__2, &t[
00233 t_offset], ldt, &t[*nb * t_dim1 + 1], &c__1);
00234
00235
00236
00237 i__2 = *n - *k - i__ + 1;
00238 i__3 = i__ - 1;
00239 q__1.r = -1.f, q__1.i = -0.f;
00240 cgemv_("NO TRANSPOSE", &i__2, &i__3, &q__1, &a[*k + i__ + a_dim1],
00241 lda, &t[*nb * t_dim1 + 1], &c__1, &c_b2, &a[*k + i__ +
00242 i__ * a_dim1], &c__1);
00243
00244
00245
00246 i__2 = i__ - 1;
00247 ctrmv_("Lower", "NO TRANSPOSE", "UNIT", &i__2, &a[*k + 1 + a_dim1]
00248 , lda, &t[*nb * t_dim1 + 1], &c__1);
00249 i__2 = i__ - 1;
00250 q__1.r = -1.f, q__1.i = -0.f;
00251 caxpy_(&i__2, &q__1, &t[*nb * t_dim1 + 1], &c__1, &a[*k + 1 + i__
00252 * a_dim1], &c__1);
00253
00254 i__2 = *k + i__ - 1 + (i__ - 1) * a_dim1;
00255 a[i__2].r = ei.r, a[i__2].i = ei.i;
00256 }
00257
00258
00259
00260
00261 i__2 = *n - *k - i__ + 1;
00262
00263 i__3 = *k + i__ + 1;
00264 clarfg_(&i__2, &a[*k + i__ + i__ * a_dim1], &a[min(i__3, *n)+ i__ *
00265 a_dim1], &c__1, &tau[i__]);
00266 i__2 = *k + i__ + i__ * a_dim1;
00267 ei.r = a[i__2].r, ei.i = a[i__2].i;
00268 i__2 = *k + i__ + i__ * a_dim1;
00269 a[i__2].r = 1.f, a[i__2].i = 0.f;
00270
00271
00272
00273 i__2 = *n - *k;
00274 i__3 = *n - *k - i__ + 1;
00275 cgemv_("NO TRANSPOSE", &i__2, &i__3, &c_b2, &a[*k + 1 + (i__ + 1) *
00276 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &y[*
00277 k + 1 + i__ * y_dim1], &c__1);
00278 i__2 = *n - *k - i__ + 1;
00279 i__3 = i__ - 1;
00280 cgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ +
00281 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &t[
00282 i__ * t_dim1 + 1], &c__1);
00283 i__2 = *n - *k;
00284 i__3 = i__ - 1;
00285 q__1.r = -1.f, q__1.i = -0.f;
00286 cgemv_("NO TRANSPOSE", &i__2, &i__3, &q__1, &y[*k + 1 + y_dim1], ldy,
00287 &t[i__ * t_dim1 + 1], &c__1, &c_b2, &y[*k + 1 + i__ * y_dim1],
00288 &c__1);
00289 i__2 = *n - *k;
00290 cscal_(&i__2, &tau[i__], &y[*k + 1 + i__ * y_dim1], &c__1);
00291
00292
00293
00294 i__2 = i__ - 1;
00295 i__3 = i__;
00296 q__1.r = -tau[i__3].r, q__1.i = -tau[i__3].i;
00297 cscal_(&i__2, &q__1, &t[i__ * t_dim1 + 1], &c__1);
00298 i__2 = i__ - 1;
00299 ctrmv_("Upper", "No Transpose", "NON-UNIT", &i__2, &t[t_offset], ldt,
00300 &t[i__ * t_dim1 + 1], &c__1)
00301 ;
00302 i__2 = i__ + i__ * t_dim1;
00303 i__3 = i__;
00304 t[i__2].r = tau[i__3].r, t[i__2].i = tau[i__3].i;
00305
00306
00307 }
00308 i__1 = *k + *nb + *nb * a_dim1;
00309 a[i__1].r = ei.r, a[i__1].i = ei.i;
00310
00311
00312
00313 clacpy_("ALL", k, nb, &a[(a_dim1 << 1) + 1], lda, &y[y_offset], ldy);
00314 ctrmm_("RIGHT", "Lower", "NO TRANSPOSE", "UNIT", k, nb, &c_b2, &a[*k + 1
00315 + a_dim1], lda, &y[y_offset], ldy);
00316 if (*n > *k + *nb) {
00317 i__1 = *n - *k - *nb;
00318 cgemm_("NO TRANSPOSE", "NO TRANSPOSE", k, nb, &i__1, &c_b2, &a[(*nb +
00319 2) * a_dim1 + 1], lda, &a[*k + 1 + *nb + a_dim1], lda, &c_b2,
00320 &y[y_offset], ldy);
00321 }
00322 ctrmm_("RIGHT", "Upper", "NO TRANSPOSE", "NON-UNIT", k, nb, &c_b2, &t[
00323 t_offset], ldt, &y[y_offset], ldy);
00324
00325 return 0;
00326
00327
00328
00329 }