template_lapack_latrs.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027  
00028  /* This file belongs to the template_lapack part of the Ergo source 
00029   * code. The source files in the template_lapack directory are modified
00030   * versions of files originally distributed as CLAPACK, see the
00031   * Copyright/license notice in the file template_lapack/COPYING.
00032   */
00033  
00034 
00035 #ifndef TEMPLATE_LAPACK_LATRS_HEADER
00036 #define TEMPLATE_LAPACK_LATRS_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *
00041         normin, const integer *n, const Treal *a, const integer *lda, Treal *x, 
00042         Treal *scale, Treal *cnorm, integer *info)
00043 {
00044 /*  -- LAPACK auxiliary routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        June 30, 1992   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DLATRS solves one of the triangular systems   
00054 
00055        A *x = s*b  or  A'*x = s*b   
00056 
00057     with scaling to prevent overflow.  Here A is an upper or lower   
00058     triangular matrix, A' denotes the transpose of A, x and b are   
00059     n-element vectors, and s is a scaling factor, usually less than   
00060     or equal to 1, chosen so that the components of x will be less than   
00061     the overflow threshold.  If the unscaled problem will not cause   
00062     overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A   
00063     is singular (A(j,j) = 0 for some j), then s is set to 0 and a   
00064     non-trivial solution to A*x = 0 is returned.   
00065 
00066     Arguments   
00067     =========   
00068 
00069     UPLO    (input) CHARACTER*1   
00070             Specifies whether the matrix A is upper or lower triangular.   
00071             = 'U':  Upper triangular   
00072             = 'L':  Lower triangular   
00073 
00074     TRANS   (input) CHARACTER*1   
00075             Specifies the operation applied to A.   
00076             = 'N':  Solve A * x = s*b  (No transpose)   
00077             = 'T':  Solve A'* x = s*b  (Transpose)   
00078             = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)   
00079 
00080     DIAG    (input) CHARACTER*1   
00081             Specifies whether or not the matrix A is unit triangular.   
00082             = 'N':  Non-unit triangular   
00083             = 'U':  Unit triangular   
00084 
00085     NORMIN  (input) CHARACTER*1   
00086             Specifies whether CNORM has been set or not.   
00087             = 'Y':  CNORM contains the column norms on entry   
00088             = 'N':  CNORM is not set on entry.  On exit, the norms will   
00089                     be computed and stored in CNORM.   
00090 
00091     N       (input) INTEGER   
00092             The order of the matrix A.  N >= 0.   
00093 
00094     A       (input) DOUBLE PRECISION array, dimension (LDA,N)   
00095             The triangular matrix A.  If UPLO = 'U', the leading n by n   
00096             upper triangular part of the array A contains the upper   
00097             triangular matrix, and the strictly lower triangular part of   
00098             A is not referenced.  If UPLO = 'L', the leading n by n lower   
00099             triangular part of the array A contains the lower triangular   
00100             matrix, and the strictly upper triangular part of A is not   
00101             referenced.  If DIAG = 'U', the diagonal elements of A are   
00102             also not referenced and are assumed to be 1.   
00103 
00104     LDA     (input) INTEGER   
00105             The leading dimension of the array A.  LDA >= max (1,N).   
00106 
00107     X       (input/output) DOUBLE PRECISION array, dimension (N)   
00108             On entry, the right hand side b of the triangular system.   
00109             On exit, X is overwritten by the solution vector x.   
00110 
00111     SCALE   (output) DOUBLE PRECISION   
00112             The scaling factor s for the triangular system   
00113                A * x = s*b  or  A'* x = s*b.   
00114             If SCALE = 0, the matrix A is singular or badly scaled, and   
00115             the vector x is an exact or approximate solution to A*x = 0.   
00116 
00117     CNORM   (input or output) DOUBLE PRECISION array, dimension (N)   
00118 
00119             If NORMIN = 'Y', CNORM is an input argument and CNORM(j)   
00120             contains the norm of the off-diagonal part of the j-th column   
00121             of A.  If TRANS = 'N', CNORM(j) must be greater than or equal   
00122             to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)   
00123             must be greater than or equal to the 1-norm.   
00124 
00125             If NORMIN = 'N', CNORM is an output argument and CNORM(j)   
00126             returns the 1-norm of the offdiagonal part of the j-th column   
00127             of A.   
00128 
00129     INFO    (output) INTEGER   
00130             = 0:  successful exit   
00131             < 0:  if INFO = -k, the k-th argument had an illegal value   
00132 
00133     Further Details   
00134     ======= =======   
00135 
00136     A rough bound on x is computed; if that is less than overflow, DTRSV   
00137     is called, otherwise, specific code is used which checks for possible   
00138     overflow or divide-by-zero at every operation.   
00139 
00140     A columnwise scheme is used for solving A*x = b.  The basic algorithm   
00141     if A is lower triangular is   
00142 
00143          x[1:n] := b[1:n]   
00144          for j = 1, ..., n   
00145               x(j) := x(j) / A(j,j)   
00146               x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]   
00147          end   
00148 
00149     Define bounds on the components of x after j iterations of the loop:   
00150        M(j) = bound on x[1:j]   
00151        G(j) = bound on x[j+1:n]   
00152     Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.   
00153 
00154     Then for iteration j+1 we have   
00155        M(j+1) <= G(j) / | A(j+1,j+1) |   
00156        G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |   
00157               <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )   
00158 
00159     where CNORM(j+1) is greater than or equal to the infinity-norm of   
00160     column j+1 of A, not counting the diagonal.  Hence   
00161 
00162        G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )   
00163                     1<=i<=j   
00164     and   
00165 
00166        |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )   
00167                                      1<=i< j   
00168 
00169     Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the   
00170     reciprocal of the largest M(j), j=1,..,n, is larger than   
00171     max(underflow, 1/overflow).   
00172 
00173     The bound on x(j) is also used to determine when a step in the   
00174     columnwise method can be performed without fear of overflow.  If   
00175     the computed bound is greater than a large constant, x is scaled to   
00176     prevent overflow, but if the bound overflows, x is set to 0, x(j) to   
00177     1, and scale to 0, and a non-trivial solution to A*x = 0 is found.   
00178 
00179     Similarly, a row-wise scheme is used to solve A'*x = b.  The basic   
00180     algorithm for A upper triangular is   
00181 
00182          for j = 1, ..., n   
00183               x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)   
00184          end   
00185 
00186     We simultaneously compute two bounds   
00187          G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j   
00188          M(j) = bound on x(i), 1<=i<=j   
00189 
00190     The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we   
00191     add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.   
00192     Then the bound on x(j) is   
00193 
00194          M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |   
00195 
00196               <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )   
00197                         1<=i<=j   
00198 
00199     and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater   
00200     than max(underflow, 1/overflow).   
00201 
00202     =====================================================================   
00203 
00204 
00205        Parameter adjustments */
00206     /* Table of constant values */
00207      integer c__1 = 1;
00208      Treal c_b36 = .5;
00209     
00210     /* System generated locals */
00211     integer a_dim1, a_offset, i__1, i__2, i__3;
00212     Treal d__1, d__2, d__3;
00213     /* Local variables */
00214      integer jinc;
00215      Treal xbnd;
00216      integer imax;
00217      Treal tmax, tjjs, xmax, grow, sumj;
00218      integer i__, j;
00219      Treal tscal, uscal;
00220      integer jlast;
00221      logical upper;
00222      Treal xj;
00223      Treal bignum;
00224      logical notran;
00225      integer jfirst;
00226      Treal smlnum;
00227      logical nounit;
00228      Treal rec, tjj;
00229 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00230 
00231 
00232     a_dim1 = *lda;
00233     a_offset = 1 + a_dim1 * 1;
00234     a -= a_offset;
00235     --x;
00236     --cnorm;
00237 
00238     /* Function Body */
00239     *info = 0;
00240     upper = template_blas_lsame(uplo, "U");
00241     notran = template_blas_lsame(trans, "N");
00242     nounit = template_blas_lsame(diag, "N");
00243 
00244 /*     Test the input parameters. */
00245 
00246     if (! upper && ! template_blas_lsame(uplo, "L")) {
00247         *info = -1;
00248     } else if (! notran && ! template_blas_lsame(trans, "T") && ! 
00249             template_blas_lsame(trans, "C")) {
00250         *info = -2;
00251     } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00252         *info = -3;
00253     } else if (! template_blas_lsame(normin, "Y") && ! template_blas_lsame(normin,
00254              "N")) {
00255         *info = -4;
00256     } else if (*n < 0) {
00257         *info = -5;
00258     } else if (*lda < maxMACRO(1,*n)) {
00259         *info = -7;
00260     }
00261     if (*info != 0) {
00262         i__1 = -(*info);
00263         template_blas_erbla("LATRS ", &i__1);
00264         return 0;
00265     }
00266 
00267 /*     Quick return if possible */
00268 
00269     if (*n == 0) {
00270         return 0;
00271     }
00272 
00273 /*     Determine machine dependent parameters to control overflow. */
00274 
00275     smlnum = template_lapack_lamch("Safe minimum", (Treal)0) / template_lapack_lamch("Precision", (Treal)0);
00276     bignum = 1. / smlnum;
00277     *scale = 1.;
00278 
00279     if (template_blas_lsame(normin, "N")) {
00280 
00281 /*        Compute the 1-norm of each column, not including the diagonal. */
00282 
00283         if (upper) {
00284 
00285 /*           A is upper triangular. */
00286 
00287             i__1 = *n;
00288             for (j = 1; j <= i__1; ++j) {
00289                 i__2 = j - 1;
00290                 cnorm[j] = template_blas_asum(&i__2, &a_ref(1, j), &c__1);
00291 /* L10: */
00292             }
00293         } else {
00294 
00295 /*           A is lower triangular. */
00296 
00297             i__1 = *n - 1;
00298             for (j = 1; j <= i__1; ++j) {
00299                 i__2 = *n - j;
00300                 cnorm[j] = template_blas_asum(&i__2, &a_ref(j + 1, j), &c__1);
00301 /* L20: */
00302             }
00303             cnorm[*n] = 0.;
00304         }
00305     }
00306 
00307 /*     Scale the column norms by TSCAL if the maximum element in CNORM is   
00308        greater than BIGNUM. */
00309 
00310     imax = template_blas_idamax(n, &cnorm[1], &c__1);
00311     tmax = cnorm[imax];
00312     if (tmax <= bignum) {
00313         tscal = 1.;
00314     } else {
00315         tscal = 1. / (smlnum * tmax);
00316         dscal_(n, &tscal, &cnorm[1], &c__1);
00317     }
00318 
00319 /*     Compute a bound on the computed solution vector to see if the   
00320        Level 2 BLAS routine DTRSV can be used. */
00321 
00322     j = template_blas_idamax(n, &x[1], &c__1);
00323     xmax = (d__1 = x[j], absMACRO(d__1));
00324     xbnd = xmax;
00325     if (notran) {
00326 
00327 /*        Compute the growth in A * x = b. */
00328 
00329         if (upper) {
00330             jfirst = *n;
00331             jlast = 1;
00332             jinc = -1;
00333         } else {
00334             jfirst = 1;
00335             jlast = *n;
00336             jinc = 1;
00337         }
00338 
00339         if (tscal != 1.) {
00340             grow = 0.;
00341             goto L50;
00342         }
00343 
00344         if (nounit) {
00345 
00346 /*           A is non-unit triangular.   
00347 
00348              Compute GROW = 1/G(j) and XBND = 1/M(j).   
00349              Initially, G(0) = max{x(i), i=1,...,n}. */
00350 
00351             grow = 1. / maxMACRO(xbnd,smlnum);
00352             xbnd = grow;
00353             i__1 = jlast;
00354             i__2 = jinc;
00355             for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00356 
00357 /*              Exit the loop if the growth factor is too small. */
00358 
00359                 if (grow <= smlnum) {
00360                     goto L50;
00361                 }
00362 
00363 /*              M(j) = G(j-1) / absMACRO(A(j,j)) */
00364 
00365                 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00366 /* Computing MIN */
00367                 d__1 = xbnd, d__2 = minMACRO(1.,tjj) * grow;
00368                 xbnd = minMACRO(d__1,d__2);
00369                 if (tjj + cnorm[j] >= smlnum) {
00370 
00371 /*                 G(j) = G(j-1)*( 1 + CNORM(j) / absMACRO(A(j,j)) ) */
00372 
00373                     grow *= tjj / (tjj + cnorm[j]);
00374                 } else {
00375 
00376 /*                 G(j) could overflow, set GROW to 0. */
00377 
00378                     grow = 0.;
00379                 }
00380 /* L30: */
00381             }
00382             grow = xbnd;
00383         } else {
00384 
00385 /*           A is unit triangular.   
00386 
00387              Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.   
00388 
00389    Computing MIN */
00390             d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00391             grow = minMACRO(d__1,d__2);
00392             i__2 = jlast;
00393             i__1 = jinc;
00394             for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00395 
00396 /*              Exit the loop if the growth factor is too small. */
00397 
00398                 if (grow <= smlnum) {
00399                     goto L50;
00400                 }
00401 
00402 /*              G(j) = G(j-1)*( 1 + CNORM(j) ) */
00403 
00404                 grow *= 1. / (cnorm[j] + 1.);
00405 /* L40: */
00406             }
00407         }
00408 L50:
00409 
00410         ;
00411     } else {
00412 
00413 /*        Compute the growth in A' * x = b. */
00414 
00415         if (upper) {
00416             jfirst = 1;
00417             jlast = *n;
00418             jinc = 1;
00419         } else {
00420             jfirst = *n;
00421             jlast = 1;
00422             jinc = -1;
00423         }
00424 
00425         if (tscal != 1.) {
00426             grow = 0.;
00427             goto L80;
00428         }
00429 
00430         if (nounit) {
00431 
00432 /*           A is non-unit triangular.   
00433 
00434              Compute GROW = 1/G(j) and XBND = 1/M(j).   
00435              Initially, M(0) = max{x(i), i=1,...,n}. */
00436 
00437             grow = 1. / maxMACRO(xbnd,smlnum);
00438             xbnd = grow;
00439             i__1 = jlast;
00440             i__2 = jinc;
00441             for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00442 
00443 /*              Exit the loop if the growth factor is too small. */
00444 
00445                 if (grow <= smlnum) {
00446                     goto L80;
00447                 }
00448 
00449 /*              G(j) = maxMACRO( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */
00450 
00451                 xj = cnorm[j] + 1.;
00452 /* Computing MIN */
00453                 d__1 = grow, d__2 = xbnd / xj;
00454                 grow = minMACRO(d__1,d__2);
00455 
00456 /*              M(j) = M(j-1)*( 1 + CNORM(j) ) / absMACRO(A(j,j)) */
00457 
00458                 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00459                 if (xj > tjj) {
00460                     xbnd *= tjj / xj;
00461                 }
00462 /* L60: */
00463             }
00464             grow = minMACRO(grow,xbnd);
00465         } else {
00466 
00467 /*           A is unit triangular.   
00468 
00469              Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.   
00470 
00471    Computing MIN */
00472             d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00473             grow = minMACRO(d__1,d__2);
00474             i__2 = jlast;
00475             i__1 = jinc;
00476             for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00477 
00478 /*              Exit the loop if the growth factor is too small. */
00479 
00480                 if (grow <= smlnum) {
00481                     goto L80;
00482                 }
00483 
00484 /*              G(j) = ( 1 + CNORM(j) )*G(j-1) */
00485 
00486                 xj = cnorm[j] + 1.;
00487                 grow /= xj;
00488 /* L70: */
00489             }
00490         }
00491 L80:
00492         ;
00493     }
00494 
00495     if (grow * tscal > smlnum) {
00496 
00497 /*        Use the Level 2 BLAS solve if the reciprocal of the bound on   
00498           elements of X is not too small. */
00499 
00500         template_blas_trsv(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1);
00501     } else {
00502 
00503 /*        Use a Level 1 BLAS solve, scaling intermediate results. */
00504 
00505         if (xmax > bignum) {
00506 
00507 /*           Scale X so that its components are less than or equal to   
00508              BIGNUM in absolute value. */
00509 
00510             *scale = bignum / xmax;
00511             dscal_(n, scale, &x[1], &c__1);
00512             xmax = bignum;
00513         }
00514 
00515         if (notran) {
00516 
00517 /*           Solve A * x = b */
00518 
00519             i__1 = jlast;
00520             i__2 = jinc;
00521             for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00522 
00523 /*              Compute x(j) = b(j) / A(j,j), scaling x if necessary. */
00524 
00525                 xj = (d__1 = x[j], absMACRO(d__1));
00526                 if (nounit) {
00527                     tjjs = a_ref(j, j) * tscal;
00528                 } else {
00529                     tjjs = tscal;
00530                     if (tscal == 1.) {
00531                         goto L100;
00532                     }
00533                 }
00534                 tjj = absMACRO(tjjs);
00535                 if (tjj > smlnum) {
00536 
00537 /*                    absMACRO(A(j,j)) > SMLNUM: */
00538 
00539                     if (tjj < 1.) {
00540                         if (xj > tjj * bignum) {
00541 
00542 /*                          Scale x by 1/b(j). */
00543 
00544                             rec = 1. / xj;
00545                             dscal_(n, &rec, &x[1], &c__1);
00546                             *scale *= rec;
00547                             xmax *= rec;
00548                         }
00549                     }
00550                     x[j] /= tjjs;
00551                     xj = (d__1 = x[j], absMACRO(d__1));
00552                 } else if (tjj > 0.) {
00553 
00554 /*                    0 < absMACRO(A(j,j)) <= SMLNUM: */
00555 
00556                     if (xj > tjj * bignum) {
00557 
00558 /*                       Scale x by (1/absMACRO(x(j)))*absMACRO(A(j,j))*BIGNUM   
00559                          to avoid overflow when dividing by A(j,j). */
00560 
00561                         rec = tjj * bignum / xj;
00562                         if (cnorm[j] > 1.) {
00563 
00564 /*                          Scale by 1/CNORM(j) to avoid overflow when   
00565                             multiplying x(j) times column j. */
00566 
00567                             rec /= cnorm[j];
00568                         }
00569                         dscal_(n, &rec, &x[1], &c__1);
00570                         *scale *= rec;
00571                         xmax *= rec;
00572                     }
00573                     x[j] /= tjjs;
00574                     xj = (d__1 = x[j], absMACRO(d__1));
00575                 } else {
00576 
00577 /*                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and   
00578                       scale = 0, and compute a solution to A*x = 0. */
00579 
00580                     i__3 = *n;
00581                     for (i__ = 1; i__ <= i__3; ++i__) {
00582                         x[i__] = 0.;
00583 /* L90: */
00584                     }
00585                     x[j] = 1.;
00586                     xj = 1.;
00587                     *scale = 0.;
00588                     xmax = 0.;
00589                 }
00590 L100:
00591 
00592 /*              Scale x if necessary to avoid overflow when adding a   
00593                 multiple of column j of A. */
00594 
00595                 if (xj > 1.) {
00596                     rec = 1. / xj;
00597                     if (cnorm[j] > (bignum - xmax) * rec) {
00598 
00599 /*                    Scale x by 1/(2*absMACRO(x(j))). */
00600 
00601                         rec *= .5;
00602                         dscal_(n, &rec, &x[1], &c__1);
00603                         *scale *= rec;
00604                     }
00605                 } else if (xj * cnorm[j] > bignum - xmax) {
00606 
00607 /*                 Scale x by 1/2. */
00608 
00609                     dscal_(n, &c_b36, &x[1], &c__1);
00610                     *scale *= .5;
00611                 }
00612 
00613                 if (upper) {
00614                     if (j > 1) {
00615 
00616 /*                    Compute the update   
00617                          x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) */
00618 
00619                         i__3 = j - 1;
00620                         d__1 = -x[j] * tscal;
00621                         daxpy_(&i__3, &d__1, &a_ref(1, j), &c__1, &x[1], &
00622                                 c__1);
00623                         i__3 = j - 1;
00624                         i__ = template_blas_idamax(&i__3, &x[1], &c__1);
00625                         xmax = (d__1 = x[i__], absMACRO(d__1));
00626                     }
00627                 } else {
00628                     if (j < *n) {
00629 
00630 /*                    Compute the update   
00631                          x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) */
00632 
00633                         i__3 = *n - j;
00634                         d__1 = -x[j] * tscal;
00635                         daxpy_(&i__3, &d__1, &a_ref(j + 1, j), &c__1, &x[j + 
00636                                 1], &c__1);
00637                         i__3 = *n - j;
00638                         i__ = j + template_blas_idamax(&i__3, &x[j + 1], &c__1);
00639                         xmax = (d__1 = x[i__], absMACRO(d__1));
00640                     }
00641                 }
00642 /* L110: */
00643             }
00644 
00645         } else {
00646 
00647 /*           Solve A' * x = b */
00648 
00649             i__2 = jlast;
00650             i__1 = jinc;
00651             for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00652 
00653 /*              Compute x(j) = b(j) - sum A(k,j)*x(k).   
00654                                       k<>j */
00655 
00656                 xj = (d__1 = x[j], absMACRO(d__1));
00657                 uscal = tscal;
00658                 rec = 1. / maxMACRO(xmax,1.);
00659                 if (cnorm[j] > (bignum - xj) * rec) {
00660 
00661 /*                 If x(j) could overflow, scale x by 1/(2*XMAX). */
00662 
00663                     rec *= .5;
00664                     if (nounit) {
00665                         tjjs = a_ref(j, j) * tscal;
00666                     } else {
00667                         tjjs = tscal;
00668                     }
00669                     tjj = absMACRO(tjjs);
00670                     if (tjj > 1.) {
00671 
00672 /*                       Divide by A(j,j) when scaling x if A(j,j) > 1.   
00673 
00674    Computing MIN */
00675                         d__1 = 1., d__2 = rec * tjj;
00676                         rec = minMACRO(d__1,d__2);
00677                         uscal /= tjjs;
00678                     }
00679                     if (rec < 1.) {
00680                         dscal_(n, &rec, &x[1], &c__1);
00681                         *scale *= rec;
00682                         xmax *= rec;
00683                     }
00684                 }
00685 
00686                 sumj = 0.;
00687                 if (uscal == 1.) {
00688 
00689 /*                 If the scaling needed for A in the dot product is 1,   
00690                    call DDOT to perform the dot product. */
00691 
00692                     if (upper) {
00693                         i__3 = j - 1;
00694                         sumj = ddot_(&i__3, &a_ref(1, j), &c__1, &x[1], &c__1)
00695                                 ;
00696                     } else if (j < *n) {
00697                         i__3 = *n - j;
00698                         sumj = ddot_(&i__3, &a_ref(j + 1, j), &c__1, &x[j + 1]
00699                                 , &c__1);
00700                     }
00701                 } else {
00702 
00703 /*                 Otherwise, use in-line code for the dot product. */
00704 
00705                     if (upper) {
00706                         i__3 = j - 1;
00707                         for (i__ = 1; i__ <= i__3; ++i__) {
00708                             sumj += a_ref(i__, j) * uscal * x[i__];
00709 /* L120: */
00710                         }
00711                     } else if (j < *n) {
00712                         i__3 = *n;
00713                         for (i__ = j + 1; i__ <= i__3; ++i__) {
00714                             sumj += a_ref(i__, j) * uscal * x[i__];
00715 /* L130: */
00716                         }
00717                     }
00718                 }
00719 
00720                 if (uscal == tscal) {
00721 
00722 /*                 Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)   
00723                    was not used to scale the dotproduct. */
00724 
00725                     x[j] -= sumj;
00726                     xj = (d__1 = x[j], absMACRO(d__1));
00727                     if (nounit) {
00728                         tjjs = a_ref(j, j) * tscal;
00729                     } else {
00730                         tjjs = tscal;
00731                         if (tscal == 1.) {
00732                             goto L150;
00733                         }
00734                     }
00735 
00736 /*                    Compute x(j) = x(j) / A(j,j), scaling if necessary. */
00737 
00738                     tjj = absMACRO(tjjs);
00739                     if (tjj > smlnum) {
00740 
00741 /*                       absMACRO(A(j,j)) > SMLNUM: */
00742 
00743                         if (tjj < 1.) {
00744                             if (xj > tjj * bignum) {
00745 
00746 /*                             Scale X by 1/absMACRO(x(j)). */
00747 
00748                                 rec = 1. / xj;
00749                                 dscal_(n, &rec, &x[1], &c__1);
00750                                 *scale *= rec;
00751                                 xmax *= rec;
00752                             }
00753                         }
00754                         x[j] /= tjjs;
00755                     } else if (tjj > 0.) {
00756 
00757 /*                       0 < absMACRO(A(j,j)) <= SMLNUM: */
00758 
00759                         if (xj > tjj * bignum) {
00760 
00761 /*                          Scale x by (1/absMACRO(x(j)))*absMACRO(A(j,j))*BIGNUM. */
00762 
00763                             rec = tjj * bignum / xj;
00764                             dscal_(n, &rec, &x[1], &c__1);
00765                             *scale *= rec;
00766                             xmax *= rec;
00767                         }
00768                         x[j] /= tjjs;
00769                     } else {
00770 
00771 /*                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and   
00772                          scale = 0, and compute a solution to A'*x = 0. */
00773 
00774                         i__3 = *n;
00775                         for (i__ = 1; i__ <= i__3; ++i__) {
00776                             x[i__] = 0.;
00777 /* L140: */
00778                         }
00779                         x[j] = 1.;
00780                         *scale = 0.;
00781                         xmax = 0.;
00782                     }
00783 L150:
00784                     ;
00785                 } else {
00786 
00787 /*                 Compute x(j) := x(j) / A(j,j)  - sumj if the dot   
00788                    product has already been divided by 1/A(j,j). */
00789 
00790                     x[j] = x[j] / tjjs - sumj;
00791                 }
00792 /* Computing MAX */
00793                 d__2 = xmax, d__3 = (d__1 = x[j], absMACRO(d__1));
00794                 xmax = maxMACRO(d__2,d__3);
00795 /* L160: */
00796             }
00797         }
00798         *scale /= tscal;
00799     }
00800 
00801 /*     Scale the column norms by 1/TSCAL for return. */
00802 
00803     if (tscal != 1.) {
00804         d__1 = 1. / tscal;
00805         dscal_(n, &d__1, &cnorm[1], &c__1);
00806     }
00807 
00808     return 0;
00809 
00810 /*     End of DLATRS */
00811 
00812 } /* dlatrs_ */
00813 
00814 #undef a_ref
00815 
00816 
00817 #endif

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7