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 integer c__3 = 3;
00019 static integer c__1 = 1;
00020 static real c_b12 = 0.f;
00021 static real c_b19 = -1.f;
00022 static real c_b26 = 1.f;
00023
00024 int slagsy_(integer *n, integer *k, real *d__, real *a,
00025 integer *lda, integer *iseed, real *work, integer *info)
00026 {
00027
00028 integer a_dim1, a_offset, i__1, i__2, i__3;
00029 real r__1;
00030
00031
00032 double r_sign(real *, real *);
00033
00034
00035 integer i__, j;
00036 real wa, wb, wn, tau;
00037 extern int sger_(integer *, integer *, real *, real *,
00038 integer *, real *, integer *, real *, integer *);
00039 extern doublereal sdot_(integer *, real *, integer *, real *, integer *),
00040 snrm2_(integer *, real *, integer *);
00041 extern int ssyr2_(char *, integer *, real *, real *,
00042 integer *, real *, integer *, real *, integer *);
00043 real alpha;
00044 extern int sscal_(integer *, real *, real *, integer *),
00045 sgemv_(char *, integer *, integer *, real *, real *, integer *,
00046 real *, integer *, real *, real *, integer *), saxpy_(
00047 integer *, real *, real *, integer *, real *, integer *), ssymv_(
00048 char *, integer *, real *, real *, integer *, real *, integer *,
00049 real *, real *, integer *), xerbla_(char *, integer *), slarnv_(integer *, integer *, integer *, real *);
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 --d__;
00119 a_dim1 = *lda;
00120 a_offset = 1 + a_dim1;
00121 a -= a_offset;
00122 --iseed;
00123 --work;
00124
00125
00126 *info = 0;
00127 if (*n < 0) {
00128 *info = -1;
00129 } else if (*k < 0 || *k > *n - 1) {
00130 *info = -2;
00131 } else if (*lda < max(1,*n)) {
00132 *info = -5;
00133 }
00134 if (*info < 0) {
00135 i__1 = -(*info);
00136 xerbla_("SLAGSY", &i__1);
00137 return 0;
00138 }
00139
00140
00141
00142 i__1 = *n;
00143 for (j = 1; j <= i__1; ++j) {
00144 i__2 = *n;
00145 for (i__ = j + 1; i__ <= i__2; ++i__) {
00146 a[i__ + j * a_dim1] = 0.f;
00147
00148 }
00149
00150 }
00151 i__1 = *n;
00152 for (i__ = 1; i__ <= i__1; ++i__) {
00153 a[i__ + i__ * a_dim1] = d__[i__];
00154
00155 }
00156
00157
00158
00159 for (i__ = *n - 1; i__ >= 1; --i__) {
00160
00161
00162
00163 i__1 = *n - i__ + 1;
00164 slarnv_(&c__3, &iseed[1], &i__1, &work[1]);
00165 i__1 = *n - i__ + 1;
00166 wn = snrm2_(&i__1, &work[1], &c__1);
00167 wa = r_sign(&wn, &work[1]);
00168 if (wn == 0.f) {
00169 tau = 0.f;
00170 } else {
00171 wb = work[1] + wa;
00172 i__1 = *n - i__;
00173 r__1 = 1.f / wb;
00174 sscal_(&i__1, &r__1, &work[2], &c__1);
00175 work[1] = 1.f;
00176 tau = wb / wa;
00177 }
00178
00179
00180
00181
00182
00183
00184 i__1 = *n - i__ + 1;
00185 ssymv_("Lower", &i__1, &tau, &a[i__ + i__ * a_dim1], lda, &work[1], &
00186 c__1, &c_b12, &work[*n + 1], &c__1);
00187
00188
00189
00190 i__1 = *n - i__ + 1;
00191 alpha = tau * -.5f * sdot_(&i__1, &work[*n + 1], &c__1, &work[1], &
00192 c__1);
00193 i__1 = *n - i__ + 1;
00194 saxpy_(&i__1, &alpha, &work[1], &c__1, &work[*n + 1], &c__1);
00195
00196
00197
00198 i__1 = *n - i__ + 1;
00199 ssyr2_("Lower", &i__1, &c_b19, &work[1], &c__1, &work[*n + 1], &c__1,
00200 &a[i__ + i__ * a_dim1], lda);
00201
00202 }
00203
00204
00205
00206 i__1 = *n - 1 - *k;
00207 for (i__ = 1; i__ <= i__1; ++i__) {
00208
00209
00210
00211 i__2 = *n - *k - i__ + 1;
00212 wn = snrm2_(&i__2, &a[*k + i__ + i__ * a_dim1], &c__1);
00213 wa = r_sign(&wn, &a[*k + i__ + i__ * a_dim1]);
00214 if (wn == 0.f) {
00215 tau = 0.f;
00216 } else {
00217 wb = a[*k + i__ + i__ * a_dim1] + wa;
00218 i__2 = *n - *k - i__;
00219 r__1 = 1.f / wb;
00220 sscal_(&i__2, &r__1, &a[*k + i__ + 1 + i__ * a_dim1], &c__1);
00221 a[*k + i__ + i__ * a_dim1] = 1.f;
00222 tau = wb / wa;
00223 }
00224
00225
00226
00227 i__2 = *n - *k - i__ + 1;
00228 i__3 = *k - 1;
00229 sgemv_("Transpose", &i__2, &i__3, &c_b26, &a[*k + i__ + (i__ + 1) *
00230 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b12, &
00231 work[1], &c__1);
00232 i__2 = *n - *k - i__ + 1;
00233 i__3 = *k - 1;
00234 r__1 = -tau;
00235 sger_(&i__2, &i__3, &r__1, &a[*k + i__ + i__ * a_dim1], &c__1, &work[
00236 1], &c__1, &a[*k + i__ + (i__ + 1) * a_dim1], lda);
00237
00238
00239
00240
00241
00242 i__2 = *n - *k - i__ + 1;
00243 ssymv_("Lower", &i__2, &tau, &a[*k + i__ + (*k + i__) * a_dim1], lda,
00244 &a[*k + i__ + i__ * a_dim1], &c__1, &c_b12, &work[1], &c__1);
00245
00246
00247
00248 i__2 = *n - *k - i__ + 1;
00249 alpha = tau * -.5f * sdot_(&i__2, &work[1], &c__1, &a[*k + i__ + i__ *
00250 a_dim1], &c__1);
00251 i__2 = *n - *k - i__ + 1;
00252 saxpy_(&i__2, &alpha, &a[*k + i__ + i__ * a_dim1], &c__1, &work[1], &
00253 c__1);
00254
00255
00256
00257 i__2 = *n - *k - i__ + 1;
00258 ssyr2_("Lower", &i__2, &c_b19, &a[*k + i__ + i__ * a_dim1], &c__1, &
00259 work[1], &c__1, &a[*k + i__ + (*k + i__) * a_dim1], lda);
00260
00261 a[*k + i__ + i__ * a_dim1] = -wa;
00262 i__2 = *n;
00263 for (j = *k + i__ + 1; j <= i__2; ++j) {
00264 a[j + i__ * a_dim1] = 0.f;
00265
00266 }
00267
00268 }
00269
00270
00271
00272 i__1 = *n;
00273 for (j = 1; j <= i__1; ++j) {
00274 i__2 = *n;
00275 for (i__ = j + 1; i__ <= i__2; ++i__) {
00276 a[j + i__ * a_dim1] = a[i__ + j * a_dim1];
00277
00278 }
00279
00280 }
00281 return 0;
00282
00283
00284
00285 }