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 real c_b15 = 1.f;
00019 static integer c__1 = 1;
00020
00021 int slavsy_(char *uplo, char *trans, char *diag, integer *n,
00022 integer *nrhs, real *a, integer *lda, integer *ipiv, real *b, integer
00023 *ldb, integer *info)
00024 {
00025
00026 integer a_dim1, a_offset, b_dim1, b_offset, i__1;
00027
00028
00029 integer j, k;
00030 real t1, t2, d11, d12, d21, d22;
00031 integer kp;
00032 extern int sger_(integer *, integer *, real *, real *,
00033 integer *, real *, integer *, real *, integer *);
00034 extern logical lsame_(char *, char *);
00035 extern int sscal_(integer *, real *, real *, integer *),
00036 sgemv_(char *, integer *, integer *, real *, real *, integer *,
00037 real *, integer *, real *, real *, integer *), sswap_(
00038 integer *, real *, integer *, real *, integer *), xerbla_(char *,
00039 integer *);
00040 logical nounit;
00041
00042
00043
00044
00045
00046
00047
00048
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 a_dim1 = *lda;
00132 a_offset = 1 + a_dim1;
00133 a -= a_offset;
00134 --ipiv;
00135 b_dim1 = *ldb;
00136 b_offset = 1 + b_dim1;
00137 b -= b_offset;
00138
00139
00140 *info = 0;
00141 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00142 *info = -1;
00143 } else if (! lsame_(trans, "N") && ! lsame_(trans,
00144 "T") && ! lsame_(trans, "C")) {
00145 *info = -2;
00146 } else if (! lsame_(diag, "U") && ! lsame_(diag,
00147 "N")) {
00148 *info = -3;
00149 } else if (*n < 0) {
00150 *info = -4;
00151 } else if (*lda < max(1,*n)) {
00152 *info = -6;
00153 } else if (*ldb < max(1,*n)) {
00154 *info = -9;
00155 }
00156 if (*info != 0) {
00157 i__1 = -(*info);
00158 xerbla_("SLAVSY ", &i__1);
00159 return 0;
00160 }
00161
00162
00163
00164 if (*n == 0) {
00165 return 0;
00166 }
00167
00168 nounit = lsame_(diag, "N");
00169
00170
00171
00172
00173
00174 if (lsame_(trans, "N")) {
00175
00176
00177
00178
00179 if (lsame_(uplo, "U")) {
00180
00181
00182
00183 k = 1;
00184 L10:
00185 if (k > *n) {
00186 goto L30;
00187 }
00188 if (ipiv[k] > 0) {
00189
00190
00191
00192
00193
00194 if (nounit) {
00195 sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00196 }
00197
00198
00199
00200 if (k > 1) {
00201
00202
00203
00204 i__1 = k - 1;
00205 sger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
00206 k + b_dim1], ldb, &b[b_dim1 + 1], ldb);
00207
00208
00209
00210 kp = ipiv[k];
00211 if (kp != k) {
00212 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00213 ldb);
00214 }
00215 }
00216 ++k;
00217 } else {
00218
00219
00220
00221
00222
00223 if (nounit) {
00224 d11 = a[k + k * a_dim1];
00225 d22 = a[k + 1 + (k + 1) * a_dim1];
00226 d12 = a[k + (k + 1) * a_dim1];
00227 d21 = d12;
00228 i__1 = *nrhs;
00229 for (j = 1; j <= i__1; ++j) {
00230 t1 = b[k + j * b_dim1];
00231 t2 = b[k + 1 + j * b_dim1];
00232 b[k + j * b_dim1] = d11 * t1 + d12 * t2;
00233 b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
00234
00235 }
00236 }
00237
00238
00239
00240 if (k > 1) {
00241
00242
00243
00244 i__1 = k - 1;
00245 sger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
00246 k + b_dim1], ldb, &b[b_dim1 + 1], ldb);
00247 i__1 = k - 1;
00248 sger_(&i__1, nrhs, &c_b15, &a[(k + 1) * a_dim1 + 1], &
00249 c__1, &b[k + 1 + b_dim1], ldb, &b[b_dim1 + 1],
00250 ldb);
00251
00252
00253
00254 kp = (i__1 = ipiv[k], abs(i__1));
00255 if (kp != k) {
00256 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00257 ldb);
00258 }
00259 }
00260 k += 2;
00261 }
00262 goto L10;
00263 L30:
00264
00265
00266
00267
00268 ;
00269 } else {
00270
00271
00272
00273 k = *n;
00274 L40:
00275 if (k < 1) {
00276 goto L60;
00277 }
00278
00279
00280
00281
00282 if (ipiv[k] > 0) {
00283
00284
00285
00286
00287
00288 if (nounit) {
00289 sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00290 }
00291
00292
00293
00294 if (k != *n) {
00295 kp = ipiv[k];
00296
00297
00298
00299 i__1 = *n - k;
00300 sger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1,
00301 &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00302
00303
00304
00305
00306 if (kp != k) {
00307 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00308 ldb);
00309 }
00310 }
00311 --k;
00312
00313 } else {
00314
00315
00316
00317
00318
00319 if (nounit) {
00320 d11 = a[k - 1 + (k - 1) * a_dim1];
00321 d22 = a[k + k * a_dim1];
00322 d21 = a[k + (k - 1) * a_dim1];
00323 d12 = d21;
00324 i__1 = *nrhs;
00325 for (j = 1; j <= i__1; ++j) {
00326 t1 = b[k - 1 + j * b_dim1];
00327 t2 = b[k + j * b_dim1];
00328 b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
00329 b[k + j * b_dim1] = d21 * t1 + d22 * t2;
00330
00331 }
00332 }
00333
00334
00335
00336 if (k != *n) {
00337
00338
00339
00340 i__1 = *n - k;
00341 sger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1,
00342 &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00343 i__1 = *n - k;
00344 sger_(&i__1, nrhs, &c_b15, &a[k + 1 + (k - 1) * a_dim1], &
00345 c__1, &b[k - 1 + b_dim1], ldb, &b[k + 1 + b_dim1],
00346 ldb);
00347
00348
00349
00350
00351 kp = (i__1 = ipiv[k], abs(i__1));
00352 if (kp != k) {
00353 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00354 ldb);
00355 }
00356 }
00357 k += -2;
00358 }
00359 goto L40;
00360 L60:
00361 ;
00362 }
00363
00364
00365
00366
00367
00368 } else {
00369
00370
00371
00372
00373
00374 if (lsame_(uplo, "U")) {
00375
00376
00377
00378 k = *n;
00379 L70:
00380 if (k < 1) {
00381 goto L90;
00382 }
00383
00384
00385
00386 if (ipiv[k] > 0) {
00387 if (k > 1) {
00388
00389
00390
00391 kp = ipiv[k];
00392 if (kp != k) {
00393 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00394 ldb);
00395 }
00396
00397
00398
00399 i__1 = k - 1;
00400 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset],
00401 ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k +
00402 b_dim1], ldb);
00403 }
00404 if (nounit) {
00405 sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00406 }
00407 --k;
00408
00409
00410
00411 } else {
00412 if (k > 2) {
00413
00414
00415
00416 kp = (i__1 = ipiv[k], abs(i__1));
00417 if (kp != k - 1) {
00418 sswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1],
00419 ldb);
00420 }
00421
00422
00423
00424 i__1 = k - 2;
00425 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset],
00426 ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k +
00427 b_dim1], ldb);
00428 i__1 = k - 2;
00429 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset],
00430 ldb, &a[(k - 1) * a_dim1 + 1], &c__1, &c_b15, &b[
00431 k - 1 + b_dim1], ldb);
00432 }
00433
00434
00435
00436 if (nounit) {
00437 d11 = a[k - 1 + (k - 1) * a_dim1];
00438 d22 = a[k + k * a_dim1];
00439 d12 = a[k - 1 + k * a_dim1];
00440 d21 = d12;
00441 i__1 = *nrhs;
00442 for (j = 1; j <= i__1; ++j) {
00443 t1 = b[k - 1 + j * b_dim1];
00444 t2 = b[k + j * b_dim1];
00445 b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
00446 b[k + j * b_dim1] = d21 * t1 + d22 * t2;
00447
00448 }
00449 }
00450 k += -2;
00451 }
00452 goto L70;
00453 L90:
00454
00455
00456
00457
00458
00459 ;
00460 } else {
00461
00462
00463
00464 k = 1;
00465 L100:
00466 if (k > *n) {
00467 goto L120;
00468 }
00469
00470
00471
00472 if (ipiv[k] > 0) {
00473 if (k < *n) {
00474
00475
00476
00477 kp = ipiv[k];
00478 if (kp != k) {
00479 sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1],
00480 ldb);
00481 }
00482
00483
00484
00485 i__1 = *n - k;
00486 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 1 +
00487 b_dim1], ldb, &a[k + 1 + k * a_dim1], &c__1, &
00488 c_b15, &b[k + b_dim1], ldb);
00489 }
00490 if (nounit) {
00491 sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00492 }
00493 ++k;
00494
00495
00496
00497 } else {
00498 if (k < *n - 1) {
00499
00500
00501
00502 kp = (i__1 = ipiv[k], abs(i__1));
00503 if (kp != k + 1) {
00504 sswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1],
00505 ldb);
00506 }
00507
00508
00509
00510 i__1 = *n - k - 1;
00511 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 +
00512 b_dim1], ldb, &a[k + 2 + (k + 1) * a_dim1], &c__1,
00513 &c_b15, &b[k + 1 + b_dim1], ldb);
00514 i__1 = *n - k - 1;
00515 sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 +
00516 b_dim1], ldb, &a[k + 2 + k * a_dim1], &c__1, &
00517 c_b15, &b[k + b_dim1], ldb);
00518 }
00519
00520
00521
00522 if (nounit) {
00523 d11 = a[k + k * a_dim1];
00524 d22 = a[k + 1 + (k + 1) * a_dim1];
00525 d21 = a[k + 1 + k * a_dim1];
00526 d12 = d21;
00527 i__1 = *nrhs;
00528 for (j = 1; j <= i__1; ++j) {
00529 t1 = b[k + j * b_dim1];
00530 t2 = b[k + 1 + j * b_dim1];
00531 b[k + j * b_dim1] = d11 * t1 + d12 * t2;
00532 b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
00533
00534 }
00535 }
00536 k += 2;
00537 }
00538 goto L100;
00539 L120:
00540 ;
00541 }
00542
00543 }
00544 return 0;
00545
00546
00547
00548 }