35 #ifndef TEMPLATE_LAPACK_LATRS_HEADER
36 #define TEMPLATE_LAPACK_LATRS_HEADER
41 normin,
const integer *n,
const Treal *a,
const integer *lda, Treal *x,
42 Treal *scale, Treal *cnorm,
integer *info)
211 integer a_dim1, a_offset, i__1, i__2, i__3;
212 Treal d__1, d__2, d__3;
217 Treal tmax, tjjs, xmax, grow, sumj;
229 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
233 a_offset = 1 + a_dim1 * 1;
276 bignum = 1. / smlnum;
288 for (j = 1; j <= i__1; ++j) {
298 for (j = 1; j <= i__1; ++j) {
312 if (tmax <= bignum) {
315 tscal = 1. / (smlnum * tmax);
316 dscal_(n, &tscal, &cnorm[1], &c__1);
323 xmax = (d__1 = x[j],
absMACRO(d__1));
355 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
359 if (grow <= smlnum) {
367 d__1 = xbnd, d__2 =
minMACRO(1.,tjj) * grow;
369 if (tjj + cnorm[j] >= smlnum) {
373 grow *= tjj / (tjj + cnorm[j]);
390 d__1 = 1., d__2 = 1. /
maxMACRO(xbnd,smlnum);
394 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
398 if (grow <= smlnum) {
404 grow *= 1. / (cnorm[j] + 1.);
441 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
445 if (grow <= smlnum) {
453 d__1 = grow, d__2 = xbnd / xj;
472 d__1 = 1., d__2 = 1. /
maxMACRO(xbnd,smlnum);
476 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
480 if (grow <= smlnum) {
495 if (grow * tscal > smlnum) {
510 *scale = bignum / xmax;
511 dscal_(n, scale, &x[1], &c__1);
521 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
527 tjjs =
a_ref(j, j) * tscal;
540 if (xj > tjj * bignum) {
545 dscal_(n, &rec, &x[1], &c__1);
552 }
else if (tjj > 0.) {
556 if (xj > tjj * bignum) {
561 rec = tjj * bignum / xj;
569 dscal_(n, &rec, &x[1], &c__1);
581 for (i__ = 1; i__ <= i__3; ++i__) {
597 if (cnorm[j] > (bignum - xmax) * rec) {
602 dscal_(n, &rec, &x[1], &c__1);
605 }
else if (xj * cnorm[j] > bignum - xmax) {
609 dscal_(n, &c_b36, &x[1], &c__1);
620 d__1 = -x[j] * tscal;
625 xmax = (d__1 = x[i__],
absMACRO(d__1));
634 d__1 = -x[j] * tscal;
635 daxpy_(&i__3, &d__1, &
a_ref(j + 1, j), &c__1, &x[j +
639 xmax = (d__1 = x[i__],
absMACRO(d__1));
651 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
659 if (cnorm[j] > (bignum - xj) * rec) {
665 tjjs =
a_ref(j, j) * tscal;
675 d__1 = 1., d__2 = rec * tjj;
680 dscal_(n, &rec, &x[1], &c__1);
694 sumj =
ddot_(&i__3, &
a_ref(1, j), &c__1, &x[1], &c__1)
698 sumj =
ddot_(&i__3, &
a_ref(j + 1, j), &c__1, &x[j + 1]
707 for (i__ = 1; i__ <= i__3; ++i__) {
708 sumj +=
a_ref(i__, j) * uscal * x[i__];
713 for (i__ = j + 1; i__ <= i__3; ++i__) {
714 sumj +=
a_ref(i__, j) * uscal * x[i__];
720 if (uscal == tscal) {
728 tjjs =
a_ref(j, j) * tscal;
744 if (xj > tjj * bignum) {
749 dscal_(n, &rec, &x[1], &c__1);
755 }
else if (tjj > 0.) {
759 if (xj > tjj * bignum) {
763 rec = tjj * bignum / xj;
764 dscal_(n, &rec, &x[1], &c__1);
775 for (i__ = 1; i__ <= i__3; ++i__) {
790 x[j] = x[j] / tjjs - sumj;
793 d__2 = xmax, d__3 = (d__1 = x[j],
absMACRO(d__1));
805 dscal_(n, &d__1, &cnorm[1], &c__1);
#define absMACRO(x)
Definition: template_blas_common.h:45
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:40
int integer
Definition: template_blas_common.h:38
double ddot_(const int *n, const double *dx, const int *incx, const double *dy, const int *incy)
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_asum.h:40
#define minMACRO(a, b)
Definition: template_blas_common.h:44
void daxpy_(const int *n, const double *da, const double *dx, const int *incx, double *dy, const int *incy)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *normin, const integer *n, const Treal *a, const integer *lda, Treal *x, Treal *scale, Treal *cnorm, integer *info)
Definition: template_lapack_latrs.h:40
void dscal_(const int *n, const double *da, double *dx, const int *incx)
int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition: template_blas_trsv.h:40
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
bool logical
Definition: template_blas_common.h:39
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44