00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "cddefines.h"
00012 #include "thirdparty.h"
00013
00014
00015
00016 #ifdef LAPACK
00017
00018
00019
00020 extern "C" void dgetrf_(int32 *M, int32 *N, double *A, int32 *LDA, int32 *IPIV, int32 *INFO);
00021 extern "C" void dgetrs_(char *TRANS, int32 *N, int32 *NRHS, double *A, int32 *LDA, int32 *iPiv, double *B,
00022 int32 *LDB, int32 *INFO, int32 translen);
00023 extern "C" void dgtsv_(int32 *n, int32 *nrhs, double *dl, double *d, double *du, double *b, int32 *ldb, int32 *info);
00024
00025 #else
00026
00027
00028
00029
00030
00031
00032
00033
00034 STATIC void DGETRF(int32,int32,double*,int32,int32[],int32*);
00035
00036
00037 STATIC void DGETRS(int32 TRANS,int32 N,int32 NRHS,double *A,int32 LDA,int32 IPIV[],double *B,int32 LDB,int32 *INFO);
00038
00039
00040
00041
00042 #endif
00043
00044
00045
00046
00047 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
00048 {
00049 if( *info == 0 )
00050 {
00051 int32 M_loc, N_loc, lda_loc;
00052
00053 ASSERT( M < INT32_MAX && N < INT32_MAX && lda < INT32_MAX );
00054
00055 M_loc = (int32)M;
00056 N_loc = (int32)N;
00057 lda_loc = (int32)lda;
00058
00059 # ifdef LAPACK
00060
00061 dgetrf_(&M_loc, &N_loc, A , &lda_loc, ipiv, info);
00062 # else
00063
00064 DGETRF(M_loc, N_loc, A, lda_loc, ipiv, info);
00065 # endif
00066 }
00067 }
00068
00069 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv,
00070 double *B, long ldb, int32 *info)
00071 {
00072 if( *info == 0 )
00073 {
00074 int32 N_loc, nrhs_loc, lda_loc, ldb_loc;
00075
00076 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && lda < INT32_MAX && ldb < INT32_MAX );
00077
00078 N_loc = (int32)N;
00079 nrhs_loc = (int32)nrhs;
00080 lda_loc = (int32)lda;
00081 ldb_loc = (int32)ldb;
00082
00083 # ifdef LAPACK
00084
00085 dgetrs_(&trans, &N_loc, &nrhs_loc, A, &lda_loc, ipiv, B, &ldb_loc, info, sizeof(char));
00086 # else
00087
00088 DGETRS(trans, N_loc, nrhs_loc, A, lda_loc, ipiv, B, ldb_loc, info);
00089 # endif
00090 }
00091 }
00092
00093 #if 0
00094 void dgtsv_wrapper(long N, long nrhs, double *dl, double *d__, double *du, double *b, long ldb, int32 *info)
00095 {
00096 printf("Inside dgtsv\n");
00097 cdEXIT( EXIT_FAILURE );
00098 if( *info == 0 )
00099 {
00100 int32 N_loc, nrhs_loc, ldb_loc;
00101
00102 ASSERT( N < INT32_MAX && nrhs < INT32_MAX && ldb < INT32_MAX );
00103
00104 N_loc = (int32)N;
00105 nrhs_loc = (int32)nrhs;
00106 ldb_loc = (int32)ldb;
00107
00108 # ifdef LAPACK
00109
00110 dgtsv_(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
00111 # else
00112
00113
00114 (void)DGTSV(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
00115 # endif
00116 }
00117 }
00118 #endif
00119
00120
00121 #ifndef LAPACK
00122
00123 #define ONE 1.0e0
00124 #define ZERO 0.0e0
00125
00126 #ifdef AA
00127 # undef AA
00128 #endif
00129 #ifdef BB
00130 # undef BB
00131 #endif
00132 #ifdef CC
00133 # undef CC
00134 #endif
00135
00136 #define AA(I_,J_) (*(A+(I_)*(LDA)+(J_)))
00137 #define BB(I_,J_) (*(B+(I_)*(LDB)+(J_)))
00138 #define CC(I_,J_) (*(C+(I_)*(LDC)+(J_)))
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 STATIC void DGEMM(int32 TRANSA,
00165 int32 TRANSB,
00166 int32 M,
00167 int32 N,
00168 int32 K,
00169 double ALPHA,
00170 double *A,
00171 int32 LDA,
00172 double *B,
00173 int32 LDB,
00174 double BETA,
00175 double *C,
00176 int32 LDC);
00177
00178
00179 STATIC int32 LSAME(int32 CA,
00180 int32 CB);
00181
00182
00183 STATIC int32 IDAMAX(int32 n,
00184 double dx[],
00185 int32 incx);
00186
00187
00188 STATIC void DTRSM(int32 SIDE,
00189 int32 UPLO,
00190 int32 TRANSA,
00191 int32 DIAG,
00192 int32 M,
00193 int32 N,
00194 double ALPHA,
00195 double *A,
00196 int32 LDA,
00197 double *B,
00198 int32 LDB);
00199
00200
00201 STATIC int32 ILAENV(int32 ISPEC,
00202 const char *NAME,
00203
00204 int32 N1,
00205 int32 N2,
00206
00207 int32 N4);
00208
00209
00210 STATIC void DSWAP(int32 n,
00211 double dx[],
00212 int32 incx,
00213 double dy[],
00214 int32 incy);
00215
00216
00217 STATIC void DSCAL(int32 n,
00218 double da,
00219 double dx[],
00220 int32 incx);
00221
00222
00223 STATIC void DLASWP(int32 N,
00224 double *A,
00225 int32 LDA,
00226 int32 K1,
00227 int32 K2,
00228 int32 IPIV[],
00229 int32 INCX);
00230
00231
00232 STATIC void DGETF2(int32 M,
00233 int32 N,
00234 double *A,
00235 int32 LDA,
00236 int32 IPIV[],
00237 int32 *INFO);
00238
00239
00240 STATIC void DGER(int32 M,
00241 int32 N,
00242 double ALPHA,
00243 double X[],
00244 int32 INCX,
00245 double Y[],
00246 int32 INCY,
00247 double *A,
00248 int32 LDA);
00249
00250
00251 STATIC void XERBLA(const char *SRNAME,
00252 int32 INFO)
00253 {
00254
00255 DEBUG_ENTRY( "XERBLA()" );
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 fprintf( ioQQQ, " ** On entry to %6.6s parameter number %2ld had an illegal value\n",
00290 SRNAME, (long)INFO );
00291
00292 cdEXIT(EXIT_FAILURE);
00293 }
00294
00295
00296 STATIC void DGETRF(
00297
00298 int32 M,
00299
00300
00301 int32 N,
00302
00303 double *A,
00304
00305 int32 LDA,
00306
00307 int32 IPIV[],
00308
00309 int32 *INFO)
00310 {
00311
00312 char chL1,
00313 chL2,
00314 chL3,
00315 chL4;
00316 int32 I,
00317 IINFO,
00318 I_,
00319 J,
00320 JB,
00321 J_,
00322 NB,
00323
00324 limit,
00325 limit2;
00326
00327
00328 DEBUG_ENTRY( "DGETRF()" );
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
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394 *INFO = 0;
00395 if( M < 0 )
00396 {
00397 *INFO = -1;
00398 }
00399 else if( N < 0 )
00400 {
00401 *INFO = -2;
00402 }
00403 else if( LDA < MAX2(1,M) )
00404 {
00405 *INFO = -4;
00406 }
00407 if( *INFO != 0 )
00408 {
00409 XERBLA("DGETRF",-*INFO);
00410
00411 }
00412
00413
00414
00415 if( M == 0 || N == 0 )
00416 {
00417 return;
00418 }
00419
00420
00421
00422
00423 NB = ILAENV(1,"DGETRF",M,N,-1);
00424 if( NB <= 1 || NB >= MIN2(M,N) )
00425 {
00426
00427
00428 DGETF2(M,N,A,LDA,IPIV,INFO);
00429 }
00430 else
00431 {
00432
00433
00434
00435 limit = MIN2(M,N);
00436
00437
00438 for( J=1; J<=limit; J += NB )
00439 {
00440 J_ = J - 1;
00441 JB = MIN2(MIN2(M,N)-J+1,NB);
00442
00443
00444
00445
00446 DGETF2(M-J+1,JB,&AA(J_,J_),LDA,&IPIV[J_],&IINFO);
00447
00448
00449
00450 if( *INFO == 0 && IINFO > 0 )
00451 *INFO = IINFO + J - 1;
00452 limit2 = MIN2(M,J+JB-1);
00453 for( I=J; I <= limit2; I++ )
00454 {
00455 I_ = I - 1;
00456 IPIV[I_] += J - 1;
00457 }
00458
00459
00460
00461 DLASWP(J-1,A,LDA,J,J+JB-1,IPIV,1);
00462
00463 if( J + JB <= N )
00464 {
00465
00466
00467
00468 DLASWP(N-J-JB+1,&AA(J_+JB,0),LDA,J,J+JB-1,IPIV,1);
00469
00470
00471
00472 chL1 = 'L';
00473 chL2 = 'L';
00474 chL3 = 'N';
00475 chL4 = 'U';
00476
00477 DTRSM(chL1,chL2,chL3,chL4,JB,N-J-JB+1,ONE,&AA(J_,J_),
00478 LDA,&AA(J_+JB,J_),LDA);
00479 if( J + JB <= M )
00480 {
00481
00482
00483
00484 chL1 = 'N';
00485 chL2 = 'N';
00486
00487 DGEMM(chL1,chL2,M-J-JB+1,N-J-JB+1,JB,-ONE,&AA(J_,J_+JB),
00488 LDA,&AA(J_+JB,J_),LDA,ONE,&AA(J_+JB,J_+JB),LDA);
00489 }
00490 }
00491 }
00492 }
00493 return;
00494
00495
00496
00497 #undef A
00498 }
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 STATIC void DGETRS(
00528
00529 int32 TRANS,
00530
00531 int32 N,
00532
00533 int32 NRHS,
00534
00535 double *A,
00536
00537 int32 LDA,
00538
00539 int32 IPIV[],
00540
00541 double *B,
00542
00543 int32 LDB,
00544
00545 int32 *INFO)
00546 {
00547
00548
00549 int32 NOTRAN;
00550 char chL1,
00551 chL2,
00552 chL3,
00553 chL4;
00554
00555 DEBUG_ENTRY( "DGETRS()" );
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
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 *INFO = 0;
00626 NOTRAN = LSAME(TRANS,'N');
00627 if( (!NOTRAN && !LSAME(TRANS,'T')) && !LSAME(TRANS,'C') )
00628 {
00629 *INFO = -1;
00630 }
00631 else if( N < 0 )
00632 {
00633 *INFO = -2;
00634 }
00635 else if( NRHS < 0 )
00636 {
00637 *INFO = -3;
00638 }
00639 else if( LDA < MAX2(1,N) )
00640 {
00641 *INFO = -5;
00642 }
00643 else if( LDB < MAX2(1,N) )
00644 {
00645 *INFO = -8;
00646 }
00647 if( *INFO != 0 )
00648 {
00649 XERBLA("DGETRS",-*INFO);
00650
00651 }
00652
00653
00654
00655 if( N == 0 || NRHS == 0 )
00656 {
00657 return;
00658 }
00659
00660 if( NOTRAN )
00661 {
00662
00663
00664
00665
00666
00667 DLASWP(NRHS,B,LDB,1,N,IPIV,1);
00668
00669
00670
00671 chL1 = 'L';
00672 chL2 = 'L';
00673 chL3 = 'N';
00674 chL4 = 'U';
00675
00676 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00677
00678
00679
00680 chL1 = 'L';
00681 chL2 = 'U';
00682 chL3 = 'N';
00683 chL4 = 'N';
00684
00685 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00686 }
00687 else
00688 {
00689
00690
00691
00692
00693
00694 chL1 = 'L';
00695 chL2 = 'U';
00696 chL3 = 'T';
00697 chL4 = 'N';
00698
00699 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00700
00701
00702
00703 chL1 = 'L';
00704 chL2 = 'L';
00705 chL3 = 'T';
00706 chL4 = 'U';
00707
00708 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,ONE,A,LDA,B,LDB);
00709
00710
00711
00712 DLASWP(NRHS,B,LDB,1,N,IPIV,-1);
00713 }
00714
00715 return;
00716
00717
00718
00719 #undef B
00720 #undef A
00721 }
00722
00723
00724
00725 STATIC int32 LSAME(int32 CA,
00726 int32 CB)
00727 {
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737 int32 LSAME_v;
00738 int32 INTA,
00739 INTB,
00740 ZCODE;
00741
00742 DEBUG_ENTRY( "LSAME()" );
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768 LSAME_v = CA == CB;
00769 if( LSAME_v )
00770 {
00771 return LSAME_v;
00772 }
00773
00774
00775
00776 ZCODE = 'Z';
00777
00778
00779
00780
00781
00782
00783 INTA = (CA);
00784 INTB = (CB);
00785
00786 if( ZCODE == 90 || ZCODE == 122 )
00787 {
00788
00789
00790
00791
00792 if( INTA >= 97 && INTA <= 122 )
00793 INTA -= 32;
00794 if( INTB >= 97 && INTB <= 122 )
00795 INTB -= 32;
00796
00797 }
00798 else if( ZCODE == 233 || ZCODE == 169 )
00799 {
00800
00801
00802
00803
00804 if( ((INTA >= 129 && INTA <= 137) || (INTA >= 145 && INTA <=
00805 153)) || (INTA >= 162 && INTA <= 169) )
00806 INTA += 64;
00807 if( ((INTB >= 129 && INTB <= 137) || (INTB >= 145 && INTB <=
00808 153)) || (INTB >= 162 && INTB <= 169) )
00809 INTB += 64;
00810
00811 }
00812 else if( ZCODE == 218 || ZCODE == 250 )
00813 {
00814
00815
00816
00817
00818 if( INTA >= 225 && INTA <= 250 )
00819 INTA -= 32;
00820 if( INTB >= 225 && INTB <= 250 )
00821 INTB -= 32;
00822 }
00823 LSAME_v = INTA == INTB;
00824
00825
00826
00827
00828
00829 return LSAME_v;
00830 }
00831
00832
00833
00834 STATIC int32 IDAMAX(int32 n,
00835 double dx[],
00836 int32 incx)
00837 {
00838
00839
00840
00841
00842
00843
00844
00845 int32 IDAMAX_v,
00846 i,
00847 ix;
00848 double dmax;
00849
00850 DEBUG_ENTRY( "IDAMAX()" );
00851
00852 IDAMAX_v = 0;
00853
00854 if( n < 1 || incx <= 0 )
00855 {
00856 return IDAMAX_v;
00857 }
00858
00859 IDAMAX_v = 1;
00860
00861 if( n == 1 )
00862 {
00863 return IDAMAX_v;
00864 }
00865
00866 if( incx == 1 )
00867 goto L_20;
00868
00869
00870
00871 ix = 1;
00872 dmax = fabs(dx[0]);
00873 ix += incx;
00874 for( i=2; i <= n; i++ )
00875 {
00876
00877 if( fabs(dx[ix-1]) > dmax )
00878 {
00879 IDAMAX_v = i;
00880 dmax = fabs(dx[ix-1]);
00881 }
00882 ix += incx;
00883 }
00884 return IDAMAX_v;
00885
00886
00887
00888 L_20:
00889 dmax = fabs(dx[0]);
00890 for( i=1; i < n; i++ )
00891 {
00892
00893 if( fabs(dx[i]) > dmax )
00894 {
00895 IDAMAX_v = i+1;
00896 dmax = fabs(dx[i]);
00897 }
00898 }
00899 return IDAMAX_v;
00900 }
00901
00902
00903 STATIC void DTRSM(int32 SIDE,
00904 int32 UPLO,
00905 int32 TRANSA,
00906 int32 DIAG,
00907 int32 M,
00908 int32 N,
00909 double ALPHA,
00910 double *A,
00911 int32 LDA,
00912 double *B,
00913 int32 LDB)
00914 {
00915 int32 LSIDE,
00916 NOUNIT,
00917 UPPER;
00918 int32 I,
00919 INFO,
00920 I_,
00921 J,
00922 J_,
00923 K,
00924 K_,
00925 NROWA;
00926 double TEMP;
00927
00928 DEBUG_ENTRY( "DTRSM()" );
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
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 LSIDE = LSAME(SIDE,'L');
01063 if( LSIDE )
01064 {
01065 NROWA = M;
01066 }
01067 else
01068 {
01069 NROWA = N;
01070 }
01071 NOUNIT = LSAME(DIAG,'N');
01072 UPPER = LSAME(UPLO,'U');
01073
01074 INFO = 0;
01075 if( (!LSIDE) && (!LSAME(SIDE,'R')) )
01076 {
01077 INFO = 1;
01078 }
01079 else if( (!UPPER) && (!LSAME(UPLO,'L')) )
01080 {
01081 INFO = 2;
01082 }
01083 else if( ((!LSAME(TRANSA,'N')) && (!LSAME(TRANSA,'T'))) && (!LSAME(TRANSA,
01084 'C')) )
01085 {
01086 INFO = 3;
01087 }
01088 else if( (!LSAME(DIAG,'U')) && (!LSAME(DIAG,'N')) )
01089 {
01090 INFO = 4;
01091 }
01092 else if( M < 0 )
01093 {
01094 INFO = 5;
01095 }
01096 else if( N < 0 )
01097 {
01098 INFO = 6;
01099 }
01100 else if( LDA < MAX2(1,NROWA) )
01101 {
01102 INFO = 9;
01103 }
01104 else if( LDB < MAX2(1,M) )
01105 {
01106 INFO = 11;
01107 }
01108 if( INFO != 0 )
01109 {
01110 XERBLA("DTRSM ",INFO);
01111
01112 }
01113
01114
01115
01116 if( N == 0 )
01117 {
01118 return;}
01119
01120
01121
01122 if( ALPHA == ZERO )
01123 {
01124 for( J=1; J <= N; J++ )
01125 {
01126 J_ = J - 1;
01127 for( I=1; I <= M; I++ )
01128 {
01129 I_ = I - 1;
01130 BB(J_,I_) = ZERO;
01131 }
01132 }
01133 return;
01134 }
01135
01136
01137
01138 if( LSIDE )
01139 {
01140 if( LSAME(TRANSA,'N') )
01141 {
01142
01143
01144
01145 if( UPPER )
01146 {
01147 for( J=1; J <= N; J++ )
01148 {
01149 J_ = J - 1;
01150 if( ALPHA != ONE )
01151 {
01152 for( I=1; I <= M; I++ )
01153 {
01154 I_ = I - 1;
01155 BB(J_,I_) *= ALPHA;
01156 }
01157 }
01158 for( K=M; K >= 1; K-- )
01159 {
01160 K_ = K - 1;
01161 if( BB(J_,K_) != ZERO )
01162 {
01163 if( NOUNIT )
01164 BB(J_,K_) /= AA(K_,K_);
01165 for( I=1; I <= (K - 1); I++ )
01166 {
01167 I_ = I - 1;
01168 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
01169 }
01170 }
01171 }
01172 }
01173 }
01174 else
01175 {
01176 for( J=1; J <= N; J++ )
01177 {
01178 J_ = J - 1;
01179 if( ALPHA != ONE )
01180 {
01181 for( I=1; I <= M; I++ )
01182 {
01183 I_ = I - 1;
01184 BB(J_,I_) *= ALPHA;
01185 }
01186 }
01187 for( K=1; K <= M; K++ )
01188 {
01189 K_ = K - 1;
01190 if( BB(J_,K_) != ZERO )
01191 {
01192 if( NOUNIT )
01193 BB(J_,K_) /= AA(K_,K_);
01194 for( I=K + 1; I <= M; I++ )
01195 {
01196 I_ = I - 1;
01197 BB(J_,I_) += -BB(J_,K_)*AA(K_,I_);
01198 }
01199 }
01200 }
01201 }
01202 }
01203 }
01204 else
01205 {
01206
01207
01208
01209 if( UPPER )
01210 {
01211 for( J=1; J <= N; J++ )
01212 {
01213 J_ = J - 1;
01214 for( I=1; I <= M; I++ )
01215 {
01216 I_ = I - 1;
01217 TEMP = ALPHA*BB(J_,I_);
01218 for( K=1; K <= (I - 1); K++ )
01219 {
01220 K_ = K - 1;
01221 TEMP += -AA(I_,K_)*BB(J_,K_);
01222 }
01223 if( NOUNIT )
01224 TEMP /= AA(I_,I_);
01225 BB(J_,I_) = TEMP;
01226 }
01227 }
01228 }
01229 else
01230 {
01231 for( J=1; J <= N; J++ )
01232 {
01233 J_ = J - 1;
01234 for( I=M; I >= 1; I-- )
01235 {
01236 I_ = I - 1;
01237 TEMP = ALPHA*BB(J_,I_);
01238 for( K=I + 1; K <= M; K++ )
01239 {
01240 K_ = K - 1;
01241 TEMP += -AA(I_,K_)*BB(J_,K_);
01242 }
01243 if( NOUNIT )
01244 TEMP /= AA(I_,I_);
01245 BB(J_,I_) = TEMP;
01246 }
01247 }
01248 }
01249 }
01250 }
01251 else
01252 {
01253 if( LSAME(TRANSA,'N') )
01254 {
01255
01256
01257
01258 if( UPPER )
01259 {
01260 for( J=1; J <= N; J++ )
01261 {
01262 J_ = J - 1;
01263 if( ALPHA != ONE )
01264 {
01265 for( I=1; I <= M; I++ )
01266 {
01267 I_ = I - 1;
01268 BB(J_,I_) *= ALPHA;
01269 }
01270 }
01271 for( K=1; K <= (J - 1); K++ )
01272 {
01273 K_ = K - 1;
01274 if( AA(J_,K_) != ZERO )
01275 {
01276 for( I=1; I <= M; I++ )
01277 {
01278 I_ = I - 1;
01279 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
01280 }
01281 }
01282 }
01283 if( NOUNIT )
01284 {
01285 TEMP = ONE/AA(J_,J_);
01286 for( I=1; I <= M; I++ )
01287 {
01288 I_ = I - 1;
01289 BB(J_,I_) *= TEMP;
01290 }
01291 }
01292 }
01293 }
01294 else
01295 {
01296 for( J=N; J >= 1; J-- )
01297 {
01298 J_ = J - 1;
01299 if( ALPHA != ONE )
01300 {
01301 for( I=1; I <= M; I++ )
01302 {
01303 I_ = I - 1;
01304 BB(J_,I_) *= ALPHA;
01305 }
01306 }
01307 for( K=J + 1; K <= N; K++ )
01308 {
01309 K_ = K - 1;
01310 if( AA(J_,K_) != ZERO )
01311 {
01312 for( I=1; I <= M; I++ )
01313 {
01314 I_ = I - 1;
01315 BB(J_,I_) += -AA(J_,K_)*BB(K_,I_);
01316 }
01317 }
01318 }
01319 if( NOUNIT )
01320 {
01321 TEMP = ONE/AA(J_,J_);
01322 for( I=1; I <= M; I++ )
01323 {
01324 I_ = I - 1;
01325 BB(J_,I_) *= TEMP;
01326 }
01327 }
01328 }
01329 }
01330 }
01331 else
01332 {
01333
01334
01335
01336 if( UPPER )
01337 {
01338 for( K=N; K >= 1; K-- )
01339 {
01340 K_ = K - 1;
01341 if( NOUNIT )
01342 {
01343 TEMP = ONE/AA(K_,K_);
01344 for( I=1; I <= M; I++ )
01345 {
01346 I_ = I - 1;
01347 BB(K_,I_) *= TEMP;
01348 }
01349 }
01350 for( J=1; J <= (K - 1); J++ )
01351 {
01352 J_ = J - 1;
01353 if( AA(K_,J_) != ZERO )
01354 {
01355 TEMP = AA(K_,J_);
01356 for( I=1; I <= M; I++ )
01357 {
01358 I_ = I - 1;
01359 BB(J_,I_) += -TEMP*BB(K_,I_);
01360 }
01361 }
01362 }
01363 if( ALPHA != ONE )
01364 {
01365 for( I=1; I <= M; I++ )
01366 {
01367 I_ = I - 1;
01368 BB(K_,I_) *= ALPHA;
01369 }
01370 }
01371 }
01372 }
01373 else
01374 {
01375 for( K=1; K <= N; K++ )
01376 {
01377 K_ = K - 1;
01378 if( NOUNIT )
01379 {
01380 TEMP = ONE/AA(K_,K_);
01381 for( I=1; I <= M; I++ )
01382 {
01383 I_ = I - 1;
01384 BB(K_,I_) *= TEMP;
01385 }
01386 }
01387 for( J=K + 1; J <= N; J++ )
01388 {
01389 J_ = J - 1;
01390 if( AA(K_,J_) != ZERO )
01391 {
01392 TEMP = AA(K_,J_);
01393 for( I=1; I <= M; I++ )
01394 {
01395 I_ = I - 1;
01396 BB(J_,I_) += -TEMP*BB(K_,I_);
01397 }
01398 }
01399 }
01400 if( ALPHA != ONE )
01401 {
01402 for( I=1; I <= M; I++ )
01403 {
01404 I_ = I - 1;
01405 BB(K_,I_) *= ALPHA;
01406 }
01407 }
01408 }
01409 }
01410 }
01411 }
01412
01413 return;
01414
01415
01416
01417 #undef B
01418 #undef A
01419 }
01420
01421
01422
01423
01424 STATIC int32 ILAENV(int32 ISPEC,
01425 const char *NAME,
01426
01427 int32 N1,
01428 int32 N2,
01429
01430 int32 N4)
01431 {
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 char C2[3],
01442 C3[4],
01443 C4[3],
01444 SUBNAM[7];
01445 int32 CNAME,
01446 SNAME;
01447 char C1;
01448 int32 I,
01449 IC,
01450 ILAENV_v,
01451 IZ,
01452 NB,
01453 NBMIN,
01454 NX;
01455
01456 DEBUG_ENTRY( "ILAENV()" );
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550 switch( ISPEC )
01551 {
01552 case 1:
01553 {
01554 goto L_100;
01555 }
01556 case 2:
01557 {
01558 goto L_100;
01559 }
01560 case 3:
01561 {
01562 goto L_100;
01563 }
01564 case 4:
01565 {
01566 goto L_400;
01567 }
01568 case 5:
01569 {
01570 goto L_500;
01571 }
01572 case 6:
01573 {
01574 goto L_600;
01575 }
01576 case 7:
01577 {
01578 goto L_700;
01579 }
01580 case 8:
01581 {
01582 goto L_800;
01583 }
01584
01585 default:
01586 {
01587
01588
01589 ILAENV_v = -1;
01590 return ILAENV_v;
01591 }
01592 }
01593
01594 L_100:
01595
01596
01597
01598 ILAENV_v = 1;
01599 strncpy( SUBNAM, NAME, 6 );
01600 SUBNAM[6] = '\0';
01601 IC = (SUBNAM[0]);
01602 IZ = 'Z';
01603 if( IZ == 90 || IZ == 122 )
01604 {
01605
01606
01607
01608 if( IC >= 97 && IC <= 122 )
01609 {
01610 SUBNAM[0] = (char)(IC - 32);
01611 for( I=2; I <= 6; I++ )
01612 {
01613 IC = (SUBNAM[I-1]);
01614 if( IC >= 97 && IC <= 122 )
01615 SUBNAM[I - 1] = (char)(IC - 32);
01616 }
01617 }
01618
01619 }
01620 else if( IZ == 233 || IZ == 169 )
01621 {
01622
01623
01624
01625 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) ||
01626 (IC >= 162 && IC <= 169) )
01627 {
01628 SUBNAM[0] = (char)(IC + 64);
01629 for( I=2; I <= 6; I++ )
01630 {
01631 IC = (SUBNAM[I-1]);
01632 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <=
01633 153)) || (IC >= 162 && IC <= 169) )
01634 SUBNAM[I - 1] = (char)(IC + 64);
01635 }
01636 }
01637
01638 }
01639 else if( IZ == 218 || IZ == 250 )
01640 {
01641
01642
01643
01644 if( IC >= 225 && IC <= 250 )
01645 {
01646 SUBNAM[0] = (char)(IC - 32);
01647 for( I=2; I <= 6; I++ )
01648 {
01649 IC = (SUBNAM[I-1]);
01650 if( IC >= 225 && IC <= 250 )
01651 SUBNAM[I - 1] = (char)(IC - 32);
01652 }
01653 }
01654 }
01655
01656 C1 = SUBNAM[0];
01657 SNAME = C1 == 'S' || C1 == 'D';
01658 CNAME = C1 == 'C' || C1 == 'Z';
01659 if( !(CNAME || SNAME) )
01660 {
01661 return ILAENV_v;
01662 }
01663
01664 # if 0
01665 strncpy( C2, SUBNAM+1, 2 );
01666 strncpy( C3, SUBNAM+3, 3 );
01667 strncpy( C4, C3+1, 2 );
01668 # endif
01669
01670
01671
01672 strncpy( C2, SUBNAM+1, 2 );
01673 C2[2] = '\0';
01674 strncpy( C3, SUBNAM+3, 3 );
01675 C3[3] = '\0';
01676 strncpy( C4, C3+1, 2 );
01677 C4[2] = '\0';
01678
01679 switch( ISPEC )
01680 {
01681 case 1: goto L_110;
01682 case 2: goto L_200;
01683 case 3: goto L_300;
01684 }
01685
01686 L_110:
01687
01688
01689
01690
01691
01692
01693
01694 NB = 1;
01695
01696 if( strcmp(C2,"GE") == 0 )
01697 {
01698 if( strcmp(C3,"TRF") == 0 )
01699 {
01700 if( SNAME )
01701 {
01702 NB = 64;
01703 }
01704 else
01705 {
01706 NB = 64;
01707 }
01708 }
01709 else if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) ||
01710 strcmp(C3,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01711 {
01712 if( SNAME )
01713 {
01714 NB = 32;
01715 }
01716 else
01717 {
01718 NB = 32;
01719 }
01720 }
01721 else if( strcmp(C3,"HRD") == 0 )
01722 {
01723 if( SNAME )
01724 {
01725 NB = 32;
01726 }
01727 else
01728 {
01729 NB = 32;
01730 }
01731 }
01732 else if( strcmp(C3,"BRD") == 0 )
01733 {
01734 if( SNAME )
01735 {
01736 NB = 32;
01737 }
01738 else
01739 {
01740 NB = 32;
01741 }
01742 }
01743 else if( strcmp(C3,"TRI") == 0 )
01744 {
01745 if( SNAME )
01746 {
01747 NB = 64;
01748 }
01749 else
01750 {
01751 NB = 64;
01752 }
01753 }
01754 }
01755 else if( strcmp(C2,"PO") == 0 )
01756 {
01757 if( strcmp(C3,"TRF") == 0 )
01758 {
01759 if( SNAME )
01760 {
01761 NB = 64;
01762 }
01763 else
01764 {
01765 NB = 64;
01766 }
01767 }
01768 }
01769 else if( strcmp(C2,"SY") == 0 )
01770 {
01771 if( strcmp(C3,"TRF") == 0 )
01772 {
01773 if( SNAME )
01774 {
01775 NB = 64;
01776 }
01777 else
01778 {
01779 NB = 64;
01780 }
01781 }
01782 else if( SNAME && strcmp(C3,"TRD") == 0 )
01783 {
01784 NB = 1;
01785 }
01786 else if( SNAME && strcmp(C3,"GST") == 0 )
01787 {
01788 NB = 64;
01789 }
01790 }
01791 else if( CNAME && strcmp(C2,"HE") == 0 )
01792 {
01793 if( strcmp(C3,"TRF") == 0 )
01794 {
01795 NB = 64;
01796 }
01797 else if( strcmp(C3,"TRD") == 0 )
01798 {
01799 NB = 1;
01800 }
01801 else if( strcmp(C3,"GST") == 0 )
01802 {
01803 NB = 64;
01804 }
01805 }
01806 else if( SNAME && strcmp(C2,"OR") == 0 )
01807 {
01808 if( C3[0] == 'G' )
01809 {
01810 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01811 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01812 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01813 0 )
01814 {
01815 NB = 32;
01816 }
01817 }
01818 else if( C3[0] == 'M' )
01819 {
01820 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01821 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01822 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01823 0 )
01824 {
01825 NB = 32;
01826 }
01827 }
01828 }
01829 else if( CNAME && strcmp(C2,"UN") == 0 )
01830 {
01831 if( C3[0] == 'G' )
01832 {
01833 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01834 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01835 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01836 0 )
01837 {
01838 NB = 32;
01839 }
01840 }
01841 else if( C3[0] == 'M' )
01842 {
01843 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01844 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01845 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01846 0 )
01847 {
01848 NB = 32;
01849 }
01850 }
01851 }
01852 else if( strcmp(C2,"GB") == 0 )
01853 {
01854 if( strcmp(C3,"TRF") == 0 )
01855 {
01856 if( SNAME )
01857 {
01858 if( N4 <= 64 )
01859 {
01860 NB = 1;
01861 }
01862 else
01863 {
01864 NB = 32;
01865 }
01866 }
01867 else
01868 {
01869 if( N4 <= 64 )
01870 {
01871 NB = 1;
01872 }
01873 else
01874 {
01875 NB = 32;
01876 }
01877 }
01878 }
01879 }
01880 else if( strcmp(C2,"PB") == 0 )
01881 {
01882 if( strcmp(C3,"TRF") == 0 )
01883 {
01884 if( SNAME )
01885 {
01886 if( N2 <= 64 )
01887 {
01888 NB = 1;
01889 }
01890 else
01891 {
01892 NB = 32;
01893 }
01894 }
01895 else
01896 {
01897 if( N2 <= 64 )
01898 {
01899 NB = 1;
01900 }
01901 else
01902 {
01903 NB = 32;
01904 }
01905 }
01906 }
01907 }
01908 else if( strcmp(C2,"TR") == 0 )
01909 {
01910 if( strcmp(C3,"TRI") == 0 )
01911 {
01912 if( SNAME )
01913 {
01914 NB = 64;
01915 }
01916 else
01917 {
01918 NB = 64;
01919 }
01920 }
01921 }
01922 else if( strcmp(C2,"LA") == 0 )
01923 {
01924 if( strcmp(C3,"UUM") == 0 )
01925 {
01926 if( SNAME )
01927 {
01928 NB = 64;
01929 }
01930 else
01931 {
01932 NB = 64;
01933 }
01934 }
01935 }
01936 else if( SNAME && strcmp(C2,"ST") == 0 )
01937 {
01938 if( strcmp(C3,"EBZ") == 0 )
01939 {
01940 NB = 1;
01941 }
01942 }
01943 ILAENV_v = NB;
01944 return ILAENV_v;
01945
01946 L_200:
01947
01948
01949
01950 NBMIN = 2;
01951 if( strcmp(C2,"GE") == 0 )
01952 {
01953 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
01954 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01955 {
01956 if( SNAME )
01957 {
01958 NBMIN = 2;
01959 }
01960 else
01961 {
01962 NBMIN = 2;
01963 }
01964 }
01965 else if( strcmp(C3,"HRD") == 0 )
01966 {
01967 if( SNAME )
01968 {
01969 NBMIN = 2;
01970 }
01971 else
01972 {
01973 NBMIN = 2;
01974 }
01975 }
01976 else if( strcmp(C3,"BRD") == 0 )
01977 {
01978 if( SNAME )
01979 {
01980 NBMIN = 2;
01981 }
01982 else
01983 {
01984 NBMIN = 2;
01985 }
01986 }
01987 else if( strcmp(C3,"TRI") == 0 )
01988 {
01989 if( SNAME )
01990 {
01991 NBMIN = 2;
01992 }
01993 else
01994 {
01995 NBMIN = 2;
01996 }
01997 }
01998 }
01999 else if( strcmp(C2,"SY") == 0 )
02000 {
02001 if( strcmp(C3,"TRF") == 0 )
02002 {
02003 if( SNAME )
02004 {
02005 NBMIN = 8;
02006 }
02007 else
02008 {
02009 NBMIN = 8;
02010 }
02011 }
02012 else if( SNAME && strcmp(C3,"TRD") == 0 )
02013 {
02014 NBMIN = 2;
02015 }
02016 }
02017 else if( CNAME && strcmp(C2,"HE") == 0 )
02018 {
02019 if( strcmp(C3,"TRD") == 0 )
02020 {
02021 NBMIN = 2;
02022 }
02023 }
02024 else if( SNAME && strcmp(C2,"OR") == 0 )
02025 {
02026 if( C3[0] == 'G' )
02027 {
02028 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02029 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02030 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02031 0 )
02032 {
02033 NBMIN = 2;
02034 }
02035 }
02036 else if( C3[0] == 'M' )
02037 {
02038 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02039 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02040 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02041 0 )
02042 {
02043 NBMIN = 2;
02044 }
02045 }
02046 }
02047 else if( CNAME && strcmp(C2,"UN") == 0 )
02048 {
02049 if( C3[0] == 'G' )
02050 {
02051 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02052 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02053 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02054 0 )
02055 {
02056 NBMIN = 2;
02057 }
02058 }
02059 else if( C3[0] == 'M' )
02060 {
02061 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02062 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02063 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02064 0 )
02065 {
02066 NBMIN = 2;
02067 }
02068 }
02069 }
02070 ILAENV_v = NBMIN;
02071 return ILAENV_v;
02072
02073 L_300:
02074
02075
02076
02077 NX = 0;
02078 if( strcmp(C2,"GE") == 0 )
02079 {
02080 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
02081 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
02082 {
02083 if( SNAME )
02084 {
02085 NX = 128;
02086 }
02087 else
02088 {
02089 NX = 128;
02090 }
02091 }
02092 else if( strcmp(C3,"HRD") == 0 )
02093 {
02094 if( SNAME )
02095 {
02096 NX = 128;
02097 }
02098 else
02099 {
02100 NX = 128;
02101 }
02102 }
02103 else if( strcmp(C3,"BRD") == 0 )
02104 {
02105 if( SNAME )
02106 {
02107 NX = 128;
02108 }
02109 else
02110 {
02111 NX = 128;
02112 }
02113 }
02114 }
02115 else if( strcmp(C2,"SY") == 0 )
02116 {
02117 if( SNAME && strcmp(C3,"TRD") == 0 )
02118 {
02119 NX = 1;
02120 }
02121 }
02122 else if( CNAME && strcmp(C2,"HE") == 0 )
02123 {
02124 if( strcmp(C3,"TRD") == 0 )
02125 {
02126 NX = 1;
02127 }
02128 }
02129 else if( SNAME && strcmp(C2,"OR") == 0 )
02130 {
02131 if( C3[0] == 'G' )
02132 {
02133 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02134 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02135 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02136 0 )
02137 {
02138 NX = 128;
02139 }
02140 }
02141 }
02142 else if( CNAME && strcmp(C2,"UN") == 0 )
02143 {
02144 if( C3[0] == 'G' )
02145 {
02146 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02147 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02148 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02149 0 )
02150 {
02151 NX = 128;
02152 }
02153 }
02154 }
02155 ILAENV_v = NX;
02156 return ILAENV_v;
02157
02158 L_400:
02159
02160
02161
02162 ILAENV_v = 6;
02163 return ILAENV_v;
02164
02165 L_500:
02166
02167
02168
02169 ILAENV_v = 2;
02170 return ILAENV_v;
02171
02172 L_600:
02173
02174
02175
02176 ILAENV_v = (int32)((realnum)(MIN2(N1,N2))*1.6e0);
02177 return ILAENV_v;
02178
02179 L_700:
02180
02181
02182
02183 ILAENV_v = 1;
02184 return ILAENV_v;
02185
02186 L_800:
02187
02188
02189
02190 ILAENV_v = 50;
02191 return ILAENV_v;
02192
02193
02194
02195 }
02196
02197
02198
02199 STATIC void DSWAP(int32 n,
02200 double dx[],
02201 int32 incx,
02202 double dy[],
02203 int32 incy)
02204 {
02205 int32 i,
02206 ix,
02207 iy,
02208 m;
02209 double dtemp;
02210
02211 DEBUG_ENTRY( "DSWAP()" );
02212
02213
02214
02215
02216
02217
02218
02219 if( n <= 0 )
02220 {
02221 return;}
02222 if( incx == 1 && incy == 1 )
02223 goto L_20;
02224
02225
02226
02227
02228 ix = 1;
02229 iy = 1;
02230
02231 if( incx < 0 )
02232 ix = (-n + 1)*incx + 1;
02233
02234 if( incy < 0 )
02235 iy = (-n + 1)*incy + 1;
02236
02237 for( i=0; i < n; i++ )
02238 {
02239 dtemp = dx[ix-1];
02240 dx[ix-1] = dy[iy-1];
02241 dy[iy-1] = dtemp;
02242 ix += incx;
02243 iy += incy;
02244 }
02245 return;
02246
02247
02248
02249
02250
02251
02252 L_20:
02253 m = n%3;
02254 if( m == 0 )
02255 goto L_40;
02256
02257 for( i=0; i < m; i++ )
02258 {
02259 dtemp = dx[i];
02260 dx[i] = dy[i];
02261 dy[i] = dtemp;
02262 }
02263
02264 if( n < 3 )
02265 {
02266 return;
02267 }
02268
02269 L_40:
02270 for( i=m; i < n; i += 3 )
02271 {
02272 dtemp = dx[i];
02273 dx[i] = dy[i];
02274 dy[i] = dtemp;
02275 dtemp = dx[i+1];
02276 dx[i+1] = dy[i+1];
02277 dy[i+1] = dtemp;
02278 dtemp = dx[i+2];
02279 dx[i+2] = dy[i+2];
02280 dy[i+2] = dtemp;
02281 }
02282 return;
02283 }
02284
02285
02286
02287 STATIC void DSCAL(int32 n,
02288 double da,
02289 double dx[],
02290 int32 incx)
02291 {
02292 int32 i,
02293 m,
02294 nincx;
02295
02296 DEBUG_ENTRY( "DSCAL()" );
02297
02298
02299
02300
02301
02302
02303
02304
02305 if( n <= 0 || incx <= 0 )
02306 {
02307 return;}
02308 if( incx == 1 )
02309 goto L_20;
02310
02311
02312
02313 nincx = n*incx;
02314
02315 for( i=0; i<nincx; i = i + incx)
02316 {
02317 dx[i] *= da;
02318 }
02319 return;
02320
02321
02322
02323
02324
02325
02326 L_20:
02327 m = n%5;
02328 if( m == 0 )
02329 goto L_40;
02330
02331 for( i=0; i < m; i++ )
02332 {
02333 dx[i] *= da;
02334 }
02335
02336 if( n < 5 )
02337 {
02338 return;
02339 }
02340
02341 L_40:
02342
02343 for( i=m; i < n; i += 5 )
02344 {
02345 dx[i] *= da;
02346 dx[i+1] *= da;
02347 dx[i+2] *= da;
02348 dx[i+3] *= da;
02349 dx[i+4] *= da;
02350 }
02351 return;
02352 }
02353
02354
02355
02356 STATIC void DLASWP(int32 N,
02357 double *A,
02358 int32 LDA,
02359 int32 K1,
02360 int32 K2,
02361 int32 IPIV[],
02362 int32 INCX)
02363 {
02364 int32 I,
02365 IP,
02366 IX,
02367 I_;
02368
02369 DEBUG_ENTRY( "DLASWP()" );
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428 if( INCX == 0 )
02429 {
02430 return;}
02431 if( INCX > 0 )
02432 {
02433 IX = K1;
02434 }
02435 else
02436 {
02437 IX = 1 + (1 - K2)*INCX;
02438 }
02439 if( INCX == 1 )
02440 {
02441 for( I=K1; I <= K2; I++ )
02442 {
02443 I_ = I - 1;
02444 IP = IPIV[I_];
02445 if( IP != I )
02446 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02447 }
02448 }
02449 else if( INCX > 1 )
02450 {
02451 for( I=K1; I <= K2; I++ )
02452 {
02453 I_ = I - 1;
02454 IP = IPIV[IX-1];
02455 if( IP != I )
02456 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02457 IX += INCX;
02458 }
02459 }
02460 else if( INCX < 0 )
02461 {
02462 for( I=K2; I >= K1; I-- )
02463 {
02464 I_ = I - 1;
02465 IP = IPIV[IX-1];
02466 if( IP != I )
02467 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02468 IX += INCX;
02469 }
02470 }
02471
02472 return;
02473
02474
02475
02476 #undef A
02477 }
02478
02479
02480
02481 STATIC void DGETF2(int32 M,
02482 int32 N,
02483 double *A,
02484 int32 LDA,
02485 int32 IPIV[],
02486 int32 *INFO)
02487 {
02488 int32 J,
02489 JP,
02490 J_,
02491 limit;
02492
02493 DEBUG_ENTRY( "DGETF2()" );
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564 *INFO = 0;
02565 if( M < 0 )
02566 {
02567 *INFO = -1;
02568 }
02569 else if( N < 0 )
02570 {
02571 *INFO = -2;
02572 }
02573 else if( LDA < MAX2(1,M) )
02574 {
02575 *INFO = -4;
02576 }
02577 if( *INFO != 0 )
02578 {
02579 XERBLA("DGETF2",-*INFO);
02580
02581 }
02582
02583
02584
02585 if( M == 0 || N == 0 )
02586 {
02587 return;}
02588
02589 limit = MIN2(M,N);
02590 for( J=1; J <= limit; J++ )
02591 {
02592 J_ = J - 1;
02593
02594
02595
02596 JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
02597 IPIV[J_] = JP;
02598 if( AA(J_,JP-1) != ZERO )
02599 {
02600
02601
02602 if( JP != J )
02603 DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
02604
02605
02606
02607 if( J < M )
02608 DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
02609 }
02610 else if( *INFO == 0 )
02611 {
02612 *INFO = J;
02613 }
02614
02615 if( J < MIN2(M,N) )
02616 {
02617
02618
02619 DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
02620 LDA);
02621 }
02622 }
02623 return;
02624
02625
02626
02627 #undef A
02628 }
02629
02630
02631
02632 STATIC void DGER(int32 M,
02633 int32 N,
02634 double ALPHA,
02635 double X[],
02636 int32 INCX,
02637 double Y[],
02638 int32 INCY,
02639 double *A,
02640 int32 LDA)
02641 {
02642 int32 I,
02643 INFO,
02644 IX,
02645 I_,
02646 J,
02647 JY,
02648 J_,
02649 KX;
02650 double TEMP;
02651
02652 DEBUG_ENTRY( "DGER()" );
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734
02735
02736 INFO = 0;
02737 if( M < 0 )
02738 {
02739 INFO = 1;
02740 }
02741 else if( N < 0 )
02742 {
02743 INFO = 2;
02744 }
02745 else if( INCX == 0 )
02746 {
02747 INFO = 5;
02748 }
02749 else if( INCY == 0 )
02750 {
02751 INFO = 7;
02752 }
02753 else if( LDA < MAX2(1,M) )
02754 {
02755 INFO = 9;
02756 }
02757 if( INFO != 0 )
02758 {
02759 XERBLA("DGER ",INFO);
02760
02761 }
02762
02763
02764
02765 if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
02766 {
02767 return;}
02768
02769
02770
02771
02772 if( INCY > 0 )
02773 {
02774 JY = 1;
02775 }
02776 else
02777 {
02778 JY = 1 - (N - 1)*INCY;
02779 }
02780 if( INCX == 1 )
02781 {
02782 for( J=1; J <= N; J++ )
02783 {
02784 J_ = J - 1;
02785 if( Y[JY-1] != ZERO )
02786 {
02787 TEMP = ALPHA*Y[JY-1];
02788 for( I=1; I <= M; I++ )
02789 {
02790 I_ = I - 1;
02791 AA(J_,I_) += X[I_]*TEMP;
02792 }
02793 }
02794 JY += INCY;
02795 }
02796 }
02797 else
02798 {
02799 if( INCX > 0 )
02800 {
02801 KX = 1;
02802 }
02803 else
02804 {
02805 KX = 1 - (M - 1)*INCX;
02806 }
02807 for( J=1; J <= N; J++ )
02808 {
02809 J_ = J - 1;
02810 if( Y[JY-1] != ZERO )
02811 {
02812 TEMP = ALPHA*Y[JY-1];
02813 IX = KX;
02814 for( I=1; I <= M; I++ )
02815 {
02816 I_ = I - 1;
02817 AA(J_,I_) += X[IX-1]*TEMP;
02818 IX += INCX;
02819 }
02820 }
02821 JY += INCY;
02822 }
02823 }
02824
02825 return;
02826
02827
02828
02829 #undef A
02830 }
02831
02832
02833
02834 STATIC void DGEMM(int32 TRANSA,
02835 int32 TRANSB,
02836 int32 M,
02837 int32 N,
02838 int32 K,
02839 double ALPHA,
02840 double *A,
02841 int32 LDA,
02842 double *B,
02843 int32 LDB,
02844 double BETA,
02845 double *C,
02846 int32 LDC)
02847 {
02848 int32 NOTA,
02849 NOTB;
02850 int32 I,
02851 INFO,
02852 I_,
02853 J,
02854 J_,
02855 L,
02856 L_,
02857 NROWA,
02858 NROWB;
02859 double TEMP;
02860
02861 DEBUG_ENTRY( "DGEMM()" );
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926
02927
02928
02929
02930
02931
02932
02933
02934
02935
02936
02937
02938
02939
02940
02941
02942
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990
02991
02992
02993
02994
02995
02996
02997
02998 NOTA = LSAME(TRANSA,'N');
02999 NOTB = LSAME(TRANSB,'N');
03000
03001 if( NOTA )
03002 {
03003 NROWA = M;
03004 }
03005 else
03006 {
03007 NROWA = K;
03008 }
03009
03010 if( NOTB )
03011 {
03012 NROWB = K;
03013 }
03014 else
03015 {
03016 NROWB = N;
03017 }
03018
03019
03020
03021 INFO = 0;
03022 if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
03023 {
03024 INFO = 1;
03025 }
03026 else if(
03027 ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
03028 {
03029 INFO = 2;
03030 }
03031
03032 else if( M < 0 )
03033 {
03034 INFO = 3;
03035 }
03036
03037 else if( N < 0 )
03038 {
03039 INFO = 4;
03040 }
03041
03042 else if( K < 0 )
03043 {
03044 INFO = 5;
03045 }
03046
03047 else if( LDA < MAX2(1,NROWA) )
03048 {
03049 INFO = 8;
03050 }
03051
03052 else if( LDB < MAX2(1,NROWB) )
03053 {
03054 INFO = 10;
03055 }
03056
03057 else if( LDC < MAX2(1,M) )
03058 {
03059 INFO = 13;
03060 }
03061
03062 if( INFO != 0 )
03063 {
03064 XERBLA("DGEMM ",INFO);
03065
03066 }
03067
03068
03069
03070 if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) &&
03071 (BETA == ONE)) )
03072 {
03073 return;
03074 }
03075
03076
03077 if( ALPHA == ZERO )
03078 {
03079 if( BETA == ZERO )
03080 {
03081 for( J=1; J <= N; J++ )
03082 {
03083 J_ = J - 1;
03084 for( I=1; I <= M; I++ )
03085 {
03086 I_ = I - 1;
03087 CC(J_,I_) = ZERO;
03088 }
03089 }
03090 }
03091
03092 else
03093 {
03094 for( J=1; J <= N; J++ )
03095 {
03096 J_ = J - 1;
03097 for( I=1; I <= M; I++ )
03098 {
03099 I_ = I - 1;
03100 CC(J_,I_) *= BETA;
03101 }
03102 }
03103 }
03104 return;
03105 }
03106
03107
03108 if( NOTB )
03109 {
03110
03111 if( NOTA )
03112 {
03113
03114
03115
03116 for( J=1; J <= N; J++ )
03117 {
03118 J_ = J - 1;
03119 if( BETA == ZERO )
03120 {
03121 for( I=1; I <= M; I++ )
03122 {
03123 I_ = I - 1;
03124 CC(J_,I_) = ZERO;
03125 }
03126 }
03127
03128 else if( BETA != ONE )
03129 {
03130 for( I=1; I <= M; I++ )
03131 {
03132 I_ = I - 1;
03133 CC(J_,I_) *= BETA;
03134 }
03135 }
03136
03137 for( L=1; L <= K; L++ )
03138 {
03139 L_ = L - 1;
03140 if( BB(J_,L_) != ZERO )
03141 {
03142 TEMP = ALPHA*BB(J_,L_);
03143 for( I=1; I <= M; I++ )
03144 {
03145 I_ = I - 1;
03146 CC(J_,I_) += TEMP*AA(L_,I_);
03147 }
03148 }
03149 }
03150 }
03151 }
03152 else
03153 {
03154
03155
03156 for( J=1; J <= N; J++ )
03157 {
03158 J_ = J - 1;
03159 for( I=1; I <= M; I++ )
03160 {
03161 I_ = I - 1;
03162 TEMP = ZERO;
03163 for( L=1; L <= K; L++ )
03164 {
03165 L_ = L - 1;
03166 TEMP += AA(I_,L_)*BB(J_,L_);
03167 }
03168
03169 if( BETA == ZERO )
03170 {
03171 CC(J_,I_) = ALPHA*TEMP;
03172 }
03173 else
03174 {
03175 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03176 }
03177 }
03178 }
03179 }
03180 }
03181 else
03182 {
03183 if( NOTA )
03184 {
03185
03186
03187
03188 for( J=1; J <= N; J++ )
03189 {
03190 J_ = J - 1;
03191 if( BETA == ZERO )
03192 {
03193 for( I=1; I <= M; I++ )
03194 {
03195 I_ = I - 1;
03196 CC(J_,I_) = ZERO;
03197 }
03198 }
03199
03200 else if( BETA != ONE )
03201 {
03202 for( I=1; I <= M; I++ )
03203 {
03204 I_ = I - 1;
03205 CC(J_,I_) *= BETA;
03206 }
03207 }
03208
03209 for( L=1; L <= K; L++ )
03210 {
03211 L_ = L - 1;
03212 if( BB(L_,J_) != ZERO )
03213 {
03214 TEMP = ALPHA*BB(L_,J_);
03215 for( I=1; I <= M; I++ )
03216 {
03217 I_ = I - 1;
03218 CC(J_,I_) += TEMP*AA(L_,I_);
03219 }
03220 }
03221 }
03222 }
03223 }
03224
03225 else
03226 {
03227
03228
03229 for( J=1; J <= N; J++ )
03230 {
03231 J_ = J - 1;
03232
03233 for( I=1; I <= M; I++ )
03234 {
03235 I_ = I - 1;
03236 TEMP = ZERO;
03237
03238 for( L=1; L <= K; L++ )
03239 {
03240 L_ = L - 1;
03241 TEMP += AA(I_,L_)*BB(L_,J_);
03242 }
03243
03244 if( BETA == ZERO )
03245 {
03246 CC(J_,I_) = ALPHA*TEMP;
03247 }
03248 else
03249 {
03250 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03251 }
03252
03253 }
03254 }
03255 }
03256 }
03257
03258 return;
03259
03260
03261 #undef C
03262 #undef B
03263 #undef A
03264 }
03265 #undef AA
03266 #undef BB
03267 #undef CC
03268
03269
03270 #if 0
03271 STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
03272 double *d__, double *du, double *b, int32 *ldb, int32
03273 *info)
03274 {
03275
03276 int32 b_dim1, b_offset, i__1, i__2;
03277
03278
03279 double fact, temp;
03280 int32 i__, j;
03281
03282
03283 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311
03312
03313
03314
03315
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327
03328
03329
03330
03331
03332
03333
03334
03335
03336
03337
03338
03339
03340
03341
03342
03343
03344
03345
03346
03347
03348
03349
03350
03351
03352
03353 --dl;
03354 --d__;
03355 --du;
03356 b_dim1 = *ldb;
03357 b_offset = 1 + b_dim1 * 1;
03358 b -= b_offset;
03359
03360
03361 *info = 0;
03362 if(*n < 0) {
03363 *info = -1;
03364 } else if(*nrhs < 0) {
03365 *info = -2;
03366 } else if(*ldb < *n && *ldb < 1) {
03367 *info = -7;
03368 }
03369 if(*info != 0) {
03370 i__1 = -(*info);
03371 XERBLA("DGTSV ", i__1);
03372
03373 }
03374
03375 if(*n == 0) {
03376 return 0;
03377 }
03378
03379 if(*nrhs == 1) {
03380 i__1 = *n - 2;
03381 for(i__ = 1; i__ <= i__1; ++i__) {
03382 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03383
03384
03385
03386 if(d__[i__] != 0.) {
03387 fact = dl[i__] / d__[i__];
03388 d__[i__ + 1] -= fact * du[i__];
03389 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
03390 1);
03391 } else {
03392 *info = i__;
03393 return 0;
03394 }
03395 dl[i__] = 0.;
03396 } else {
03397
03398
03399
03400 fact = d__[i__] / dl[i__];
03401 d__[i__] = dl[i__];
03402 temp = d__[i__ + 1];
03403 d__[i__ + 1] = du[i__] - fact * temp;
03404 dl[i__] = du[i__ + 1];
03405 du[i__ + 1] = -fact * dl[i__];
03406 du[i__] = temp;
03407 temp = b_ref(i__, 1);
03408 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03409 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03410 }
03411
03412 }
03413 if(*n > 1) {
03414 i__ = *n - 1;
03415 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03416 if(d__[i__] != 0.) {
03417 fact = dl[i__] / d__[i__];
03418 d__[i__ + 1] -= fact * du[i__];
03419 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
03420 1);
03421 } else {
03422 *info = i__;
03423 return 0;
03424 }
03425 } else {
03426 fact = d__[i__] / dl[i__];
03427 d__[i__] = dl[i__];
03428 temp = d__[i__ + 1];
03429 d__[i__ + 1] = du[i__] - fact * temp;
03430 du[i__] = temp;
03431 temp = b_ref(i__, 1);
03432 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03433 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03434 }
03435 }
03436 if(d__[*n] == 0.) {
03437 *info = *n;
03438 return 0;
03439 }
03440 } else {
03441 i__1 = *n - 2;
03442 for(i__ = 1; i__ <= i__1; ++i__) {
03443 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03444
03445
03446
03447 if(d__[i__] != 0.) {
03448 fact = dl[i__] / d__[i__];
03449 d__[i__ + 1] -= fact * du[i__];
03450 i__2 = *nrhs;
03451 for(j = 1; j <= i__2; ++j) {
03452 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03453 i__, j);
03454
03455 }
03456 } else {
03457 *info = i__;
03458 return 0;
03459 }
03460 dl[i__] = 0.;
03461 } else {
03462
03463
03464
03465 fact = d__[i__] / dl[i__];
03466 d__[i__] = dl[i__];
03467 temp = d__[i__ + 1];
03468 d__[i__ + 1] = du[i__] - fact * temp;
03469 dl[i__] = du[i__ + 1];
03470 du[i__ + 1] = -fact * dl[i__];
03471 du[i__] = temp;
03472 i__2 = *nrhs;
03473 for(j = 1; j <= i__2; ++j) {
03474 temp = b_ref(i__, j);
03475 b_ref(i__, j) = b_ref(i__ + 1, j);
03476 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03477
03478 }
03479 }
03480
03481 }
03482 if(*n > 1) {
03483 i__ = *n - 1;
03484 if( fabs(d__[i__]) >= fabs(dl[i__]))
03485 {
03486 if(d__[i__] != 0.)
03487 {
03488 fact = dl[i__] / d__[i__];
03489 d__[i__ + 1] -= fact * du[i__];
03490 i__1 = *nrhs;
03491 for(j = 1; j <= i__1; ++j) {
03492 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03493 i__, j);
03494
03495 }
03496 }
03497 else
03498 {
03499 *info = i__;
03500 return 0;
03501 }
03502 } else {
03503 fact = d__[i__] / dl[i__];
03504 d__[i__] = dl[i__];
03505 temp = d__[i__ + 1];
03506 d__[i__ + 1] = du[i__] - fact * temp;
03507 du[i__] = temp;
03508 i__1 = *nrhs;
03509 for(j = 1; j <= i__1; ++j) {
03510 temp = b_ref(i__, j);
03511 b_ref(i__, j) = b_ref(i__ + 1, j);
03512 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03513
03514 }
03515 }
03516 }
03517 if(d__[*n] == 0.) {
03518 *info = *n;
03519 return 0;
03520 }
03521 }
03522
03523
03524
03525 if(*nrhs <= 2) {
03526 j = 1;
03527 L70:
03528 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03529 if(*n > 1) {
03530 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
03531 / d__[*n - 1];
03532 }
03533 for(i__ = *n - 2; i__ >= 1; --i__) {
03534 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
03535 i__] * b_ref(i__ + 2, j)) / d__[i__];
03536
03537 }
03538 if(j < *nrhs) {
03539 ++j;
03540 goto L70;
03541 }
03542 } else {
03543 i__1 = *nrhs;
03544 for(j = 1; j <= i__1; ++j) {
03545 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03546 if(*n > 1) {
03547 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
03548 j)) / d__[*n - 1];
03549 }
03550 for(i__ = *n - 2; i__ >= 1; --i__) {
03551 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
03552 - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
03553
03554 }
03555
03556 }
03557 }
03558
03559 return 0;
03560
03561
03562
03563 }
03564 #endif
03565 #undef b_ref
03566 #endif
03567
03568