00001
00002
00003
00004
00005 #include "f2c.h"
00006
00007 #ifdef HAVE_CONFIG
00008 #include "config.h"
00009 #else
00010 extern doublereal slamch_(char *);
00011 #define EPSILON slamch_("Epsilon")
00012 #define SAFEMINIMUM slamch_("Safe minimum")
00013 #define PRECISION slamch_("Precision")
00014 #define BASE slamch_("Base")
00015 #endif
00016
00017
00018 extern doublereal slapy2_(real *, real *);
00019
00020
00021
00022
00023
00024 static integer c__1 = 1;
00025
00026 logical lsame_(char *ca, char *cb)
00027 {
00028
00029 logical ret_val;
00030
00031
00032 static integer inta, intb, zcode;
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
00062 if (ret_val) {
00063 return ret_val;
00064 }
00065
00066
00067
00068 zcode = 'Z';
00069
00070
00071
00072
00073
00074
00075
00076
00077 inta = *(unsigned char *)ca;
00078 intb = *(unsigned char *)cb;
00079
00080 if (zcode == 90 || zcode == 122) {
00081
00082
00083
00084
00085
00086
00087 if (inta >= 97 && inta <= 122) {
00088 inta += -32;
00089 }
00090 if (intb >= 97 && intb <= 122) {
00091 intb += -32;
00092 }
00093
00094 } else if (zcode == 233 || zcode == 169) {
00095
00096
00097
00098
00099
00100
00101 if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta
00102 >= 162 && inta <= 169) {
00103 inta += 64;
00104 }
00105 if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb
00106 >= 162 && intb <= 169) {
00107 intb += 64;
00108 }
00109
00110 } else if (zcode == 218 || zcode == 250) {
00111
00112
00113
00114
00115
00116
00117 if (inta >= 225 && inta <= 250) {
00118 inta += -32;
00119 }
00120 if (intb >= 225 && intb <= 250) {
00121 intb += -32;
00122 }
00123 }
00124 ret_val = inta == intb;
00125
00126
00127
00128
00129
00130
00131
00132 return ret_val;
00133 }
00134
00135 doublereal sdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
00136 {
00137
00138 integer i__1;
00139 real ret_val;
00140
00141
00142 static integer i__, m, ix, iy, mp1;
00143 static real stemp;
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 --sy;
00156 --sx;
00157
00158
00159 stemp = 0.f;
00160 ret_val = 0.f;
00161 if (*n <= 0) {
00162 return ret_val;
00163 }
00164 if (*incx == 1 && *incy == 1) {
00165 goto L20;
00166 }
00167
00168
00169
00170
00171
00172
00173 ix = 1;
00174 iy = 1;
00175 if (*incx < 0) {
00176 ix = (-(*n) + 1) * *incx + 1;
00177 }
00178 if (*incy < 0) {
00179 iy = (-(*n) + 1) * *incy + 1;
00180 }
00181 i__1 = *n;
00182 for (i__ = 1; i__ <= i__1; ++i__) {
00183 stemp += sx[ix] * sy[iy];
00184 ix += *incx;
00185 iy += *incy;
00186
00187 }
00188 ret_val = stemp;
00189 return ret_val;
00190
00191
00192
00193
00194
00195
00196
00197
00198 L20:
00199 m = *n % 5;
00200 if (m == 0) {
00201 goto L40;
00202 }
00203 i__1 = m;
00204 for (i__ = 1; i__ <= i__1; ++i__) {
00205 stemp += sx[i__] * sy[i__];
00206
00207 }
00208 if (*n < 5) {
00209 goto L60;
00210 }
00211 L40:
00212 mp1 = m + 1;
00213 i__1 = *n;
00214 for (i__ = mp1; i__ <= i__1; i__ += 5) {
00215 stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[
00216 i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ +
00217 4] * sy[i__ + 4];
00218
00219 }
00220 L60:
00221 ret_val = stemp;
00222 return ret_val;
00223 }
00224
00225 int sgemm_(char *transa, char *transb, integer *m, integer *
00226 n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
00227 ldb, real *beta, real *c__, integer *ldc)
00228 {
00229
00230 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
00231 i__3;
00232
00233
00234 static integer i__, j, l, info;
00235 static logical nota, notb;
00236 static real temp;
00237 static integer ncola;
00238 extern logical lsame_(char *, char *);
00239 static integer nrowa, nrowb;
00240 extern int xerbla_(char *, integer *);
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370 a_dim1 = *lda;
00371 a_offset = 1 + a_dim1;
00372 a -= a_offset;
00373 b_dim1 = *ldb;
00374 b_offset = 1 + b_dim1;
00375 b -= b_offset;
00376 c_dim1 = *ldc;
00377 c_offset = 1 + c_dim1;
00378 c__ -= c_offset;
00379
00380
00381 nota = lsame_(transa, "N");
00382 notb = lsame_(transb, "N");
00383 if (nota) {
00384 nrowa = *m;
00385 ncola = *k;
00386 } else {
00387 nrowa = *k;
00388 ncola = *m;
00389 }
00390 if (notb) {
00391 nrowb = *k;
00392 } else {
00393 nrowb = *n;
00394 }
00395
00396
00397
00398 info = 0;
00399 if (! nota && ! lsame_(transa, "C") && ! lsame_(
00400 transa, "T")) {
00401 info = 1;
00402 } else if (! notb && ! lsame_(transb, "C") && !
00403 lsame_(transb, "T")) {
00404 info = 2;
00405 } else if (*m < 0) {
00406 info = 3;
00407 } else if (*n < 0) {
00408 info = 4;
00409 } else if (*k < 0) {
00410 info = 5;
00411 } else if (*lda < max(1,nrowa)) {
00412 info = 8;
00413 } else if (*ldb < max(1,nrowb)) {
00414 info = 10;
00415 } else if (*ldc < max(1,*m)) {
00416 info = 13;
00417 }
00418 if (info != 0) {
00419 xerbla_("SGEMM ", &info);
00420 return 0;
00421 }
00422
00423
00424
00425 if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
00426 return 0;
00427 }
00428
00429
00430
00431 if (*alpha == 0.f) {
00432 if (*beta == 0.f) {
00433 i__1 = *n;
00434 for (j = 1; j <= i__1; ++j) {
00435 i__2 = *m;
00436 for (i__ = 1; i__ <= i__2; ++i__) {
00437 c__[i__ + j * c_dim1] = 0.f;
00438
00439 }
00440
00441 }
00442 } else {
00443 i__1 = *n;
00444 for (j = 1; j <= i__1; ++j) {
00445 i__2 = *m;
00446 for (i__ = 1; i__ <= i__2; ++i__) {
00447 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00448
00449 }
00450
00451 }
00452 }
00453 return 0;
00454 }
00455
00456
00457
00458 if (notb) {
00459 if (nota) {
00460
00461
00462
00463 i__1 = *n;
00464 for (j = 1; j <= i__1; ++j) {
00465 if (*beta == 0.f) {
00466 i__2 = *m;
00467 for (i__ = 1; i__ <= i__2; ++i__) {
00468 c__[i__ + j * c_dim1] = 0.f;
00469
00470 }
00471 } else if (*beta != 1.f) {
00472 i__2 = *m;
00473 for (i__ = 1; i__ <= i__2; ++i__) {
00474 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00475
00476 }
00477 }
00478 i__2 = *k;
00479 for (l = 1; l <= i__2; ++l) {
00480 if (b[l + j * b_dim1] != 0.f) {
00481 temp = *alpha * b[l + j * b_dim1];
00482 i__3 = *m;
00483 for (i__ = 1; i__ <= i__3; ++i__) {
00484 c__[i__ + j * c_dim1] += temp * a[i__ + l *
00485 a_dim1];
00486
00487 }
00488 }
00489
00490 }
00491
00492 }
00493 } else {
00494
00495
00496
00497 i__1 = *n;
00498 for (j = 1; j <= i__1; ++j) {
00499 i__2 = *m;
00500 for (i__ = 1; i__ <= i__2; ++i__) {
00501 temp = 0.f;
00502 i__3 = *k;
00503 for (l = 1; l <= i__3; ++l) {
00504 temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
00505
00506 }
00507 if (*beta == 0.f) {
00508 c__[i__ + j * c_dim1] = *alpha * temp;
00509 } else {
00510 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
00511 i__ + j * c_dim1];
00512 }
00513
00514 }
00515
00516 }
00517 }
00518 } else {
00519 if (nota) {
00520
00521
00522
00523 i__1 = *n;
00524 for (j = 1; j <= i__1; ++j) {
00525 if (*beta == 0.f) {
00526 i__2 = *m;
00527 for (i__ = 1; i__ <= i__2; ++i__) {
00528 c__[i__ + j * c_dim1] = 0.f;
00529
00530 }
00531 } else if (*beta != 1.f) {
00532 i__2 = *m;
00533 for (i__ = 1; i__ <= i__2; ++i__) {
00534 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
00535
00536 }
00537 }
00538 i__2 = *k;
00539 for (l = 1; l <= i__2; ++l) {
00540 if (b[j + l * b_dim1] != 0.f) {
00541 temp = *alpha * b[j + l * b_dim1];
00542 i__3 = *m;
00543 for (i__ = 1; i__ <= i__3; ++i__) {
00544 c__[i__ + j * c_dim1] += temp * a[i__ + l *
00545 a_dim1];
00546
00547 }
00548 }
00549
00550 }
00551
00552 }
00553 } else {
00554
00555
00556
00557 i__1 = *n;
00558 for (j = 1; j <= i__1; ++j) {
00559 i__2 = *m;
00560 for (i__ = 1; i__ <= i__2; ++i__) {
00561 temp = 0.f;
00562 i__3 = *k;
00563 for (l = 1; l <= i__3; ++l) {
00564 temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
00565
00566 }
00567 if (*beta == 0.f) {
00568 c__[i__ + j * c_dim1] = *alpha * temp;
00569 } else {
00570 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
00571 i__ + j * c_dim1];
00572 }
00573
00574 }
00575
00576 }
00577 }
00578 }
00579
00580 return 0;
00581
00582
00583
00584 }
00585
00586 int sgemv_(char *trans, integer *m, integer *n, real *alpha,
00587 real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
00588 integer *incy)
00589 {
00590
00591 integer a_dim1, a_offset, i__1, i__2;
00592
00593
00594 static integer i__, j, ix, iy, jx, jy, kx, ky, info;
00595 static real temp;
00596 static integer lenx, leny;
00597 extern logical lsame_(char *, char *);
00598 extern int xerbla_(char *, integer *);
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697 a_dim1 = *lda;
00698 a_offset = 1 + a_dim1;
00699 a -= a_offset;
00700 --x;
00701 --y;
00702
00703
00704 info = 0;
00705 if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
00706 ) {
00707 info = 1;
00708 } else if (*m < 0) {
00709 info = 2;
00710 } else if (*n < 0) {
00711 info = 3;
00712 } else if (*lda < max(1,*m)) {
00713 info = 6;
00714 } else if (*incx == 0) {
00715 info = 8;
00716 } else if (*incy == 0) {
00717 info = 11;
00718 }
00719 if (info != 0) {
00720 xerbla_("SGEMV ", &info);
00721 return 0;
00722 }
00723
00724
00725
00726 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
00727 return 0;
00728 }
00729
00730
00731
00732
00733
00734
00735 if (lsame_(trans, "N")) {
00736 lenx = *n;
00737 leny = *m;
00738 } else {
00739 lenx = *m;
00740 leny = *n;
00741 }
00742 if (*incx > 0) {
00743 kx = 1;
00744 } else {
00745 kx = 1 - (lenx - 1) * *incx;
00746 }
00747 if (*incy > 0) {
00748 ky = 1;
00749 } else {
00750 ky = 1 - (leny - 1) * *incy;
00751 }
00752
00753
00754
00755
00756
00757
00758
00759
00760 if (*beta != 1.f) {
00761 if (*incy == 1) {
00762 if (*beta == 0.f) {
00763 i__1 = leny;
00764 for (i__ = 1; i__ <= i__1; ++i__) {
00765 y[i__] = 0.f;
00766
00767 }
00768 } else {
00769 i__1 = leny;
00770 for (i__ = 1; i__ <= i__1; ++i__) {
00771 y[i__] = *beta * y[i__];
00772
00773 }
00774 }
00775 } else {
00776 iy = ky;
00777 if (*beta == 0.f) {
00778 i__1 = leny;
00779 for (i__ = 1; i__ <= i__1; ++i__) {
00780 y[iy] = 0.f;
00781 iy += *incy;
00782
00783 }
00784 } else {
00785 i__1 = leny;
00786 for (i__ = 1; i__ <= i__1; ++i__) {
00787 y[iy] = *beta * y[iy];
00788 iy += *incy;
00789
00790 }
00791 }
00792 }
00793 }
00794 if (*alpha == 0.f) {
00795 return 0;
00796 }
00797 if (lsame_(trans, "N")) {
00798
00799
00800
00801 jx = kx;
00802 if (*incy == 1) {
00803 i__1 = *n;
00804 for (j = 1; j <= i__1; ++j) {
00805 if (x[jx] != 0.f) {
00806 temp = *alpha * x[jx];
00807 i__2 = *m;
00808 for (i__ = 1; i__ <= i__2; ++i__) {
00809 y[i__] += temp * a[i__ + j * a_dim1];
00810
00811 }
00812 }
00813 jx += *incx;
00814
00815 }
00816 } else {
00817 i__1 = *n;
00818 for (j = 1; j <= i__1; ++j) {
00819 if (x[jx] != 0.f) {
00820 temp = *alpha * x[jx];
00821 iy = ky;
00822 i__2 = *m;
00823 for (i__ = 1; i__ <= i__2; ++i__) {
00824 y[iy] += temp * a[i__ + j * a_dim1];
00825 iy += *incy;
00826
00827 }
00828 }
00829 jx += *incx;
00830
00831 }
00832 }
00833 } else {
00834
00835
00836
00837 jy = ky;
00838 if (*incx == 1) {
00839 i__1 = *n;
00840 for (j = 1; j <= i__1; ++j) {
00841 temp = 0.f;
00842 i__2 = *m;
00843 for (i__ = 1; i__ <= i__2; ++i__) {
00844 temp += a[i__ + j * a_dim1] * x[i__];
00845
00846 }
00847 y[jy] += *alpha * temp;
00848 jy += *incy;
00849
00850 }
00851 } else {
00852 i__1 = *n;
00853 for (j = 1; j <= i__1; ++j) {
00854 temp = 0.f;
00855 ix = kx;
00856 i__2 = *m;
00857 for (i__ = 1; i__ <= i__2; ++i__) {
00858 temp += a[i__ + j * a_dim1] * x[ix];
00859 ix += *incx;
00860
00861 }
00862 y[jy] += *alpha * temp;
00863 jy += *incy;
00864
00865 }
00866 }
00867 }
00868
00869 return 0;
00870
00871
00872
00873 }
00874
00875 int sscal_(integer *n, real *sa, real *sx, integer *incx)
00876 {
00877
00878 integer i__1, i__2;
00879
00880
00881 static integer i__, m, mp1, nincx;
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894 --sx;
00895
00896
00897 if (*n <= 0 || *incx <= 0) {
00898 return 0;
00899 }
00900 if (*incx == 1) {
00901 goto L20;
00902 }
00903
00904
00905
00906 nincx = *n * *incx;
00907 i__1 = nincx;
00908 i__2 = *incx;
00909 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
00910 sx[i__] = *sa * sx[i__];
00911
00912 }
00913 return 0;
00914
00915
00916
00917
00918
00919
00920
00921
00922 L20:
00923 m = *n % 5;
00924 if (m == 0) {
00925 goto L40;
00926 }
00927 i__2 = m;
00928 for (i__ = 1; i__ <= i__2; ++i__) {
00929 sx[i__] = *sa * sx[i__];
00930
00931 }
00932 if (*n < 5) {
00933 return 0;
00934 }
00935 L40:
00936 mp1 = m + 1;
00937 i__2 = *n;
00938 for (i__ = mp1; i__ <= i__2; i__ += 5) {
00939 sx[i__] = *sa * sx[i__];
00940 sx[i__ + 1] = *sa * sx[i__ + 1];
00941 sx[i__ + 2] = *sa * sx[i__ + 2];
00942 sx[i__ + 3] = *sa * sx[i__ + 3];
00943 sx[i__ + 4] = *sa * sx[i__ + 4];
00944
00945 }
00946 return 0;
00947 }
00948
00949 int ssymm_(char *side, char *uplo, integer *m, integer *n,
00950 real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta,
00951 real *c__, integer *ldc)
00952 {
00953
00954 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
00955 i__3;
00956
00957
00958 static integer i__, j, k, info;
00959 static real temp1, temp2;
00960 extern logical lsame_(char *, char *);
00961 static integer nrowa;
00962 static logical upper;
00963 extern int xerbla_(char *, integer *);
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095 a_dim1 = *lda;
01096 a_offset = 1 + a_dim1;
01097 a -= a_offset;
01098 b_dim1 = *ldb;
01099 b_offset = 1 + b_dim1;
01100 b -= b_offset;
01101 c_dim1 = *ldc;
01102 c_offset = 1 + c_dim1;
01103 c__ -= c_offset;
01104
01105
01106 if (lsame_(side, "L")) {
01107 nrowa = *m;
01108 } else {
01109 nrowa = *n;
01110 }
01111 upper = lsame_(uplo, "U");
01112
01113
01114
01115 info = 0;
01116 if (! lsame_(side, "L") && ! lsame_(side, "R")) {
01117 info = 1;
01118 } else if (! upper && ! lsame_(uplo, "L")) {
01119 info = 2;
01120 } else if (*m < 0) {
01121 info = 3;
01122 } else if (*n < 0) {
01123 info = 4;
01124 } else if (*lda < max(1,nrowa)) {
01125 info = 7;
01126 } else if (*ldb < max(1,*m)) {
01127 info = 9;
01128 } else if (*ldc < max(1,*m)) {
01129 info = 12;
01130 }
01131 if (info != 0) {
01132 xerbla_("SSYMM ", &info);
01133 return 0;
01134 }
01135
01136
01137
01138 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
01139 return 0;
01140 }
01141
01142
01143
01144 if (*alpha == 0.f) {
01145 if (*beta == 0.f) {
01146 i__1 = *n;
01147 for (j = 1; j <= i__1; ++j) {
01148 i__2 = *m;
01149 for (i__ = 1; i__ <= i__2; ++i__) {
01150 c__[i__ + j * c_dim1] = 0.f;
01151
01152 }
01153
01154 }
01155 } else {
01156 i__1 = *n;
01157 for (j = 1; j <= i__1; ++j) {
01158 i__2 = *m;
01159 for (i__ = 1; i__ <= i__2; ++i__) {
01160 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01161
01162 }
01163
01164 }
01165 }
01166 return 0;
01167 }
01168
01169
01170
01171 if (lsame_(side, "L")) {
01172
01173
01174
01175 if (upper) {
01176 i__1 = *n;
01177 for (j = 1; j <= i__1; ++j) {
01178 i__2 = *m;
01179 for (i__ = 1; i__ <= i__2; ++i__) {
01180 temp1 = *alpha * b[i__ + j * b_dim1];
01181 temp2 = 0.f;
01182 i__3 = i__ - 1;
01183 for (k = 1; k <= i__3; ++k) {
01184 c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
01185 temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
01186
01187 }
01188 if (*beta == 0.f) {
01189 c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
01190 + *alpha * temp2;
01191 } else {
01192 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
01193 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
01194 temp2;
01195 }
01196
01197 }
01198
01199 }
01200 } else {
01201 i__1 = *n;
01202 for (j = 1; j <= i__1; ++j) {
01203 for (i__ = *m; i__ >= 1; --i__) {
01204 temp1 = *alpha * b[i__ + j * b_dim1];
01205 temp2 = 0.f;
01206 i__2 = *m;
01207 for (k = i__ + 1; k <= i__2; ++k) {
01208 c__[k + j * c_dim1] += temp1 * a[k + i__ * a_dim1];
01209 temp2 += b[k + j * b_dim1] * a[k + i__ * a_dim1];
01210
01211 }
01212 if (*beta == 0.f) {
01213 c__[i__ + j * c_dim1] = temp1 * a[i__ + i__ * a_dim1]
01214 + *alpha * temp2;
01215 } else {
01216 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]
01217 + temp1 * a[i__ + i__ * a_dim1] + *alpha *
01218 temp2;
01219 }
01220
01221 }
01222
01223 }
01224 }
01225 } else {
01226
01227
01228
01229 i__1 = *n;
01230 for (j = 1; j <= i__1; ++j) {
01231 temp1 = *alpha * a[j + j * a_dim1];
01232 if (*beta == 0.f) {
01233 i__2 = *m;
01234 for (i__ = 1; i__ <= i__2; ++i__) {
01235 c__[i__ + j * c_dim1] = temp1 * b[i__ + j * b_dim1];
01236
01237 }
01238 } else {
01239 i__2 = *m;
01240 for (i__ = 1; i__ <= i__2; ++i__) {
01241 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1] +
01242 temp1 * b[i__ + j * b_dim1];
01243
01244 }
01245 }
01246 i__2 = j - 1;
01247 for (k = 1; k <= i__2; ++k) {
01248 if (upper) {
01249 temp1 = *alpha * a[k + j * a_dim1];
01250 } else {
01251 temp1 = *alpha * a[j + k * a_dim1];
01252 }
01253 i__3 = *m;
01254 for (i__ = 1; i__ <= i__3; ++i__) {
01255 c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
01256
01257 }
01258
01259 }
01260 i__2 = *n;
01261 for (k = j + 1; k <= i__2; ++k) {
01262 if (upper) {
01263 temp1 = *alpha * a[j + k * a_dim1];
01264 } else {
01265 temp1 = *alpha * a[k + j * a_dim1];
01266 }
01267 i__3 = *m;
01268 for (i__ = 1; i__ <= i__3; ++i__) {
01269 c__[i__ + j * c_dim1] += temp1 * b[i__ + k * b_dim1];
01270
01271 }
01272
01273 }
01274
01275 }
01276 }
01277
01278 return 0;
01279
01280
01281
01282 }
01283
01284 int ssyrk_(char *uplo, char *trans, integer *n, integer *k,
01285 real *alpha, real *a, integer *lda, real *beta, real *c__, integer *
01286 ldc)
01287 {
01288
01289 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
01290
01291
01292 static integer i__, j, l, info;
01293 static real temp;
01294 extern logical lsame_(char *, char *);
01295 static integer nrowa;
01296 static logical upper;
01297 extern int xerbla_(char *, integer *);
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413 a_dim1 = *lda;
01414 a_offset = 1 + a_dim1;
01415 a -= a_offset;
01416 c_dim1 = *ldc;
01417 c_offset = 1 + c_dim1;
01418 c__ -= c_offset;
01419
01420
01421 if (lsame_(trans, "N")) {
01422 nrowa = *n;
01423 } else {
01424 nrowa = *k;
01425 }
01426 upper = lsame_(uplo, "U");
01427
01428 info = 0;
01429 if (! upper && ! lsame_(uplo, "L")) {
01430 info = 1;
01431 } else if (! lsame_(trans, "N") && ! lsame_(trans,
01432 "T") && ! lsame_(trans, "C")) {
01433 info = 2;
01434 } else if (*n < 0) {
01435 info = 3;
01436 } else if (*k < 0) {
01437 info = 4;
01438 } else if (*lda < max(1,nrowa)) {
01439 info = 7;
01440 } else if (*ldc < max(1,*n)) {
01441 info = 10;
01442 }
01443 if (info != 0) {
01444 xerbla_("SSYRK ", &info);
01445 return 0;
01446 }
01447
01448
01449
01450 if (*n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
01451 return 0;
01452 }
01453
01454
01455
01456 if (*alpha == 0.f) {
01457 if (upper) {
01458 if (*beta == 0.f) {
01459 i__1 = *n;
01460 for (j = 1; j <= i__1; ++j) {
01461 i__2 = j;
01462 for (i__ = 1; i__ <= i__2; ++i__) {
01463 c__[i__ + j * c_dim1] = 0.f;
01464
01465 }
01466
01467 }
01468 } else {
01469 i__1 = *n;
01470 for (j = 1; j <= i__1; ++j) {
01471 i__2 = j;
01472 for (i__ = 1; i__ <= i__2; ++i__) {
01473 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01474
01475 }
01476
01477 }
01478 }
01479 } else {
01480 if (*beta == 0.f) {
01481 i__1 = *n;
01482 for (j = 1; j <= i__1; ++j) {
01483 i__2 = *n;
01484 for (i__ = j; i__ <= i__2; ++i__) {
01485 c__[i__ + j * c_dim1] = 0.f;
01486
01487 }
01488
01489 }
01490 } else {
01491 i__1 = *n;
01492 for (j = 1; j <= i__1; ++j) {
01493 i__2 = *n;
01494 for (i__ = j; i__ <= i__2; ++i__) {
01495 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01496
01497 }
01498
01499 }
01500 }
01501 }
01502 return 0;
01503 }
01504
01505
01506
01507 if (lsame_(trans, "N")) {
01508
01509
01510
01511 if (upper) {
01512 i__1 = *n;
01513 for (j = 1; j <= i__1; ++j) {
01514 if (*beta == 0.f) {
01515 i__2 = j;
01516 for (i__ = 1; i__ <= i__2; ++i__) {
01517 c__[i__ + j * c_dim1] = 0.f;
01518
01519 }
01520 } else if (*beta != 1.f) {
01521 i__2 = j;
01522 for (i__ = 1; i__ <= i__2; ++i__) {
01523 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01524
01525 }
01526 }
01527 i__2 = *k;
01528 for (l = 1; l <= i__2; ++l) {
01529 if (a[j + l * a_dim1] != 0.f) {
01530 temp = *alpha * a[j + l * a_dim1];
01531 i__3 = j;
01532 for (i__ = 1; i__ <= i__3; ++i__) {
01533 c__[i__ + j * c_dim1] += temp * a[i__ + l *
01534 a_dim1];
01535
01536 }
01537 }
01538
01539 }
01540
01541 }
01542 } else {
01543 i__1 = *n;
01544 for (j = 1; j <= i__1; ++j) {
01545 if (*beta == 0.f) {
01546 i__2 = *n;
01547 for (i__ = j; i__ <= i__2; ++i__) {
01548 c__[i__ + j * c_dim1] = 0.f;
01549
01550 }
01551 } else if (*beta != 1.f) {
01552 i__2 = *n;
01553 for (i__ = j; i__ <= i__2; ++i__) {
01554 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
01555
01556 }
01557 }
01558 i__2 = *k;
01559 for (l = 1; l <= i__2; ++l) {
01560 if (a[j + l * a_dim1] != 0.f) {
01561 temp = *alpha * a[j + l * a_dim1];
01562 i__3 = *n;
01563 for (i__ = j; i__ <= i__3; ++i__) {
01564 c__[i__ + j * c_dim1] += temp * a[i__ + l *
01565 a_dim1];
01566
01567 }
01568 }
01569
01570 }
01571
01572 }
01573 }
01574 } else {
01575
01576
01577
01578 if (upper) {
01579 i__1 = *n;
01580 for (j = 1; j <= i__1; ++j) {
01581 i__2 = j;
01582 for (i__ = 1; i__ <= i__2; ++i__) {
01583 temp = 0.f;
01584 i__3 = *k;
01585 for (l = 1; l <= i__3; ++l) {
01586 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
01587
01588 }
01589 if (*beta == 0.f) {
01590 c__[i__ + j * c_dim1] = *alpha * temp;
01591 } else {
01592 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
01593 i__ + j * c_dim1];
01594 }
01595
01596 }
01597
01598 }
01599 } else {
01600 i__1 = *n;
01601 for (j = 1; j <= i__1; ++j) {
01602 i__2 = *n;
01603 for (i__ = j; i__ <= i__2; ++i__) {
01604 temp = 0.f;
01605 i__3 = *k;
01606 for (l = 1; l <= i__3; ++l) {
01607 temp += a[l + i__ * a_dim1] * a[l + j * a_dim1];
01608
01609 }
01610 if (*beta == 0.f) {
01611 c__[i__ + j * c_dim1] = *alpha * temp;
01612 } else {
01613 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
01614 i__ + j * c_dim1];
01615 }
01616
01617 }
01618
01619 }
01620 }
01621 }
01622
01623 return 0;
01624
01625
01626
01627 }
01628
01629 int strsm_(char *side, char *uplo, char *transa, char *diag,
01630 integer *m, integer *n, real *alpha, real *a, integer *lda, real *b,
01631 integer *ldb)
01632 {
01633
01634 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
01635
01636
01637 static integer i__, j, k, info;
01638 static real temp;
01639 static logical lside;
01640 extern logical lsame_(char *, char *);
01641 static integer nrowa;
01642 static logical upper;
01643 extern int xerbla_(char *, integer *);
01644 static logical nounit;
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771 a_dim1 = *lda;
01772 a_offset = 1 + a_dim1;
01773 a -= a_offset;
01774 b_dim1 = *ldb;
01775 b_offset = 1 + b_dim1;
01776 b -= b_offset;
01777
01778
01779 lside = lsame_(side, "L");
01780 if (lside) {
01781 nrowa = *m;
01782 } else {
01783 nrowa = *n;
01784 }
01785 nounit = lsame_(diag, "N");
01786 upper = lsame_(uplo, "U");
01787
01788 info = 0;
01789 if (! lside && ! lsame_(side, "R")) {
01790 info = 1;
01791 } else if (! upper && ! lsame_(uplo, "L")) {
01792 info = 2;
01793 } else if (! lsame_(transa, "N") && ! lsame_(transa,
01794 "T") && ! lsame_(transa, "C")) {
01795 info = 3;
01796 } else if (! lsame_(diag, "U") && ! lsame_(diag,
01797 "N")) {
01798 info = 4;
01799 } else if (*m < 0) {
01800 info = 5;
01801 } else if (*n < 0) {
01802 info = 6;
01803 } else if (*lda < max(1,nrowa)) {
01804 info = 9;
01805 } else if (*ldb < max(1,*m)) {
01806 info = 11;
01807 }
01808 if (info != 0) {
01809 xerbla_("STRSM ", &info);
01810 return 0;
01811 }
01812
01813
01814
01815 if (*n == 0) {
01816 return 0;
01817 }
01818
01819
01820
01821 if (*alpha == 0.f) {
01822 i__1 = *n;
01823 for (j = 1; j <= i__1; ++j) {
01824 i__2 = *m;
01825 for (i__ = 1; i__ <= i__2; ++i__) {
01826 b[i__ + j * b_dim1] = 0.f;
01827
01828 }
01829
01830 }
01831 return 0;
01832 }
01833
01834
01835
01836 if (lside) {
01837 if (lsame_(transa, "N")) {
01838
01839
01840
01841 if (upper) {
01842 i__1 = *n;
01843 for (j = 1; j <= i__1; ++j) {
01844 if (*alpha != 1.f) {
01845 i__2 = *m;
01846 for (i__ = 1; i__ <= i__2; ++i__) {
01847 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01848 ;
01849
01850 }
01851 }
01852 for (k = *m; k >= 1; --k) {
01853 if (b[k + j * b_dim1] != 0.f) {
01854 if (nounit) {
01855 b[k + j * b_dim1] /= a[k + k * a_dim1];
01856 }
01857 i__2 = k - 1;
01858 for (i__ = 1; i__ <= i__2; ++i__) {
01859 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
01860 i__ + k * a_dim1];
01861
01862 }
01863 }
01864
01865 }
01866
01867 }
01868 } else {
01869 i__1 = *n;
01870 for (j = 1; j <= i__1; ++j) {
01871 if (*alpha != 1.f) {
01872 i__2 = *m;
01873 for (i__ = 1; i__ <= i__2; ++i__) {
01874 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01875 ;
01876
01877 }
01878 }
01879 i__2 = *m;
01880 for (k = 1; k <= i__2; ++k) {
01881 if (b[k + j * b_dim1] != 0.f) {
01882 if (nounit) {
01883 b[k + j * b_dim1] /= a[k + k * a_dim1];
01884 }
01885 i__3 = *m;
01886 for (i__ = k + 1; i__ <= i__3; ++i__) {
01887 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
01888 i__ + k * a_dim1];
01889
01890 }
01891 }
01892
01893 }
01894
01895 }
01896 }
01897 } else {
01898
01899
01900
01901 if (upper) {
01902 i__1 = *n;
01903 for (j = 1; j <= i__1; ++j) {
01904 i__2 = *m;
01905 for (i__ = 1; i__ <= i__2; ++i__) {
01906 temp = *alpha * b[i__ + j * b_dim1];
01907 i__3 = i__ - 1;
01908 for (k = 1; k <= i__3; ++k) {
01909 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
01910
01911 }
01912 if (nounit) {
01913 temp /= a[i__ + i__ * a_dim1];
01914 }
01915 b[i__ + j * b_dim1] = temp;
01916
01917 }
01918
01919 }
01920 } else {
01921 i__1 = *n;
01922 for (j = 1; j <= i__1; ++j) {
01923 for (i__ = *m; i__ >= 1; --i__) {
01924 temp = *alpha * b[i__ + j * b_dim1];
01925 i__2 = *m;
01926 for (k = i__ + 1; k <= i__2; ++k) {
01927 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
01928
01929 }
01930 if (nounit) {
01931 temp /= a[i__ + i__ * a_dim1];
01932 }
01933 b[i__ + j * b_dim1] = temp;
01934
01935 }
01936
01937 }
01938 }
01939 }
01940 } else {
01941 if (lsame_(transa, "N")) {
01942
01943
01944
01945 if (upper) {
01946 i__1 = *n;
01947 for (j = 1; j <= i__1; ++j) {
01948 if (*alpha != 1.f) {
01949 i__2 = *m;
01950 for (i__ = 1; i__ <= i__2; ++i__) {
01951 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01952 ;
01953
01954 }
01955 }
01956 i__2 = j - 1;
01957 for (k = 1; k <= i__2; ++k) {
01958 if (a[k + j * a_dim1] != 0.f) {
01959 i__3 = *m;
01960 for (i__ = 1; i__ <= i__3; ++i__) {
01961 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
01962 i__ + k * b_dim1];
01963
01964 }
01965 }
01966
01967 }
01968 if (nounit) {
01969 temp = 1.f / a[j + j * a_dim1];
01970 i__2 = *m;
01971 for (i__ = 1; i__ <= i__2; ++i__) {
01972 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
01973
01974 }
01975 }
01976
01977 }
01978 } else {
01979 for (j = *n; j >= 1; --j) {
01980 if (*alpha != 1.f) {
01981 i__1 = *m;
01982 for (i__ = 1; i__ <= i__1; ++i__) {
01983 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
01984 ;
01985
01986 }
01987 }
01988 i__1 = *n;
01989 for (k = j + 1; k <= i__1; ++k) {
01990 if (a[k + j * a_dim1] != 0.f) {
01991 i__2 = *m;
01992 for (i__ = 1; i__ <= i__2; ++i__) {
01993 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
01994 i__ + k * b_dim1];
01995
01996 }
01997 }
01998
01999 }
02000 if (nounit) {
02001 temp = 1.f / a[j + j * a_dim1];
02002 i__1 = *m;
02003 for (i__ = 1; i__ <= i__1; ++i__) {
02004 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
02005
02006 }
02007 }
02008
02009 }
02010 }
02011 } else {
02012
02013
02014
02015 if (upper) {
02016 for (k = *n; k >= 1; --k) {
02017 if (nounit) {
02018 temp = 1.f / a[k + k * a_dim1];
02019 i__1 = *m;
02020 for (i__ = 1; i__ <= i__1; ++i__) {
02021 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
02022
02023 }
02024 }
02025 i__1 = k - 1;
02026 for (j = 1; j <= i__1; ++j) {
02027 if (a[j + k * a_dim1] != 0.f) {
02028 temp = a[j + k * a_dim1];
02029 i__2 = *m;
02030 for (i__ = 1; i__ <= i__2; ++i__) {
02031 b[i__ + j * b_dim1] -= temp * b[i__ + k *
02032 b_dim1];
02033
02034 }
02035 }
02036
02037 }
02038 if (*alpha != 1.f) {
02039 i__1 = *m;
02040 for (i__ = 1; i__ <= i__1; ++i__) {
02041 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
02042 ;
02043
02044 }
02045 }
02046
02047 }
02048 } else {
02049 i__1 = *n;
02050 for (k = 1; k <= i__1; ++k) {
02051 if (nounit) {
02052 temp = 1.f / a[k + k * a_dim1];
02053 i__2 = *m;
02054 for (i__ = 1; i__ <= i__2; ++i__) {
02055 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
02056
02057 }
02058 }
02059 i__2 = *n;
02060 for (j = k + 1; j <= i__2; ++j) {
02061 if (a[j + k * a_dim1] != 0.f) {
02062 temp = a[j + k * a_dim1];
02063 i__3 = *m;
02064 for (i__ = 1; i__ <= i__3; ++i__) {
02065 b[i__ + j * b_dim1] -= temp * b[i__ + k *
02066 b_dim1];
02067
02068 }
02069 }
02070
02071 }
02072 if (*alpha != 1.f) {
02073 i__2 = *m;
02074 for (i__ = 1; i__ <= i__2; ++i__) {
02075 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
02076 ;
02077
02078 }
02079 }
02080
02081 }
02082 }
02083 }
02084 }
02085
02086 return 0;
02087
02088
02089
02090 }
02091
02092 int xerbla_(char *srname, integer *info)
02093 {
02094
02095 static char fmt_9999[] = "(\002 ** On entry to \002,a6,\002 parameter nu"
02096 "mber \002,i2,\002 had \002,\002an illegal value\002)";
02097
02098
02099 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
02100 int s_stop(char *, ftnlen);
02101
02102
02103 static cilist io___60 = { 0, 6, 0, fmt_9999, 0 };
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135 s_wsfe(&io___60);
02136 do_fio(&c__1, srname, (ftnlen)6);
02137 do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer));
02138 e_wsfe();
02139
02140 s_stop("", (ftnlen)0);
02141
02142
02143
02144
02145 return 0;
02146 }
02147