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