📄 dsptrs.c
字号:
/*< KC = 1 >*/
kc = 1;
/*< 40 CONTINUE >*/
L40:
/* If K > N, exit from loop. */
/*< >*/
if (k > *n) {
goto L50;
}
/*< IF( IPIV( K ).GT.0 ) THEN >*/
if (ipiv[k] > 0) {
/* 1 x 1 diagonal block */
/* Multiply by inv(U'(K)), where U(K) is the transformation */
/* stored in column K of A. */
/*< >*/
i__1 = k - 1;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc]
, &c__1, &c_b19, &b[k + b_dim1], ldb, (ftnlen)9);
/* Interchange rows K and IPIV(K). */
/*< KP = IPIV( K ) >*/
kp = ipiv[k];
/*< >*/
if (kp != k) {
dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/*< KC = KC + K >*/
kc += k;
/*< K = K + 1 >*/
++k;
/*< ELSE >*/
} else {
/* 2 x 2 diagonal block */
/* Multiply by inv(U'(K+1)), where U(K+1) is the transformation */
/* stored in columns K and K+1 of A. */
/*< >*/
i__1 = k - 1;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc]
, &c__1, &c_b19, &b[k + b_dim1], ldb, (ftnlen)9);
/*< >*/
i__1 = k - 1;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &ap[kc
+ k], &c__1, &c_b19, &b[k + 1 + b_dim1], ldb, (ftnlen)9);
/* Interchange rows K and -IPIV(K). */
/*< KP = -IPIV( K ) >*/
kp = -ipiv[k];
/*< >*/
if (kp != k) {
dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/*< KC = KC + 2*K + 1 >*/
kc += (k << 1) + 1;
/*< K = K + 2 >*/
k += 2;
/*< END IF >*/
}
/*< GO TO 40 >*/
goto L40;
/*< 50 CONTINUE >*/
L50:
/*< ELSE >*/
;
} else {
/* Solve A*X = B, where A = L*D*L'. */
/* First solve L*D*X = B, overwriting B with X. */
/* K is the main loop index, increasing from 1 to N in steps of */
/* 1 or 2, depending on the size of the diagonal blocks. */
/*< K = 1 >*/
k = 1;
/*< KC = 1 >*/
kc = 1;
/*< 60 CONTINUE >*/
L60:
/* If K > N, exit from loop. */
/*< >*/
if (k > *n) {
goto L80;
}
/*< IF( IPIV( K ).GT.0 ) THEN >*/
if (ipiv[k] > 0) {
/* 1 x 1 diagonal block */
/* Interchange rows K and IPIV(K). */
/*< KP = IPIV( K ) >*/
kp = ipiv[k];
/*< >*/
if (kp != k) {
dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/* Multiply by inv(L(K)), where L(K) is the transformation */
/* stored in column K of A. */
/*< >*/
if (k < *n) {
i__1 = *n - k;
dger_(&i__1, nrhs, &c_b7, &ap[kc + 1], &c__1, &b[k + b_dim1],
ldb, &b[k + 1 + b_dim1], ldb);
}
/* Multiply by the inverse of the diagonal block. */
/*< CALL DSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB ) >*/
d__1 = 1. / ap[kc];
dscal_(nrhs, &d__1, &b[k + b_dim1], ldb);
/*< KC = KC + N - K + 1 >*/
kc += *n - k + 1;
/*< K = K + 1 >*/
++k;
/*< ELSE >*/
} else {
/* 2 x 2 diagonal block */
/* Interchange rows K+1 and -IPIV(K). */
/*< KP = -IPIV( K ) >*/
kp = -ipiv[k];
/*< >*/
if (kp != k + 1) {
dswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/* Multiply by inv(L(K)), where L(K) is the transformation */
/* stored in columns K and K+1 of A. */
/*< IF( K.LT.N-1 ) THEN >*/
if (k < *n - 1) {
/*< >*/
i__1 = *n - k - 1;
dger_(&i__1, nrhs, &c_b7, &ap[kc + 2], &c__1, &b[k + b_dim1],
ldb, &b[k + 2 + b_dim1], ldb);
/*< >*/
i__1 = *n - k - 1;
dger_(&i__1, nrhs, &c_b7, &ap[kc + *n - k + 2], &c__1, &b[k +
1 + b_dim1], ldb, &b[k + 2 + b_dim1], ldb);
/*< END IF >*/
}
/* Multiply by the inverse of the diagonal block. */
/*< AKM1K = AP( KC+1 ) >*/
akm1k = ap[kc + 1];
/*< AKM1 = AP( KC ) / AKM1K >*/
akm1 = ap[kc] / akm1k;
/*< AK = AP( KC+N-K+1 ) / AKM1K >*/
ak = ap[kc + *n - k + 1] / akm1k;
/*< DENOM = AKM1*AK - ONE >*/
denom = akm1 * ak - 1.;
/*< DO 70 J = 1, NRHS >*/
i__1 = *nrhs;
for (j = 1; j <= i__1; ++j) {
/*< BKM1 = B( K, J ) / AKM1K >*/
bkm1 = b[k + j * b_dim1] / akm1k;
/*< BK = B( K+1, J ) / AKM1K >*/
bk = b[k + 1 + j * b_dim1] / akm1k;
/*< B( K, J ) = ( AK*BKM1-BK ) / DENOM >*/
b[k + j * b_dim1] = (ak * bkm1 - bk) / denom;
/*< B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM >*/
b[k + 1 + j * b_dim1] = (akm1 * bk - bkm1) / denom;
/*< 70 CONTINUE >*/
/* L70: */
}
/*< KC = KC + 2*( N-K ) + 1 >*/
kc += ((*n - k) << 1) + 1;
/*< K = K + 2 >*/
k += 2;
/*< END IF >*/
}
/*< GO TO 60 >*/
goto L60;
/*< 80 CONTINUE >*/
L80:
/* Next solve L'*X = B, overwriting B with X. */
/* K is the main loop index, decreasing from N to 1 in steps of */
/* 1 or 2, depending on the size of the diagonal blocks. */
/*< K = N >*/
k = *n;
/*< KC = N*( N+1 ) / 2 + 1 >*/
kc = *n * (*n + 1) / 2 + 1;
/*< 90 CONTINUE >*/
L90:
/* If K < 1, exit from loop. */
/*< >*/
if (k < 1) {
goto L100;
}
/*< KC = KC - ( N-K+1 ) >*/
kc -= *n - k + 1;
/*< IF( IPIV( K ).GT.0 ) THEN >*/
if (ipiv[k] > 0) {
/* 1 x 1 diagonal block */
/* Multiply by inv(L'(K)), where L(K) is the transformation */
/* stored in column K of A. */
/*< >*/
if (k < *n) {
i__1 = *n - k;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1],
ldb, &ap[kc + 1], &c__1, &c_b19, &b[k + b_dim1], ldb,
(ftnlen)9);
}
/* Interchange rows K and IPIV(K). */
/*< KP = IPIV( K ) >*/
kp = ipiv[k];
/*< >*/
if (kp != k) {
dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/*< K = K - 1 >*/
--k;
/*< ELSE >*/
} else {
/* 2 x 2 diagonal block */
/* Multiply by inv(L'(K-1)), where L(K-1) is the transformation */
/* stored in columns K-1 and K of A. */
/*< IF( K.LT.N ) THEN >*/
if (k < *n) {
/*< >*/
i__1 = *n - k;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1],
ldb, &ap[kc + 1], &c__1, &c_b19, &b[k + b_dim1], ldb,
(ftnlen)9);
/*< >*/
i__1 = *n - k;
dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[k + 1 + b_dim1],
ldb, &ap[kc - (*n - k)], &c__1, &c_b19, &b[k - 1 +
b_dim1], ldb, (ftnlen)9);
/*< END IF >*/
}
/* Interchange rows K and -IPIV(K). */
/*< KP = -IPIV( K ) >*/
kp = -ipiv[k];
/*< >*/
if (kp != k) {
dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
}
/*< KC = KC - ( N-K+2 ) >*/
kc -= *n - k + 2;
/*< K = K - 2 >*/
k += -2;
/*< END IF >*/
}
/*< GO TO 90 >*/
goto L90;
/*< 100 CONTINUE >*/
L100:
/*< END IF >*/
;
}
/*< RETURN >*/
return 0;
/* End of DSPTRS */
/*< END >*/
} /* dsptrs_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -