37 #ifndef TEMPLATE_LAPACK_LATRS_HEADER
38 #define TEMPLATE_LAPACK_LATRS_HEADER
43 normin,
const integer *n,
const Treal *a,
const integer *lda, Treal *x,
44 Treal *scale, Treal *cnorm,
integer *info)
213 integer a_dim1, a_offset, i__1, i__2, i__3;
214 Treal d__1, d__2, d__3;
219 Treal tmax, tjjs, xmax, grow, sumj;
231 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
235 a_offset = 1 + a_dim1 * 1;
278 bignum = 1. / smlnum;
290 for (j = 1; j <= i__1; ++j) {
300 for (j = 1; j <= i__1; ++j) {
314 if (tmax <= bignum) {
317 tscal = 1. / (smlnum * tmax);
318 dscal_(n, &tscal, &cnorm[1], &c__1);
325 xmax = (d__1 = x[j],
absMACRO(d__1));
357 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
361 if (grow <= smlnum) {
369 d__1 = xbnd, d__2 =
minMACRO(1.,tjj) * grow;
371 if (tjj + cnorm[j] >= smlnum) {
375 grow *= tjj / (tjj + cnorm[j]);
392 d__1 = 1., d__2 = 1. /
maxMACRO(xbnd,smlnum);
396 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
400 if (grow <= smlnum) {
406 grow *= 1. / (cnorm[j] + 1.);
443 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
447 if (grow <= smlnum) {
455 d__1 = grow, d__2 = xbnd / xj;
474 d__1 = 1., d__2 = 1. /
maxMACRO(xbnd,smlnum);
478 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
482 if (grow <= smlnum) {
497 if (grow * tscal > smlnum) {
512 *scale = bignum / xmax;
513 dscal_(n, scale, &x[1], &c__1);
523 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
529 tjjs =
a_ref(j, j) * tscal;
542 if (xj > tjj * bignum) {
547 dscal_(n, &rec, &x[1], &c__1);
554 }
else if (tjj > 0.) {
558 if (xj > tjj * bignum) {
563 rec = tjj * bignum / xj;
571 dscal_(n, &rec, &x[1], &c__1);
583 for (i__ = 1; i__ <= i__3; ++i__) {
599 if (cnorm[j] > (bignum - xmax) * rec) {
604 dscal_(n, &rec, &x[1], &c__1);
607 }
else if (xj * cnorm[j] > bignum - xmax) {
611 dscal_(n, &c_b36, &x[1], &c__1);
622 d__1 = -x[j] * tscal;
627 xmax = (d__1 = x[i__],
absMACRO(d__1));
636 d__1 = -x[j] * tscal;
637 daxpy_(&i__3, &d__1, &
a_ref(j + 1, j), &c__1, &x[j +
641 xmax = (d__1 = x[i__],
absMACRO(d__1));
653 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
661 if (cnorm[j] > (bignum - xj) * rec) {
667 tjjs =
a_ref(j, j) * tscal;
677 d__1 = 1., d__2 = rec * tjj;
682 dscal_(n, &rec, &x[1], &c__1);
696 sumj =
ddot_(&i__3, &
a_ref(1, j), &c__1, &x[1], &c__1)
700 sumj =
ddot_(&i__3, &
a_ref(j + 1, j), &c__1, &x[j + 1]
709 for (i__ = 1; i__ <= i__3; ++i__) {
710 sumj +=
a_ref(i__, j) * uscal * x[i__];
715 for (i__ = j + 1; i__ <= i__3; ++i__) {
716 sumj +=
a_ref(i__, j) * uscal * x[i__];
722 if (uscal == tscal) {
730 tjjs =
a_ref(j, j) * tscal;
746 if (xj > tjj * bignum) {
751 dscal_(n, &rec, &x[1], &c__1);
757 }
else if (tjj > 0.) {
761 if (xj > tjj * bignum) {
765 rec = tjj * bignum / xj;
766 dscal_(n, &rec, &x[1], &c__1);
777 for (i__ = 1; i__ <= i__3; ++i__) {
792 x[j] = x[j] / tjjs - sumj;
795 d__2 = xmax, d__3 = (d__1 = x[j],
absMACRO(d__1));
807 dscal_(n, &d__1, &cnorm[1], &c__1);