ergo
|
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