📄 dgelsx.c
字号:
/* Scale A, B if max elements outside range [SMLNUM,BIGNUM] */
anrm = dlange_("M", m, n, &a[a_offset], lda, &work[1]);
iascl = 0;
if (anrm > 0. && anrm < smlnum) {
/* Scale matrix norm up to SMLNUM */
lstime_1.opcnt[gelsx - 1] += (doublereal) (*m * *n);
dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
info);
iascl = 1;
} else if (anrm > bignum) {
/* Scale matrix norm down to BIGNUM */
lstime_1.opcnt[gelsx - 1] += (doublereal) (*m * *n);
dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
info);
iascl = 2;
} else if (anrm == 0.) {
/* Matrix all zero. Return zero solution. */
i__1 = max(*m,*n);
dlaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb);
*rank = 0;
goto L100;
}
bnrm = dlange_("M", m, nrhs, &b[b_offset], ldb, &work[1]);
ibscl = 0;
if (bnrm > 0. && bnrm < smlnum) {
/* Scale matrix norm up to SMLNUM */
lstime_1.opcnt[gelsx - 1] += (doublereal) (*m * *nrhs);
dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb,
info);
ibscl = 1;
} else if (bnrm > bignum) {
/* Scale matrix norm down to BIGNUM */
lstime_1.opcnt[gelsx - 1] += (doublereal) (*m * *nrhs);
dlascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb,
info);
ibscl = 2;
}
/* Compute QR factorization with column pivoting of A:
A * P = Q * R */
lstime_1.opcnt[geqpf - 1] += dopla_("DGEQPF", m, n, &c__0, &c__0, &c__0);
tim1 = dsecnd_();
dgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], info);
tim2 = dsecnd_();
lstime_1.timng[geqpf - 1] += tim2 - tim1;
/* workspace 3*N. Details of Householder rotations stored
in WORK(1:MN).
Determine RANK using incremental condition estimation */
work[ismin] = 1.;
work[ismax] = 1.;
smax = (d__1 = a_ref(1, 1), abs(d__1));
smin = smax;
if ((d__1 = a_ref(1, 1), abs(d__1)) == 0.) {
*rank = 0;
i__1 = max(*m,*n);
dlaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb);
goto L100;
} else {
*rank = 1;
}
L10:
if (*rank < mn) {
i__ = *rank + 1;
latime_1.ops = 0.;
dlaic1_(&c__2, rank, &work[ismin], &smin, &a_ref(1, i__), &a_ref(i__,
i__), &sminpr, &s1, &c1);
dlaic1_(&c__1, rank, &work[ismax], &smax, &a_ref(1, i__), &a_ref(i__,
i__), &smaxpr, &s2, &c2);
lstime_1.opcnt[gelsx - 1] = lstime_1.opcnt[gelsx - 1] + latime_1.ops
+ 1.;
if (smaxpr * *rcond <= sminpr) {
lstime_1.opcnt[gelsx - 1] += (doublereal) (*rank << 1);
i__1 = *rank;
for (i__ = 1; i__ <= i__1; ++i__) {
work[ismin + i__ - 1] = s1 * work[ismin + i__ - 1];
work[ismax + i__ - 1] = s2 * work[ismax + i__ - 1];
/* L20: */
}
work[ismin + *rank] = c1;
work[ismax + *rank] = c2;
smin = sminpr;
smax = smaxpr;
++(*rank);
goto L10;
}
}
/* Logically partition R = [ R11 R12 ]
[ 0 R22 ]
where R11 = R(1:RANK,1:RANK)
[R11,R12] = [ T11, 0 ] * Y */
if (*rank < *n) {
lstime_1.opcnt[tzrqf - 1] += dopla_("DTZRQF", rank, n, &c__0, &c__0, &
c__0);
tim1 = dsecnd_();
dtzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info);
tim2 = dsecnd_();
lstime_1.timng[tzrqf - 1] += tim2 - tim1;
}
/* Details of Householder rotations stored in WORK(MN+1:2*MN)
B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
lstime_1.opcnt[orm2r - 1] += dopla_("DORMQR", m, nrhs, &mn, &c__0, &c__0);
tim1 = dsecnd_();
dorm2r_("Left", "Transpose", m, nrhs, &mn, &a[a_offset], lda, &work[1], &
b[b_offset], ldb, &work[(mn << 1) + 1], info);
tim2 = dsecnd_();
lstime_1.timng[orm2r - 1] += tim2 - tim1;
/* workspace NRHS
B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */
lstime_1.opcnt[trsm - 1] += dopbl3_("DTRSM ", rank, nrhs, &c__0);
tim1 = dsecnd_();
dtrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b49, &
a[a_offset], lda, &b[b_offset], ldb);
tim2 = dsecnd_();
lstime_1.timng[trsm - 1] += tim2 - tim1;
i__1 = *n;
for (i__ = *rank + 1; i__ <= i__1; ++i__) {
i__2 = *nrhs;
for (j = 1; j <= i__2; ++j) {
b_ref(i__, j) = 0.;
/* L30: */
}
/* L40: */
}
/* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */
if (*rank < *n) {
lstime_1.opcnt[latzm - 1] += (doublereal) (((*n - *rank) * *nrhs + *
nrhs + (*n - *rank) * *nrhs << 1) * *rank);
tim1 = dsecnd_();
i__1 = *rank;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *n - *rank + 1;
dlatzm_("Left", &i__2, nrhs, &a_ref(i__, *rank + 1), lda, &work[
mn + i__], &b_ref(i__, 1), &b_ref(*rank + 1, 1), ldb, &
work[(mn << 1) + 1]);
/* L50: */
}
tim2 = dsecnd_();
lstime_1.timng[latzm - 1] += tim2 - tim1;
}
/* workspace NRHS
B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */
i__1 = *nrhs;
for (j = 1; j <= i__1; ++j) {
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
work[(mn << 1) + i__] = 1.;
/* L60: */
}
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
if (work[(mn << 1) + i__] == 1.) {
if (jpvt[i__] != i__) {
k = i__;
t1 = b_ref(k, j);
t2 = b_ref(jpvt[k], j);
L70:
b_ref(jpvt[k], j) = t1;
work[(mn << 1) + k] = 0.;
t1 = t2;
k = jpvt[k];
t2 = b_ref(jpvt[k], j);
if (jpvt[k] != i__) {
goto L70;
}
b_ref(i__, j) = t1;
work[(mn << 1) + k] = 0.;
}
}
/* L80: */
}
/* L90: */
}
/* Undo scaling */
if (iascl == 1) {
lstime_1.opcnt[gelsx - 1] += (doublereal) (*n * *nrhs + *rank * *rank)
;
dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb,
info);
dlascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset],
lda, info);
} else if (iascl == 2) {
lstime_1.opcnt[gelsx - 1] += (doublereal) (*n * *nrhs + *rank * *rank)
;
dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb,
info);
dlascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset],
lda, info);
}
if (ibscl == 1) {
lstime_1.opcnt[gelsx - 1] += (doublereal) (*n * *nrhs);
dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb,
info);
} else if (ibscl == 2) {
lstime_1.opcnt[gelsx - 1] += (doublereal) (*n * *nrhs);
dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb,
info);
}
L100:
return 0;
/* End of DGELSX */
} /* dgelsx_ */
#undef b_ref
#undef a_ref
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -