20 extern "C" void dgetrf_(int32 *
M, int32 *N,
double *A, int32 *LDA, int32 *IPIV, int32 *INFO);
21 extern "C" void dgetrs_(
char *TRANS, int32 *N, int32 *NRHS,
double *A, int32 *LDA, int32 *iPiv,
double *B,
22 int32 *LDB, int32 *INFO, int32 translen);
23 extern "C" void dgtsv_(int32 *n, int32 *nrhs,
double *dl,
double *d,
double *du,
double *b, int32 *ldb, int32 *info);
34 STATIC void DGETRF(int32,int32,
double*,int32,int32[],int32*);
37 STATIC void DGETRS(int32 TRANS,int32 N,int32 NRHS,
double *A,int32 LDA,int32 IPIV[],
double *B,int32 LDB,int32 *INFO);
51 int32 M_loc, N_loc, lda_loc;
61 dgetrf_(&M_loc, &N_loc, A , &lda_loc, ipiv, info);
64 DGETRF(M_loc, N_loc, A, lda_loc, ipiv, info);
70 double *B,
long ldb, int32 *info)
74 int32 N_loc, nrhs_loc, lda_loc, ldb_loc;
79 nrhs_loc = (int32)nrhs;
85 dgetrs_(&trans, &N_loc, &nrhs_loc, A, &lda_loc, ipiv, B, &ldb_loc, info,
sizeof(
char));
88 DGETRS(trans, N_loc, nrhs_loc, A, lda_loc, ipiv, B, ldb_loc, info);
94 void dgtsv_wrapper(
long N,
long nrhs,
double *dl,
double *d__,
double *du,
double *b,
long ldb, int32 *info)
96 printf(
"Inside dgtsv\n");
100 int32 N_loc, nrhs_loc, ldb_loc;
105 nrhs_loc = (int32)nrhs;
106 ldb_loc = (int32)ldb;
110 dgtsv_(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
114 (void)DGTSV(&N_loc, &nrhs_loc, dl, d__, du, b, &ldb_loc, info);
136 #define AA(I_,J_) (*(A+(I_)*(LDA)+(J_)))
137 #define BB(I_,J_) (*(B+(I_)*(LDB)+(J_)))
138 #define CC(I_,J_) (*(C+(I_)*(LDC)+(J_)))
289 fprintf(
ioQQQ,
" ** On entry to %6.6s parameter number %2ld had an illegal value\n",
290 SRNAME, (
long)INFO );
403 else if( LDA <
MAX2(1,M) )
415 if( M == 0 || N == 0 )
423 NB =
ILAENV(1,
"DGETRF",M,N,-1);
424 if( NB <= 1 || NB >=
MIN2(M,N) )
428 DGETF2(M,N,A,LDA,IPIV,INFO);
438 for( J=1; J<=limit; J += NB )
446 DGETF2(M-J+1,JB,&
AA(J_,J_),LDA,&IPIV[J_],&IINFO);
450 if( *INFO == 0 && IINFO > 0 )
451 *INFO = IINFO + J - 1;
452 limit2 =
MIN2(M,J+JB-1);
453 for( I=J; I <= limit2; I++ )
461 DLASWP(J-1,A,LDA,J,J+JB-1,IPIV,1);
468 DLASWP(N-J-JB+1,&
AA(J_+JB,0),LDA,J,J+JB-1,IPIV,1);
477 DTRSM(chL1,chL2,chL3,chL4,JB,N-J-JB+1,
ONE,&
AA(J_,J_),
478 LDA,&
AA(J_+JB,J_),LDA);
487 DGEMM(chL1,chL2,M-J-JB+1,N-J-JB+1,JB,-
ONE,&
AA(J_,J_+JB),
488 LDA,&
AA(J_+JB,J_),LDA,
ONE,&
AA(J_+JB,J_+JB),LDA);
626 NOTRAN =
LSAME(TRANS,
'N');
627 if( (!NOTRAN && !
LSAME(TRANS,
'T')) && !
LSAME(TRANS,
'C') )
639 else if( LDA <
MAX2(1,N) )
643 else if( LDB <
MAX2(1,N) )
655 if( N == 0 || NRHS == 0 )
667 DLASWP(NRHS,B,LDB,1,N,IPIV,1);
676 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,
ONE,A,LDA,B,LDB);
685 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,
ONE,A,LDA,B,LDB);
699 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,
ONE,A,LDA,B,LDB);
708 DTRSM(chL1,chL2,chL3,chL4,N,NRHS,
ONE,A,LDA,B,LDB);
712 DLASWP(NRHS,B,LDB,1,N,IPIV,-1);
786 if( ZCODE == 90 || ZCODE == 122 )
792 if( INTA >= 97 && INTA <= 122 )
794 if( INTB >= 97 && INTB <= 122 )
798 else if( ZCODE == 233 || ZCODE == 169 )
804 if( ((INTA >= 129 && INTA <= 137) || (INTA >= 145 && INTA <=
805 153)) || (INTA >= 162 && INTA <= 169) )
807 if( ((INTB >= 129 && INTB <= 137) || (INTB >= 145 && INTB <=
808 153)) || (INTB >= 162 && INTB <= 169) )
812 else if( ZCODE == 218 || ZCODE == 250 )
818 if( INTA >= 225 && INTA <= 250 )
820 if( INTB >= 225 && INTB <= 250 )
823 LSAME_v = INTA == INTB;
854 if( n < 1 || incx <= 0 )
874 for( i=2; i <= n; i++ )
877 if( fabs(dx[ix-1]) > dmax )
880 dmax = fabs(dx[ix-1]);
890 for( i=1; i < n; i++ )
893 if( fabs(dx[i]) > dmax )
1062 LSIDE =
LSAME(SIDE,
'L');
1071 NOUNIT =
LSAME(DIAG,
'N');
1072 UPPER =
LSAME(UPLO,
'U');
1075 if( (!LSIDE) && (!
LSAME(SIDE,
'R')) )
1079 else if( (!UPPER) && (!
LSAME(UPLO,
'L')) )
1083 else if( ((!
LSAME(TRANSA,
'N')) && (!
LSAME(TRANSA,
'T'))) && (!
LSAME(TRANSA,
1088 else if( (!
LSAME(DIAG,
'U')) && (!
LSAME(DIAG,
'N')) )
1100 else if( LDA <
MAX2(1,NROWA) )
1104 else if( LDB <
MAX2(1,M) )
1124 for( J=1; J <= N; J++ )
1127 for( I=1; I <=
M; I++ )
1140 if(
LSAME(TRANSA,
'N') )
1147 for( J=1; J <= N; J++ )
1152 for( I=1; I <=
M; I++ )
1158 for( K=M; K >= 1; K-- )
1164 BB(J_,K_) /=
AA(K_,K_);
1165 for( I=1; I <= (K - 1); I++ )
1168 BB(J_,I_) += -
BB(J_,K_)*
AA(K_,I_);
1176 for( J=1; J <= N; J++ )
1181 for( I=1; I <=
M; I++ )
1187 for( K=1; K <=
M; K++ )
1193 BB(J_,K_) /=
AA(K_,K_);
1194 for( I=K + 1; I <=
M; I++ )
1197 BB(J_,I_) += -
BB(J_,K_)*
AA(K_,I_);
1211 for( J=1; J <= N; J++ )
1214 for( I=1; I <=
M; I++ )
1217 TEMP = ALPHA*
BB(J_,I_);
1218 for( K=1; K <= (I - 1); K++ )
1221 TEMP += -
AA(I_,K_)*
BB(J_,K_);
1231 for( J=1; J <= N; J++ )
1234 for( I=M; I >= 1; I-- )
1237 TEMP = ALPHA*
BB(J_,I_);
1238 for( K=I + 1; K <=
M; K++ )
1241 TEMP += -
AA(I_,K_)*
BB(J_,K_);
1253 if(
LSAME(TRANSA,
'N') )
1260 for( J=1; J <= N; J++ )
1265 for( I=1; I <=
M; I++ )
1271 for( K=1; K <= (J - 1); K++ )
1276 for( I=1; I <=
M; I++ )
1279 BB(J_,I_) += -
AA(J_,K_)*
BB(K_,I_);
1285 TEMP =
ONE/
AA(J_,J_);
1286 for( I=1; I <=
M; I++ )
1296 for( J=N; J >= 1; J-- )
1301 for( I=1; I <=
M; I++ )
1307 for( K=J + 1; K <= N; K++ )
1312 for( I=1; I <=
M; I++ )
1315 BB(J_,I_) += -
AA(J_,K_)*
BB(K_,I_);
1321 TEMP =
ONE/
AA(J_,J_);
1322 for( I=1; I <=
M; I++ )
1338 for( K=N; K >= 1; K-- )
1343 TEMP =
ONE/
AA(K_,K_);
1344 for( I=1; I <=
M; I++ )
1350 for( J=1; J <= (K - 1); J++ )
1356 for( I=1; I <=
M; I++ )
1359 BB(J_,I_) += -TEMP*
BB(K_,I_);
1365 for( I=1; I <=
M; I++ )
1375 for( K=1; K <= N; K++ )
1380 TEMP =
ONE/
AA(K_,K_);
1381 for( I=1; I <=
M; I++ )
1387 for( J=K + 1; J <= N; J++ )
1393 for( I=1; I <=
M; I++ )
1396 BB(J_,I_) += -TEMP*
BB(K_,I_);
1402 for( I=1; I <=
M; I++ )
1599 strncpy( SUBNAM, NAME, 6 );
1602 if( IZ == 90 || IZ == 122 )
1607 if( IC >= 97 && IC <= 122 )
1609 SUBNAM[0] = (char)(IC - 32);
1610 for( I=2; I <= 6; I++ )
1613 if( IC >= 97 && IC <= 122 )
1614 SUBNAM[I - 1] = (char)(IC - 32);
1619 else if( IZ == 233 || IZ == 169 )
1624 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <= 153)) ||
1625 (IC >= 162 && IC <= 169) )
1627 SUBNAM[0] = (char)(IC + 64);
1628 for( I=2; I <= 6; I++ )
1631 if( ((IC >= 129 && IC <= 137) || (IC >= 145 && IC <=
1632 153)) || (IC >= 162 && IC <= 169) )
1633 SUBNAM[I - 1] = (char)(IC + 64);
1638 else if( IZ == 218 || IZ == 250 )
1643 if( IC >= 225 && IC <= 250 )
1645 SUBNAM[0] = (char)(IC - 32);
1646 for( I=2; I <= 6; I++ )
1649 if( IC >= 225 && IC <= 250 )
1650 SUBNAM[I - 1] = (char)(IC - 32);
1656 SNAME = C1 ==
'S' || C1 ==
'D';
1657 CNAME = C1 ==
'C' || C1 ==
'Z';
1658 if( !(CNAME || SNAME) )
1664 strncpy( C2, SUBNAM+1, 2 );
1665 strncpy( C3, SUBNAM+3, 3 );
1666 strncpy( C4, C3+1, 2 );
1671 strncpy( C2, SUBNAM+1, 2 );
1673 strncpy( C3, SUBNAM+3, 3 );
1675 strncpy( C4, C3+1, 2 );
1695 if( strcmp(C2,
"GE") == 0 )
1697 if( strcmp(C3,
"TRF") == 0 )
1708 else if( ((strcmp(C3,
"QRF") == 0 || strcmp(C3,
"RQF") == 0) ||
1709 strcmp(C3,
"LQF") == 0) || strcmp(C3,
"QLF") == 0 )
1720 else if( strcmp(C3,
"HRD") == 0 )
1731 else if( strcmp(C3,
"BRD") == 0 )
1742 else if( strcmp(C3,
"TRI") == 0 )
1754 else if( strcmp(C2,
"PO") == 0 )
1756 if( strcmp(C3,
"TRF") == 0 )
1768 else if( strcmp(C2,
"SY") == 0 )
1770 if( strcmp(C3,
"TRF") == 0 )
1781 else if( SNAME && strcmp(C3,
"TRD") == 0 )
1785 else if( SNAME && strcmp(C3,
"GST") == 0 )
1790 else if( CNAME && strcmp(C2,
"HE") == 0 )
1792 if( strcmp(C3,
"TRF") == 0 )
1796 else if( strcmp(C3,
"TRD") == 0 )
1800 else if( strcmp(C3,
"GST") == 0 )
1805 else if( SNAME && strcmp(C2,
"OR") == 0 )
1809 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
1810 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
1811 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
1817 else if( C3[0] ==
'M' )
1819 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
1820 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
1821 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
1828 else if( CNAME && strcmp(C2,
"UN") == 0 )
1832 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
1833 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
1834 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
1840 else if( C3[0] ==
'M' )
1842 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
1843 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
1844 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
1851 else if( strcmp(C2,
"GB") == 0 )
1853 if( strcmp(C3,
"TRF") == 0 )
1879 else if( strcmp(C2,
"PB") == 0 )
1881 if( strcmp(C3,
"TRF") == 0 )
1907 else if( strcmp(C2,
"TR") == 0 )
1909 if( strcmp(C3,
"TRI") == 0 )
1921 else if( strcmp(C2,
"LA") == 0 )
1923 if( strcmp(C3,
"UUM") == 0 )
1935 else if( SNAME && strcmp(C2,
"ST") == 0 )
1937 if( strcmp(C3,
"EBZ") == 0 )
1950 if( strcmp(C2,
"GE") == 0 )
1952 if( ((strcmp(C3,
"QRF") == 0 || strcmp(C3,
"RQF") == 0) || strcmp(C3
1953 ,
"LQF") == 0) || strcmp(C3,
"QLF") == 0 )
1964 else if( strcmp(C3,
"HRD") == 0 )
1975 else if( strcmp(C3,
"BRD") == 0 )
1986 else if( strcmp(C3,
"TRI") == 0 )
1998 else if( strcmp(C2,
"SY") == 0 )
2000 if( strcmp(C3,
"TRF") == 0 )
2011 else if( SNAME && strcmp(C3,
"TRD") == 0 )
2016 else if( CNAME && strcmp(C2,
"HE") == 0 )
2018 if( strcmp(C3,
"TRD") == 0 )
2023 else if( SNAME && strcmp(C2,
"OR") == 0 )
2027 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2028 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2029 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2035 else if( C3[0] ==
'M' )
2037 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2038 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2039 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2046 else if( CNAME && strcmp(C2,
"UN") == 0 )
2050 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2051 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2052 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2058 else if( C3[0] ==
'M' )
2060 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2061 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2062 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2077 if( strcmp(C2,
"GE") == 0 )
2079 if( ((strcmp(C3,
"QRF") == 0 || strcmp(C3,
"RQF") == 0) || strcmp(C3
2080 ,
"LQF") == 0) || strcmp(C3,
"QLF") == 0 )
2091 else if( strcmp(C3,
"HRD") == 0 )
2102 else if( strcmp(C3,
"BRD") == 0 )
2114 else if( strcmp(C2,
"SY") == 0 )
2116 if( SNAME && strcmp(C3,
"TRD") == 0 )
2121 else if( CNAME && strcmp(C2,
"HE") == 0 )
2123 if( strcmp(C3,
"TRD") == 0 )
2128 else if( SNAME && strcmp(C2,
"OR") == 0 )
2132 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2133 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2134 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2141 else if( CNAME && strcmp(C2,
"UN") == 0 )
2145 if( (((((strcmp(C4,
"QR") == 0 || strcmp(C4,
"RQ") == 0) ||
2146 strcmp(C4,
"LQ") == 0) || strcmp(C4,
"QL") == 0) || strcmp(C4
2147 ,
"HR") == 0) || strcmp(C4,
"TR") == 0) || strcmp(C4,
"BR") ==
2221 if( incx == 1 && incy == 1 )
2231 ix = (-n + 1)*incx + 1;
2234 iy = (-n + 1)*incy + 1;
2236 for( i=0; i < n; i++ )
2239 dx[ix-1] = dy[iy-1];
2256 for( i=0; i < m; i++ )
2269 for( i=m; i < n; i += 3 )
2304 if( n <= 0 || incx <= 0 )
2314 for( i=0; i<nincx; i = i + incx)
2330 for( i=0; i < m; i++ )
2342 for( i=m; i < n; i += 5 )
2436 IX = 1 + (1 - K2)*INCX;
2440 for( I=K1; I <= K2; I++ )
2450 for( I=K1; I <= K2; I++ )
2461 for( I=K2; I >= K1; I-- )
2572 else if( LDA <
MAX2(1,M) )
2584 if( M == 0 || N == 0 )
2589 for( J=1; J <= limit; J++ )
2595 JP = J - 1 +
IDAMAX(M-J+1,&
AA(J_,J_),1);
2597 if(
AA(J_,JP-1) !=
ZERO )
2609 else if( *INFO == 0 )
2618 DGER(M-J,N-J,-
ONE,&
AA(J_,J_+1),1,&
AA(J_+1,J_),LDA,&
AA(J_+1,J_+1),
2744 else if( INCX == 0 )
2748 else if( INCY == 0 )
2752 else if( LDA <
MAX2(1,M) )
2764 if( ((M == 0) || (N == 0)) || (ALPHA ==
ZERO) )
2777 JY = 1 - (N - 1)*INCY;
2781 for( J=1; J <= N; J++ )
2784 if( Y[JY-1] !=
ZERO )
2786 TEMP = ALPHA*Y[JY-1];
2787 for( I=1; I <=
M; I++ )
2790 AA(J_,I_) += X[I_]*TEMP;
2804 KX = 1 - (M - 1)*INCX;
2806 for( J=1; J <= N; J++ )
2809 if( Y[JY-1] !=
ZERO )
2811 TEMP = ALPHA*Y[JY-1];
2813 for( I=1; I <=
M; I++ )
2816 AA(J_,I_) += X[IX-1]*TEMP;
2997 NOTA =
LSAME(TRANSA,
'N');
2998 NOTB =
LSAME(TRANSB,
'N');
3021 if( ((!NOTA) && (!
LSAME(TRANSA,
'C'))) && (!
LSAME(TRANSA,
'T')) )
3026 ((!NOTB) && (!
LSAME(TRANSB,
'C'))) && (!
LSAME(TRANSB,
'T')) )
3046 else if( LDA <
MAX2(1,NROWA) )
3051 else if( LDB <
MAX2(1,NROWB) )
3056 else if( LDC <
MAX2(1,M) )
3069 if( ((M == 0) || (N == 0)) || (((ALPHA ==
ZERO) || (K == 0)) &&
3080 for( J=1; J <= N; J++ )
3083 for( I=1; I <=
M; I++ )
3093 for( J=1; J <= N; J++ )
3096 for( I=1; I <=
M; I++ )
3115 for( J=1; J <= N; J++ )
3120 for( I=1; I <=
M; I++ )
3127 else if( BETA !=
ONE )
3129 for( I=1; I <=
M; I++ )
3136 for( L=1; L <= K; L++ )
3141 TEMP = ALPHA*
BB(J_,L_);
3142 for( I=1; I <=
M; I++ )
3145 CC(J_,I_) += TEMP*
AA(L_,I_);
3155 for( J=1; J <= N; J++ )
3158 for( I=1; I <=
M; I++ )
3162 for( L=1; L <= K; L++ )
3165 TEMP +=
AA(I_,L_)*
BB(J_,L_);
3170 CC(J_,I_) = ALPHA*TEMP;
3174 CC(J_,I_) = ALPHA*TEMP + BETA*
CC(J_,I_);
3187 for( J=1; J <= N; J++ )
3192 for( I=1; I <=
M; I++ )
3199 else if( BETA !=
ONE )
3201 for( I=1; I <=
M; I++ )
3208 for( L=1; L <= K; L++ )
3213 TEMP = ALPHA*
BB(L_,J_);
3214 for( I=1; I <=
M; I++ )
3217 CC(J_,I_) += TEMP*
AA(L_,I_);
3228 for( J=1; J <= N; J++ )
3232 for( I=1; I <=
M; I++ )
3237 for( L=1; L <= K; L++ )
3240 TEMP +=
AA(I_,L_)*
BB(L_,J_);
3245 CC(J_,I_) = ALPHA*TEMP;
3249 CC(J_,I_) = ALPHA*TEMP + BETA*
CC(J_,I_);
3270 STATIC int32 DGTSV(int32 *n, int32 *nrhs,
double *dl,
3271 double *d__,
double *du,
double *b, int32 *ldb, int32
3275 int32 b_dim1, b_offset, i__1, i__2;
3282 #define b_ref(a_1,a_2) b[(a_2)*(b_dim1) + (a_1)]
3356 b_offset = 1 + b_dim1 * 1;
3363 }
else if(*nrhs < 0) {
3365 }
else if(*ldb < *n && *ldb < 1) {
3380 for(i__ = 1; i__ <= i__1; ++i__) {
3381 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3385 if(d__[i__] != 0.) {
3386 fact = dl[i__] / d__[i__];
3387 d__[i__ + 1] -= fact * du[i__];
3388 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3399 fact = d__[i__] / dl[i__];
3401 temp = d__[i__ + 1];
3402 d__[i__ + 1] = du[i__] - fact * temp;
3403 dl[i__] = du[i__ + 1];
3404 du[i__ + 1] = -fact * dl[i__];
3406 temp = b_ref(i__, 1);
3407 b_ref(i__, 1) = b_ref(i__ + 1, 1);
3408 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3414 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3415 if(d__[i__] != 0.) {
3416 fact = dl[i__] / d__[i__];
3417 d__[i__ + 1] -= fact * du[i__];
3418 b_ref(i__ + 1, 1) = b_ref(i__ + 1, 1) - fact * b_ref(i__,
3425 fact = d__[i__] / dl[i__];
3427 temp = d__[i__ + 1];
3428 d__[i__ + 1] = du[i__] - fact * temp;
3430 temp = b_ref(i__, 1);
3431 b_ref(i__, 1) = b_ref(i__ + 1, 1);
3432 b_ref(i__ + 1, 1) = temp - fact * b_ref(i__ + 1, 1);
3441 for(i__ = 1; i__ <= i__1; ++i__) {
3442 if(fabs(d__[i__]) >= fabs(dl[i__])) {
3446 if(d__[i__] != 0.) {
3447 fact = dl[i__] / d__[i__];
3448 d__[i__ + 1] -= fact * du[i__];
3450 for(j = 1; j <= i__2; ++j) {
3451 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3464 fact = d__[i__] / dl[i__];
3466 temp = d__[i__ + 1];
3467 d__[i__ + 1] = du[i__] - fact * temp;
3468 dl[i__] = du[i__ + 1];
3469 du[i__ + 1] = -fact * dl[i__];
3472 for(j = 1; j <= i__2; ++j) {
3473 temp = b_ref(i__, j);
3474 b_ref(i__, j) = b_ref(i__ + 1, j);
3475 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3483 if( fabs(d__[i__]) >= fabs(dl[i__]))
3487 fact = dl[i__] / d__[i__];
3488 d__[i__ + 1] -= fact * du[i__];
3490 for(j = 1; j <= i__1; ++j) {
3491 b_ref(i__ + 1, j) = b_ref(i__ + 1, j) - fact * b_ref(
3502 fact = d__[i__] / dl[i__];
3504 temp = d__[i__ + 1];
3505 d__[i__ + 1] = du[i__] - fact * temp;
3508 for(j = 1; j <= i__1; ++j) {
3509 temp = b_ref(i__, j);
3510 b_ref(i__, j) = b_ref(i__ + 1, j);
3511 b_ref(i__ + 1, j) = temp - fact * b_ref(i__ + 1, j);
3527 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3529 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n, j))
3532 for(i__ = *n - 2; i__ >= 1; --i__) {
3533 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j) - dl[
3534 i__] * b_ref(i__ + 2, j)) / d__[i__];
3543 for(j = 1; j <= i__1; ++j) {
3544 b_ref(*n, j) = b_ref(*n, j) / d__[*n];
3546 b_ref(*n - 1, j) = (b_ref(*n - 1, j) - du[*n - 1] * b_ref(*n,
3549 for(i__ = *n - 2; i__ >= 1; --i__) {
3550 b_ref(i__, j) = (b_ref(i__, j) - du[i__] * b_ref(i__ + 1, j)
3551 - dl[i__] * b_ref(i__ + 2, j)) / d__[i__];