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 IC = (SUBNAM[0]);
01601 IZ = 'Z';
01602 if( IZ == 90 || IZ == 122 )
01603 {
01604
01605
01606
01607 if( IC >= 97 && IC <= 122 )
01608 {
01609 SUBNAM[0] = (char)(IC - 32);
01610 for( I=2; I <= 6; I++ )
01611 {
01612 IC = (SUBNAM[I-1]);
01613 if( IC >= 97 && IC <= 122 )
01614 SUBNAM[I - 1] = (char)(IC - 32);
01615 }
01616 }
01617
01618 }
01619 else if( IZ == 233 || IZ == 169 )
01620 {
01621
01622
01623
01624 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) ||
01625 (IC >= 162 && IC <= 169) )
01626 {
01627 SUBNAM[0] = (char)(IC + 64);
01628 for( I=2; I <= 6; I++ )
01629 {
01630 IC = (SUBNAM[I-1]);
01631 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <=
01632 153)) || (IC >= 162 && IC <= 169) )
01633 SUBNAM[I - 1] = (char)(IC + 64);
01634 }
01635 }
01636
01637 }
01638 else if( IZ == 218 || IZ == 250 )
01639 {
01640
01641
01642
01643 if( IC >= 225 && IC <= 250 )
01644 {
01645 SUBNAM[0] = (char)(IC - 32);
01646 for( I=2; I <= 6; I++ )
01647 {
01648 IC = (SUBNAM[I-1]);
01649 if( IC >= 225 && IC <= 250 )
01650 SUBNAM[I - 1] = (char)(IC - 32);
01651 }
01652 }
01653 }
01654
01655 C1 = SUBNAM[0];
01656 SNAME = C1 == 'S' || C1 == 'D';
01657 CNAME = C1 == 'C' || C1 == 'Z';
01658 if( !(CNAME || SNAME) )
01659 {
01660 return ILAENV_v;
01661 }
01662
01663 # if 0
01664 strncpy( C2, SUBNAM+1, 2 );
01665 strncpy( C3, SUBNAM+3, 3 );
01666 strncpy( C4, C3+1, 2 );
01667 # endif
01668
01669
01670
01671 strncpy( C2, SUBNAM+1, 2 );
01672 C2[2] = '\0';
01673 strncpy( C3, SUBNAM+3, 3 );
01674 C3[3] = '\0';
01675 strncpy( C4, C3+1, 2 );
01676 C4[2] = '\0';
01677
01678 switch( ISPEC )
01679 {
01680 case 1: goto L_110;
01681 case 2: goto L_200;
01682 case 3: goto L_300;
01683 }
01684
01685 L_110:
01686
01687
01688
01689
01690
01691
01692
01693 NB = 1;
01694
01695 if( strcmp(C2,"GE") == 0 )
01696 {
01697 if( strcmp(C3,"TRF") == 0 )
01698 {
01699 if( SNAME )
01700 {
01701 NB = 64;
01702 }
01703 else
01704 {
01705 NB = 64;
01706 }
01707 }
01708 else if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) ||
01709 strcmp(C3,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01710 {
01711 if( SNAME )
01712 {
01713 NB = 32;
01714 }
01715 else
01716 {
01717 NB = 32;
01718 }
01719 }
01720 else if( strcmp(C3,"HRD") == 0 )
01721 {
01722 if( SNAME )
01723 {
01724 NB = 32;
01725 }
01726 else
01727 {
01728 NB = 32;
01729 }
01730 }
01731 else if( strcmp(C3,"BRD") == 0 )
01732 {
01733 if( SNAME )
01734 {
01735 NB = 32;
01736 }
01737 else
01738 {
01739 NB = 32;
01740 }
01741 }
01742 else if( strcmp(C3,"TRI") == 0 )
01743 {
01744 if( SNAME )
01745 {
01746 NB = 64;
01747 }
01748 else
01749 {
01750 NB = 64;
01751 }
01752 }
01753 }
01754 else if( strcmp(C2,"PO") == 0 )
01755 {
01756 if( strcmp(C3,"TRF") == 0 )
01757 {
01758 if( SNAME )
01759 {
01760 NB = 64;
01761 }
01762 else
01763 {
01764 NB = 64;
01765 }
01766 }
01767 }
01768 else if( strcmp(C2,"SY") == 0 )
01769 {
01770 if( strcmp(C3,"TRF") == 0 )
01771 {
01772 if( SNAME )
01773 {
01774 NB = 64;
01775 }
01776 else
01777 {
01778 NB = 64;
01779 }
01780 }
01781 else if( SNAME && strcmp(C3,"TRD") == 0 )
01782 {
01783 NB = 1;
01784 }
01785 else if( SNAME && strcmp(C3,"GST") == 0 )
01786 {
01787 NB = 64;
01788 }
01789 }
01790 else if( CNAME && strcmp(C2,"HE") == 0 )
01791 {
01792 if( strcmp(C3,"TRF") == 0 )
01793 {
01794 NB = 64;
01795 }
01796 else if( strcmp(C3,"TRD") == 0 )
01797 {
01798 NB = 1;
01799 }
01800 else if( strcmp(C3,"GST") == 0 )
01801 {
01802 NB = 64;
01803 }
01804 }
01805 else if( SNAME && strcmp(C2,"OR") == 0 )
01806 {
01807 if( C3[0] == 'G' )
01808 {
01809 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01810 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01811 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01812 0 )
01813 {
01814 NB = 32;
01815 }
01816 }
01817 else if( C3[0] == 'M' )
01818 {
01819 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01820 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01821 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01822 0 )
01823 {
01824 NB = 32;
01825 }
01826 }
01827 }
01828 else if( CNAME && strcmp(C2,"UN") == 0 )
01829 {
01830 if( C3[0] == 'G' )
01831 {
01832 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01833 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01834 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01835 0 )
01836 {
01837 NB = 32;
01838 }
01839 }
01840 else if( C3[0] == 'M' )
01841 {
01842 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
01843 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
01844 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
01845 0 )
01846 {
01847 NB = 32;
01848 }
01849 }
01850 }
01851 else if( strcmp(C2,"GB") == 0 )
01852 {
01853 if( strcmp(C3,"TRF") == 0 )
01854 {
01855 if( SNAME )
01856 {
01857 if( N4 <= 64 )
01858 {
01859 NB = 1;
01860 }
01861 else
01862 {
01863 NB = 32;
01864 }
01865 }
01866 else
01867 {
01868 if( N4 <= 64 )
01869 {
01870 NB = 1;
01871 }
01872 else
01873 {
01874 NB = 32;
01875 }
01876 }
01877 }
01878 }
01879 else if( strcmp(C2,"PB") == 0 )
01880 {
01881 if( strcmp(C3,"TRF") == 0 )
01882 {
01883 if( SNAME )
01884 {
01885 if( N2 <= 64 )
01886 {
01887 NB = 1;
01888 }
01889 else
01890 {
01891 NB = 32;
01892 }
01893 }
01894 else
01895 {
01896 if( N2 <= 64 )
01897 {
01898 NB = 1;
01899 }
01900 else
01901 {
01902 NB = 32;
01903 }
01904 }
01905 }
01906 }
01907 else if( strcmp(C2,"TR") == 0 )
01908 {
01909 if( strcmp(C3,"TRI") == 0 )
01910 {
01911 if( SNAME )
01912 {
01913 NB = 64;
01914 }
01915 else
01916 {
01917 NB = 64;
01918 }
01919 }
01920 }
01921 else if( strcmp(C2,"LA") == 0 )
01922 {
01923 if( strcmp(C3,"UUM") == 0 )
01924 {
01925 if( SNAME )
01926 {
01927 NB = 64;
01928 }
01929 else
01930 {
01931 NB = 64;
01932 }
01933 }
01934 }
01935 else if( SNAME && strcmp(C2,"ST") == 0 )
01936 {
01937 if( strcmp(C3,"EBZ") == 0 )
01938 {
01939 NB = 1;
01940 }
01941 }
01942 ILAENV_v = NB;
01943 return ILAENV_v;
01944
01945 L_200:
01946
01947
01948
01949 NBMIN = 2;
01950 if( strcmp(C2,"GE") == 0 )
01951 {
01952 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
01953 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
01954 {
01955 if( SNAME )
01956 {
01957 NBMIN = 2;
01958 }
01959 else
01960 {
01961 NBMIN = 2;
01962 }
01963 }
01964 else if( strcmp(C3,"HRD") == 0 )
01965 {
01966 if( SNAME )
01967 {
01968 NBMIN = 2;
01969 }
01970 else
01971 {
01972 NBMIN = 2;
01973 }
01974 }
01975 else if( strcmp(C3,"BRD") == 0 )
01976 {
01977 if( SNAME )
01978 {
01979 NBMIN = 2;
01980 }
01981 else
01982 {
01983 NBMIN = 2;
01984 }
01985 }
01986 else if( strcmp(C3,"TRI") == 0 )
01987 {
01988 if( SNAME )
01989 {
01990 NBMIN = 2;
01991 }
01992 else
01993 {
01994 NBMIN = 2;
01995 }
01996 }
01997 }
01998 else if( strcmp(C2,"SY") == 0 )
01999 {
02000 if( strcmp(C3,"TRF") == 0 )
02001 {
02002 if( SNAME )
02003 {
02004 NBMIN = 8;
02005 }
02006 else
02007 {
02008 NBMIN = 8;
02009 }
02010 }
02011 else if( SNAME && strcmp(C3,"TRD") == 0 )
02012 {
02013 NBMIN = 2;
02014 }
02015 }
02016 else if( CNAME && strcmp(C2,"HE") == 0 )
02017 {
02018 if( strcmp(C3,"TRD") == 0 )
02019 {
02020 NBMIN = 2;
02021 }
02022 }
02023 else if( SNAME && strcmp(C2,"OR") == 0 )
02024 {
02025 if( C3[0] == 'G' )
02026 {
02027 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02028 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02029 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02030 0 )
02031 {
02032 NBMIN = 2;
02033 }
02034 }
02035 else if( C3[0] == 'M' )
02036 {
02037 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02038 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02039 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02040 0 )
02041 {
02042 NBMIN = 2;
02043 }
02044 }
02045 }
02046 else if( CNAME && strcmp(C2,"UN") == 0 )
02047 {
02048 if( C3[0] == 'G' )
02049 {
02050 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02051 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02052 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02053 0 )
02054 {
02055 NBMIN = 2;
02056 }
02057 }
02058 else if( C3[0] == 'M' )
02059 {
02060 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02061 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02062 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02063 0 )
02064 {
02065 NBMIN = 2;
02066 }
02067 }
02068 }
02069 ILAENV_v = NBMIN;
02070 return ILAENV_v;
02071
02072 L_300:
02073
02074
02075
02076 NX = 0;
02077 if( strcmp(C2,"GE") == 0 )
02078 {
02079 if( ((strcmp(C3,"QRF") == 0 || strcmp(C3,"RQF") == 0) || strcmp(C3
02080 ,"LQF") == 0) || strcmp(C3,"QLF") == 0 )
02081 {
02082 if( SNAME )
02083 {
02084 NX = 128;
02085 }
02086 else
02087 {
02088 NX = 128;
02089 }
02090 }
02091 else if( strcmp(C3,"HRD") == 0 )
02092 {
02093 if( SNAME )
02094 {
02095 NX = 128;
02096 }
02097 else
02098 {
02099 NX = 128;
02100 }
02101 }
02102 else if( strcmp(C3,"BRD") == 0 )
02103 {
02104 if( SNAME )
02105 {
02106 NX = 128;
02107 }
02108 else
02109 {
02110 NX = 128;
02111 }
02112 }
02113 }
02114 else if( strcmp(C2,"SY") == 0 )
02115 {
02116 if( SNAME && strcmp(C3,"TRD") == 0 )
02117 {
02118 NX = 1;
02119 }
02120 }
02121 else if( CNAME && strcmp(C2,"HE") == 0 )
02122 {
02123 if( strcmp(C3,"TRD") == 0 )
02124 {
02125 NX = 1;
02126 }
02127 }
02128 else if( SNAME && strcmp(C2,"OR") == 0 )
02129 {
02130 if( C3[0] == 'G' )
02131 {
02132 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02133 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02134 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02135 0 )
02136 {
02137 NX = 128;
02138 }
02139 }
02140 }
02141 else if( CNAME && strcmp(C2,"UN") == 0 )
02142 {
02143 if( C3[0] == 'G' )
02144 {
02145 if( (((((strcmp(C4,"QR") == 0 || strcmp(C4,"RQ") == 0) ||
02146 strcmp(C4,"LQ") == 0) || strcmp(C4,"QL") == 0) || strcmp(C4
02147 ,"HR") == 0) || strcmp(C4,"TR") == 0) || strcmp(C4,"BR") ==
02148 0 )
02149 {
02150 NX = 128;
02151 }
02152 }
02153 }
02154 ILAENV_v = NX;
02155 return ILAENV_v;
02156
02157 L_400:
02158
02159
02160
02161 ILAENV_v = 6;
02162 return ILAENV_v;
02163
02164 L_500:
02165
02166
02167
02168 ILAENV_v = 2;
02169 return ILAENV_v;
02170
02171 L_600:
02172
02173
02174
02175 ILAENV_v = (int32)((realnum)(MIN2(N1,N2))*1.6e0);
02176 return ILAENV_v;
02177
02178 L_700:
02179
02180
02181
02182 ILAENV_v = 1;
02183 return ILAENV_v;
02184
02185 L_800:
02186
02187
02188
02189 ILAENV_v = 50;
02190 return ILAENV_v;
02191
02192
02193
02194 }
02195
02196
02197
02198 STATIC void DSWAP(int32 n,
02199 double dx[],
02200 int32 incx,
02201 double dy[],
02202 int32 incy)
02203 {
02204 int32 i,
02205 ix,
02206 iy,
02207 m;
02208 double dtemp;
02209
02210 DEBUG_ENTRY( "DSWAP()" );
02211
02212
02213
02214
02215
02216
02217
02218 if( n <= 0 )
02219 {
02220 return;}
02221 if( incx == 1 && incy == 1 )
02222 goto L_20;
02223
02224
02225
02226
02227 ix = 1;
02228 iy = 1;
02229
02230 if( incx < 0 )
02231 ix = (-n + 1)*incx + 1;
02232
02233 if( incy < 0 )
02234 iy = (-n + 1)*incy + 1;
02235
02236 for( i=0; i < n; i++ )
02237 {
02238 dtemp = dx[ix-1];
02239 dx[ix-1] = dy[iy-1];
02240 dy[iy-1] = dtemp;
02241 ix += incx;
02242 iy += incy;
02243 }
02244 return;
02245
02246
02247
02248
02249
02250
02251 L_20:
02252 m = n%3;
02253 if( m == 0 )
02254 goto L_40;
02255
02256 for( i=0; i < m; i++ )
02257 {
02258 dtemp = dx[i];
02259 dx[i] = dy[i];
02260 dy[i] = dtemp;
02261 }
02262
02263 if( n < 3 )
02264 {
02265 return;
02266 }
02267
02268 L_40:
02269 for( i=m; i < n; i += 3 )
02270 {
02271 dtemp = dx[i];
02272 dx[i] = dy[i];
02273 dy[i] = dtemp;
02274 dtemp = dx[i+1];
02275 dx[i+1] = dy[i+1];
02276 dy[i+1] = dtemp;
02277 dtemp = dx[i+2];
02278 dx[i+2] = dy[i+2];
02279 dy[i+2] = dtemp;
02280 }
02281 return;
02282 }
02283
02284
02285
02286 STATIC void DSCAL(int32 n,
02287 double da,
02288 double dx[],
02289 int32 incx)
02290 {
02291 int32 i,
02292 m,
02293 nincx;
02294
02295 DEBUG_ENTRY( "DSCAL()" );
02296
02297
02298
02299
02300
02301
02302
02303
02304 if( n <= 0 || incx <= 0 )
02305 {
02306 return;}
02307 if( incx == 1 )
02308 goto L_20;
02309
02310
02311
02312 nincx = n*incx;
02313
02314 for( i=0; i<nincx; i = i + incx)
02315 {
02316 dx[i] *= da;
02317 }
02318 return;
02319
02320
02321
02322
02323
02324
02325 L_20:
02326 m = n%5;
02327 if( m == 0 )
02328 goto L_40;
02329
02330 for( i=0; i < m; i++ )
02331 {
02332 dx[i] *= da;
02333 }
02334
02335 if( n < 5 )
02336 {
02337 return;
02338 }
02339
02340 L_40:
02341
02342 for( i=m; i < n; i += 5 )
02343 {
02344 dx[i] *= da;
02345 dx[i+1] *= da;
02346 dx[i+2] *= da;
02347 dx[i+3] *= da;
02348 dx[i+4] *= da;
02349 }
02350 return;
02351 }
02352
02353
02354
02355 STATIC void DLASWP(int32 N,
02356 double *A,
02357 int32 LDA,
02358 int32 K1,
02359 int32 K2,
02360 int32 IPIV[],
02361 int32 INCX)
02362 {
02363 int32 I,
02364 IP,
02365 IX,
02366 I_;
02367
02368 DEBUG_ENTRY( "DLASWP()" );
02369
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 if( INCX == 0 )
02428 {
02429 return;}
02430 if( INCX > 0 )
02431 {
02432 IX = K1;
02433 }
02434 else
02435 {
02436 IX = 1 + (1 - K2)*INCX;
02437 }
02438 if( INCX == 1 )
02439 {
02440 for( I=K1; I <= K2; I++ )
02441 {
02442 I_ = I - 1;
02443 IP = IPIV[I_];
02444 if( IP != I )
02445 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02446 }
02447 }
02448 else if( INCX > 1 )
02449 {
02450 for( I=K1; I <= K2; I++ )
02451 {
02452 I_ = I - 1;
02453 IP = IPIV[IX-1];
02454 if( IP != I )
02455 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02456 IX += INCX;
02457 }
02458 }
02459 else if( INCX < 0 )
02460 {
02461 for( I=K2; I >= K1; I-- )
02462 {
02463 I_ = I - 1;
02464 IP = IPIV[IX-1];
02465 if( IP != I )
02466 DSWAP(N,&AA(0,I_),LDA,&AA(0,IP-1),LDA);
02467 IX += INCX;
02468 }
02469 }
02470
02471 return;
02472
02473
02474
02475 #undef A
02476 }
02477
02478
02479
02480 STATIC void DGETF2(int32 M,
02481 int32 N,
02482 double *A,
02483 int32 LDA,
02484 int32 IPIV[],
02485 int32 *INFO)
02486 {
02487 int32 J,
02488 JP,
02489 J_,
02490 limit;
02491
02492 DEBUG_ENTRY( "DGETF2()" );
02493
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 *INFO = 0;
02564 if( M < 0 )
02565 {
02566 *INFO = -1;
02567 }
02568 else if( N < 0 )
02569 {
02570 *INFO = -2;
02571 }
02572 else if( LDA < MAX2(1,M) )
02573 {
02574 *INFO = -4;
02575 }
02576 if( *INFO != 0 )
02577 {
02578 XERBLA("DGETF2",-*INFO);
02579
02580 }
02581
02582
02583
02584 if( M == 0 || N == 0 )
02585 {
02586 return;}
02587
02588 limit = MIN2(M,N);
02589 for( J=1; J <= limit; J++ )
02590 {
02591 J_ = J - 1;
02592
02593
02594
02595 JP = J - 1 + IDAMAX(M-J+1,&AA(J_,J_),1);
02596 IPIV[J_] = JP;
02597 if( AA(J_,JP-1) != ZERO )
02598 {
02599
02600
02601 if( JP != J )
02602 DSWAP(N,&AA(0,J_),LDA,&AA(0,JP-1),LDA);
02603
02604
02605
02606 if( J < M )
02607 DSCAL(M-J,ONE/AA(J_,J_),&AA(J_,J_+1),1);
02608 }
02609 else if( *INFO == 0 )
02610 {
02611 *INFO = J;
02612 }
02613
02614 if( J < MIN2(M,N) )
02615 {
02616
02617
02618 DGER(M-J,N-J,-ONE,&AA(J_,J_+1),1,&AA(J_+1,J_),LDA,&AA(J_+1,J_+1),
02619 LDA);
02620 }
02621 }
02622 return;
02623
02624
02625
02626 #undef A
02627 }
02628
02629
02630
02631 STATIC void DGER(int32 M,
02632 int32 N,
02633 double ALPHA,
02634 double X[],
02635 int32 INCX,
02636 double Y[],
02637 int32 INCY,
02638 double *A,
02639 int32 LDA)
02640 {
02641 int32 I,
02642 INFO,
02643 IX,
02644 I_,
02645 J,
02646 JY,
02647 J_,
02648 KX;
02649 double TEMP;
02650
02651 DEBUG_ENTRY( "DGER()" );
02652
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 INFO = 0;
02736 if( M < 0 )
02737 {
02738 INFO = 1;
02739 }
02740 else if( N < 0 )
02741 {
02742 INFO = 2;
02743 }
02744 else if( INCX == 0 )
02745 {
02746 INFO = 5;
02747 }
02748 else if( INCY == 0 )
02749 {
02750 INFO = 7;
02751 }
02752 else if( LDA < MAX2(1,M) )
02753 {
02754 INFO = 9;
02755 }
02756 if( INFO != 0 )
02757 {
02758 XERBLA("DGER ",INFO);
02759
02760 }
02761
02762
02763
02764 if( ((M == 0) || (N == 0)) || (ALPHA == ZERO) )
02765 {
02766 return;}
02767
02768
02769
02770
02771 if( INCY > 0 )
02772 {
02773 JY = 1;
02774 }
02775 else
02776 {
02777 JY = 1 - (N - 1)*INCY;
02778 }
02779 if( INCX == 1 )
02780 {
02781 for( J=1; J <= N; J++ )
02782 {
02783 J_ = J - 1;
02784 if( Y[JY-1] != ZERO )
02785 {
02786 TEMP = ALPHA*Y[JY-1];
02787 for( I=1; I <= M; I++ )
02788 {
02789 I_ = I - 1;
02790 AA(J_,I_) += X[I_]*TEMP;
02791 }
02792 }
02793 JY += INCY;
02794 }
02795 }
02796 else
02797 {
02798 if( INCX > 0 )
02799 {
02800 KX = 1;
02801 }
02802 else
02803 {
02804 KX = 1 - (M - 1)*INCX;
02805 }
02806 for( J=1; J <= N; J++ )
02807 {
02808 J_ = J - 1;
02809 if( Y[JY-1] != ZERO )
02810 {
02811 TEMP = ALPHA*Y[JY-1];
02812 IX = KX;
02813 for( I=1; I <= M; I++ )
02814 {
02815 I_ = I - 1;
02816 AA(J_,I_) += X[IX-1]*TEMP;
02817 IX += INCX;
02818 }
02819 }
02820 JY += INCY;
02821 }
02822 }
02823
02824 return;
02825
02826
02827
02828 #undef A
02829 }
02830
02831
02832
02833 STATIC void DGEMM(int32 TRANSA,
02834 int32 TRANSB,
02835 int32 M,
02836 int32 N,
02837 int32 K,
02838 double ALPHA,
02839 double *A,
02840 int32 LDA,
02841 double *B,
02842 int32 LDB,
02843 double BETA,
02844 double *C,
02845 int32 LDC)
02846 {
02847 int32 NOTA,
02848 NOTB;
02849 int32 I,
02850 INFO,
02851 I_,
02852 J,
02853 J_,
02854 L,
02855 L_,
02856 NROWA,
02857 NROWB;
02858 double TEMP;
02859
02860 DEBUG_ENTRY( "DGEMM()" );
02861
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 NOTA = LSAME(TRANSA,'N');
02998 NOTB = LSAME(TRANSB,'N');
02999
03000 if( NOTA )
03001 {
03002 NROWA = M;
03003 }
03004 else
03005 {
03006 NROWA = K;
03007 }
03008
03009 if( NOTB )
03010 {
03011 NROWB = K;
03012 }
03013 else
03014 {
03015 NROWB = N;
03016 }
03017
03018
03019
03020 INFO = 0;
03021 if( ((!NOTA) && (!LSAME(TRANSA,'C'))) && (!LSAME(TRANSA,'T')) )
03022 {
03023 INFO = 1;
03024 }
03025 else if(
03026 ((!NOTB) && (!LSAME(TRANSB,'C'))) && (!LSAME(TRANSB,'T')) )
03027 {
03028 INFO = 2;
03029 }
03030
03031 else if( M < 0 )
03032 {
03033 INFO = 3;
03034 }
03035
03036 else if( N < 0 )
03037 {
03038 INFO = 4;
03039 }
03040
03041 else if( K < 0 )
03042 {
03043 INFO = 5;
03044 }
03045
03046 else if( LDA < MAX2(1,NROWA) )
03047 {
03048 INFO = 8;
03049 }
03050
03051 else if( LDB < MAX2(1,NROWB) )
03052 {
03053 INFO = 10;
03054 }
03055
03056 else if( LDC < MAX2(1,M) )
03057 {
03058 INFO = 13;
03059 }
03060
03061 if( INFO != 0 )
03062 {
03063 XERBLA("DGEMM ",INFO);
03064
03065 }
03066
03067
03068
03069 if( ((M == 0) || (N == 0)) || (((ALPHA == ZERO) || (K == 0)) &&
03070 (BETA == ONE)) )
03071 {
03072 return;
03073 }
03074
03075
03076 if( ALPHA == ZERO )
03077 {
03078 if( BETA == ZERO )
03079 {
03080 for( J=1; J <= N; J++ )
03081 {
03082 J_ = J - 1;
03083 for( I=1; I <= M; I++ )
03084 {
03085 I_ = I - 1;
03086 CC(J_,I_) = ZERO;
03087 }
03088 }
03089 }
03090
03091 else
03092 {
03093 for( J=1; J <= N; J++ )
03094 {
03095 J_ = J - 1;
03096 for( I=1; I <= M; I++ )
03097 {
03098 I_ = I - 1;
03099 CC(J_,I_) *= BETA;
03100 }
03101 }
03102 }
03103 return;
03104 }
03105
03106
03107 if( NOTB )
03108 {
03109
03110 if( NOTA )
03111 {
03112
03113
03114
03115 for( J=1; J <= N; J++ )
03116 {
03117 J_ = J - 1;
03118 if( BETA == ZERO )
03119 {
03120 for( I=1; I <= M; I++ )
03121 {
03122 I_ = I - 1;
03123 CC(J_,I_) = ZERO;
03124 }
03125 }
03126
03127 else if( BETA != ONE )
03128 {
03129 for( I=1; I <= M; I++ )
03130 {
03131 I_ = I - 1;
03132 CC(J_,I_) *= BETA;
03133 }
03134 }
03135
03136 for( L=1; L <= K; L++ )
03137 {
03138 L_ = L - 1;
03139 if( BB(J_,L_) != ZERO )
03140 {
03141 TEMP = ALPHA*BB(J_,L_);
03142 for( I=1; I <= M; I++ )
03143 {
03144 I_ = I - 1;
03145 CC(J_,I_) += TEMP*AA(L_,I_);
03146 }
03147 }
03148 }
03149 }
03150 }
03151 else
03152 {
03153
03154
03155 for( J=1; J <= N; J++ )
03156 {
03157 J_ = J - 1;
03158 for( I=1; I <= M; I++ )
03159 {
03160 I_ = I - 1;
03161 TEMP = ZERO;
03162 for( L=1; L <= K; L++ )
03163 {
03164 L_ = L - 1;
03165 TEMP += AA(I_,L_)*BB(J_,L_);
03166 }
03167
03168 if( BETA == ZERO )
03169 {
03170 CC(J_,I_) = ALPHA*TEMP;
03171 }
03172 else
03173 {
03174 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03175 }
03176 }
03177 }
03178 }
03179 }
03180 else
03181 {
03182 if( NOTA )
03183 {
03184
03185
03186
03187 for( J=1; J <= N; J++ )
03188 {
03189 J_ = J - 1;
03190 if( BETA == ZERO )
03191 {
03192 for( I=1; I <= M; I++ )
03193 {
03194 I_ = I - 1;
03195 CC(J_,I_) = ZERO;
03196 }
03197 }
03198
03199 else if( BETA != ONE )
03200 {
03201 for( I=1; I <= M; I++ )
03202 {
03203 I_ = I - 1;
03204 CC(J_,I_) *= BETA;
03205 }
03206 }
03207
03208 for( L=1; L <= K; L++ )
03209 {
03210 L_ = L - 1;
03211 if( BB(L_,J_) != ZERO )
03212 {
03213 TEMP = ALPHA*BB(L_,J_);
03214 for( I=1; I <= M; I++ )
03215 {
03216 I_ = I - 1;
03217 CC(J_,I_) += TEMP*AA(L_,I_);
03218 }
03219 }
03220 }
03221 }
03222 }
03223
03224 else
03225 {
03226
03227
03228 for( J=1; J <= N; J++ )
03229 {
03230 J_ = J - 1;
03231
03232 for( I=1; I <= M; I++ )
03233 {
03234 I_ = I - 1;
03235 TEMP = ZERO;
03236
03237 for( L=1; L <= K; L++ )
03238 {
03239 L_ = L - 1;
03240 TEMP += AA(I_,L_)*BB(L_,J_);
03241 }
03242
03243 if( BETA == ZERO )
03244 {
03245 CC(J_,I_) = ALPHA*TEMP;
03246 }
03247 else
03248 {
03249 CC(J_,I_) = ALPHA*TEMP + BETA*CC(J_,I_);
03250 }
03251
03252 }
03253 }
03254 }
03255 }
03256
03257 return;
03258
03259
03260 #undef C
03261 #undef B
03262 #undef A
03263 }
03264 #undef AA
03265 #undef BB
03266 #undef CC
03267
03268
03269 #if 0
03270 STATIC int32 DGTSV(int32 *n, int32 *nrhs, double *dl,
03271 double *d__, double *du, double *b, int32 *ldb, int32
03272 *info)
03273 {
03274
03275 int32 b_dim1, b_offset, i__1, i__2;
03276
03277
03278 double fact, temp;
03279 int32 i__, j;
03280
03281
03282 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
03283
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 --dl;
03353 --d__;
03354 --du;
03355 b_dim1 = *ldb;
03356 b_offset = 1 + b_dim1 * 1;
03357 b -= b_offset;
03358
03359
03360 *info = 0;
03361 if(*n < 0) {
03362 *info = -1;
03363 } else if(*nrhs < 0) {
03364 *info = -2;
03365 } else if(*ldb < *n && *ldb < 1) {
03366 *info = -7;
03367 }
03368 if(*info != 0) {
03369 i__1 = -(*info);
03370 XERBLA("DGTSV ", i__1);
03371
03372 }
03373
03374 if(*n == 0) {
03375 return 0;
03376 }
03377
03378 if(*nrhs == 1) {
03379 i__1 = *n - 2;
03380 for(i__ = 1; i__ <= i__1; ++i__) {
03381 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03382
03383
03384
03385 if(d__[i__] != 0.) {
03386 fact = dl[i__] / d__[i__];
03387 d__[i__ + 1] -= fact * du[i__];
03388 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
03389 1);
03390 } else {
03391 *info = i__;
03392 return 0;
03393 }
03394 dl[i__] = 0.;
03395 } else {
03396
03397
03398
03399 fact = d__[i__] / dl[i__];
03400 d__[i__] = dl[i__];
03401 temp = d__[i__ + 1];
03402 d__[i__ + 1] = du[i__] - fact * temp;
03403 dl[i__] = du[i__ + 1];
03404 du[i__ + 1] = -fact * dl[i__];
03405 du[i__] = temp;
03406 temp = b_ref(i__, 1);
03407 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03408 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03409 }
03410
03411 }
03412 if(*n > 1) {
03413 i__ = *n - 1;
03414 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03415 if(d__[i__] != 0.) {
03416 fact = dl[i__] / d__[i__];
03417 d__[i__ + 1] -= fact * du[i__];
03418 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
03419 1);
03420 } else {
03421 *info = i__;
03422 return 0;
03423 }
03424 } else {
03425 fact = d__[i__] / dl[i__];
03426 d__[i__] = dl[i__];
03427 temp = d__[i__ + 1];
03428 d__[i__ + 1] = du[i__] - fact * temp;
03429 du[i__] = temp;
03430 temp = b_ref(i__, 1);
03431 b_ref(i__, 1) = b_ref(i__ + 1, 1);
03432 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
03433 }
03434 }
03435 if(d__[*n] == 0.) {
03436 *info = *n;
03437 return 0;
03438 }
03439 } else {
03440 i__1 = *n - 2;
03441 for(i__ = 1; i__ <= i__1; ++i__) {
03442 if(fabs(d__[i__]) >= fabs(dl[i__])) {
03443
03444
03445
03446 if(d__[i__] != 0.) {
03447 fact = dl[i__] / d__[i__];
03448 d__[i__ + 1] -= fact * du[i__];
03449 i__2 = *nrhs;
03450 for(j = 1; j <= i__2; ++j) {
03451 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03452 i__, j);
03453
03454 }
03455 } else {
03456 *info = i__;
03457 return 0;
03458 }
03459 dl[i__] = 0.;
03460 } else {
03461
03462
03463
03464 fact = d__[i__] / dl[i__];
03465 d__[i__] = dl[i__];
03466 temp = d__[i__ + 1];
03467 d__[i__ + 1] = du[i__] - fact * temp;
03468 dl[i__] = du[i__ + 1];
03469 du[i__ + 1] = -fact * dl[i__];
03470 du[i__] = temp;
03471 i__2 = *nrhs;
03472 for(j = 1; j <= i__2; ++j) {
03473 temp = b_ref(i__, j);
03474 b_ref(i__, j) = b_ref(i__ + 1, j);
03475 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03476
03477 }
03478 }
03479
03480 }
03481 if(*n > 1) {
03482 i__ = *n - 1;
03483 if( fabs(d__[i__]) >= fabs(dl[i__]))
03484 {
03485 if(d__[i__] != 0.)
03486 {
03487 fact = dl[i__] / d__[i__];
03488 d__[i__ + 1] -= fact * du[i__];
03489 i__1 = *nrhs;
03490 for(j = 1; j <= i__1; ++j) {
03491 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
03492 i__, j);
03493
03494 }
03495 }
03496 else
03497 {
03498 *info = i__;
03499 return 0;
03500 }
03501 } else {
03502 fact = d__[i__] / dl[i__];
03503 d__[i__] = dl[i__];
03504 temp = d__[i__ + 1];
03505 d__[i__ + 1] = du[i__] - fact * temp;
03506 du[i__] = temp;
03507 i__1 = *nrhs;
03508 for(j = 1; j <= i__1; ++j) {
03509 temp = b_ref(i__, j);
03510 b_ref(i__, j) = b_ref(i__ + 1, j);
03511 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
03512
03513 }
03514 }
03515 }
03516 if(d__[*n] == 0.) {
03517 *info = *n;
03518 return 0;
03519 }
03520 }
03521
03522
03523
03524 if(*nrhs <= 2) {
03525 j = 1;
03526 L70:
03527 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03528 if(*n > 1) {
03529 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
03530 / d__[*n - 1];
03531 }
03532 for(i__ = *n - 2; i__ >= 1; --i__) {
03533 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
03534 i__] * b_ref(i__ + 2, j)) / d__[i__];
03535
03536 }
03537 if(j < *nrhs) {
03538 ++j;
03539 goto L70;
03540 }
03541 } else {
03542 i__1 = *nrhs;
03543 for(j = 1; j <= i__1; ++j) {
03544 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
03545 if(*n > 1) {
03546 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
03547 j)) / d__[*n - 1];
03548 }
03549 for(i__ = *n - 2; i__ >= 1; --i__) {
03550 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
03551 - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];
03552
03553 }
03554
03555 }
03556 }
03557
03558 return 0;
03559
03560
03561
03562 }
03563 #endif
03564 #undef b_ref
03565 #endif
03566
03567