LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clatbs.f
Go to the documentation of this file.
00001 *> \brief \b CLATBS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLATBS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatbs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatbs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatbs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
00022 *                          SCALE, CNORM, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          DIAG, NORMIN, TRANS, UPLO
00026 *       INTEGER            INFO, KD, LDAB, N
00027 *       REAL               SCALE
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       REAL               CNORM( * )
00031 *       COMPLEX            AB( LDAB, * ), X( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> CLATBS solves one of the triangular systems
00041 *>
00042 *>    A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
00043 *>
00044 *> with scaling to prevent overflow, where A is an upper or lower
00045 *> triangular band matrix.  Here A**T denotes the transpose of A, x and b
00046 *> are n-element vectors, and s is a scaling factor, usually less than
00047 *> or equal to 1, chosen so that the components of x will be less than
00048 *> the overflow threshold.  If the unscaled problem will not cause
00049 *> overflow, the Level 2 BLAS routine CTBSV is called.  If the matrix A
00050 *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
00051 *> non-trivial solution to A*x = 0 is returned.
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] UPLO
00058 *> \verbatim
00059 *>          UPLO is CHARACTER*1
00060 *>          Specifies whether the matrix A is upper or lower triangular.
00061 *>          = 'U':  Upper triangular
00062 *>          = 'L':  Lower triangular
00063 *> \endverbatim
00064 *>
00065 *> \param[in] TRANS
00066 *> \verbatim
00067 *>          TRANS is CHARACTER*1
00068 *>          Specifies the operation applied to A.
00069 *>          = 'N':  Solve A * x = s*b     (No transpose)
00070 *>          = 'T':  Solve A**T * x = s*b  (Transpose)
00071 *>          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
00072 *> \endverbatim
00073 *>
00074 *> \param[in] DIAG
00075 *> \verbatim
00076 *>          DIAG is CHARACTER*1
00077 *>          Specifies whether or not the matrix A is unit triangular.
00078 *>          = 'N':  Non-unit triangular
00079 *>          = 'U':  Unit triangular
00080 *> \endverbatim
00081 *>
00082 *> \param[in] NORMIN
00083 *> \verbatim
00084 *>          NORMIN is CHARACTER*1
00085 *>          Specifies whether CNORM has been set or not.
00086 *>          = 'Y':  CNORM contains the column norms on entry
00087 *>          = 'N':  CNORM is not set on entry.  On exit, the norms will
00088 *>                  be computed and stored in CNORM.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] N
00092 *> \verbatim
00093 *>          N is INTEGER
00094 *>          The order of the matrix A.  N >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] KD
00098 *> \verbatim
00099 *>          KD is INTEGER
00100 *>          The number of subdiagonals or superdiagonals in the
00101 *>          triangular matrix A.  KD >= 0.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] AB
00105 *> \verbatim
00106 *>          AB is COMPLEX array, dimension (LDAB,N)
00107 *>          The upper or lower triangular band matrix A, stored in the
00108 *>          first KD+1 rows of the array. The j-th column of A is stored
00109 *>          in the j-th column of the array AB as follows:
00110 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00111 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00112 *> \endverbatim
00113 *>
00114 *> \param[in] LDAB
00115 *> \verbatim
00116 *>          LDAB is INTEGER
00117 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00118 *> \endverbatim
00119 *>
00120 *> \param[in,out] X
00121 *> \verbatim
00122 *>          X is COMPLEX array, dimension (N)
00123 *>          On entry, the right hand side b of the triangular system.
00124 *>          On exit, X is overwritten by the solution vector x.
00125 *> \endverbatim
00126 *>
00127 *> \param[out] SCALE
00128 *> \verbatim
00129 *>          SCALE is REAL
00130 *>          The scaling factor s for the triangular system
00131 *>             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
00132 *>          If SCALE = 0, the matrix A is singular or badly scaled, and
00133 *>          the vector x is an exact or approximate solution to A*x = 0.
00134 *> \endverbatim
00135 *>
00136 *> \param[in,out] CNORM
00137 *> \verbatim
00138 *>          CNORM is REAL array, dimension (N)
00139 *>
00140 *>          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
00141 *>          contains the norm of the off-diagonal part of the j-th column
00142 *>          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
00143 *>          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
00144 *>          must be greater than or equal to the 1-norm.
00145 *>
00146 *>          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
00147 *>          returns the 1-norm of the offdiagonal part of the j-th column
00148 *>          of A.
00149 *> \endverbatim
00150 *>
00151 *> \param[out] INFO
00152 *> \verbatim
00153 *>          INFO is INTEGER
00154 *>          = 0:  successful exit
00155 *>          < 0:  if INFO = -k, the k-th argument had an illegal value
00156 *> \endverbatim
00157 *
00158 *  Authors:
00159 *  ========
00160 *
00161 *> \author Univ. of Tennessee 
00162 *> \author Univ. of California Berkeley 
00163 *> \author Univ. of Colorado Denver 
00164 *> \author NAG Ltd. 
00165 *
00166 *> \date April 2012
00167 *
00168 *> \ingroup complexOTHERauxiliary
00169 *
00170 *> \par Further Details:
00171 *  =====================
00172 *>
00173 *> \verbatim
00174 *>
00175 *>  A rough bound on x is computed; if that is less than overflow, CTBSV
00176 *>  is called, otherwise, specific code is used which checks for possible
00177 *>  overflow or divide-by-zero at every operation.
00178 *>
00179 *>  A columnwise scheme is used for solving A*x = b.  The basic algorithm
00180 *>  if A is lower triangular is
00181 *>
00182 *>       x[1:n] := b[1:n]
00183 *>       for j = 1, ..., n
00184 *>            x(j) := x(j) / A(j,j)
00185 *>            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
00186 *>       end
00187 *>
00188 *>  Define bounds on the components of x after j iterations of the loop:
00189 *>     M(j) = bound on x[1:j]
00190 *>     G(j) = bound on x[j+1:n]
00191 *>  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
00192 *>
00193 *>  Then for iteration j+1 we have
00194 *>     M(j+1) <= G(j) / | A(j+1,j+1) |
00195 *>     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
00196 *>            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
00197 *>
00198 *>  where CNORM(j+1) is greater than or equal to the infinity-norm of
00199 *>  column j+1 of A, not counting the diagonal.  Hence
00200 *>
00201 *>     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
00202 *>                  1<=i<=j
00203 *>  and
00204 *>
00205 *>     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
00206 *>                                   1<=i< j
00207 *>
00208 *>  Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the
00209 *>  reciprocal of the largest M(j), j=1,..,n, is larger than
00210 *>  max(underflow, 1/overflow).
00211 *>
00212 *>  The bound on x(j) is also used to determine when a step in the
00213 *>  columnwise method can be performed without fear of overflow.  If
00214 *>  the computed bound is greater than a large constant, x is scaled to
00215 *>  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
00216 *>  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
00217 *>
00218 *>  Similarly, a row-wise scheme is used to solve A**T *x = b  or
00219 *>  A**H *x = b.  The basic algorithm for A upper triangular is
00220 *>
00221 *>       for j = 1, ..., n
00222 *>            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
00223 *>       end
00224 *>
00225 *>  We simultaneously compute two bounds
00226 *>       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
00227 *>       M(j) = bound on x(i), 1<=i<=j
00228 *>
00229 *>  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
00230 *>  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
00231 *>  Then the bound on x(j) is
00232 *>
00233 *>       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
00234 *>
00235 *>            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
00236 *>                      1<=i<=j
00237 *>
00238 *>  and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater
00239 *>  than max(underflow, 1/overflow).
00240 *> \endverbatim
00241 *>
00242 *  =====================================================================
00243       SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
00244      $                   SCALE, CNORM, INFO )
00245 *
00246 *  -- LAPACK auxiliary routine (version 3.4.1) --
00247 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00248 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00249 *     April 2012
00250 *
00251 *     .. Scalar Arguments ..
00252       CHARACTER          DIAG, NORMIN, TRANS, UPLO
00253       INTEGER            INFO, KD, LDAB, N
00254       REAL               SCALE
00255 *     ..
00256 *     .. Array Arguments ..
00257       REAL               CNORM( * )
00258       COMPLEX            AB( LDAB, * ), X( * )
00259 *     ..
00260 *
00261 *  =====================================================================
00262 *
00263 *     .. Parameters ..
00264       REAL               ZERO, HALF, ONE, TWO
00265       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0,
00266      $                   TWO = 2.0E+0 )
00267 *     ..
00268 *     .. Local Scalars ..
00269       LOGICAL            NOTRAN, NOUNIT, UPPER
00270       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
00271       REAL               BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
00272      $                   XBND, XJ, XMAX
00273       COMPLEX            CSUMJ, TJJS, USCAL, ZDUM
00274 *     ..
00275 *     .. External Functions ..
00276       LOGICAL            LSAME
00277       INTEGER            ICAMAX, ISAMAX
00278       REAL               SCASUM, SLAMCH
00279       COMPLEX            CDOTC, CDOTU, CLADIV
00280       EXTERNAL           LSAME, ICAMAX, ISAMAX, SCASUM, SLAMCH, CDOTC,
00281      $                   CDOTU, CLADIV
00282 *     ..
00283 *     .. External Subroutines ..
00284       EXTERNAL           CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA
00285 *     ..
00286 *     .. Intrinsic Functions ..
00287       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
00288 *     ..
00289 *     .. Statement Functions ..
00290       REAL               CABS1, CABS2
00291 *     ..
00292 *     .. Statement Function definitions ..
00293       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00294       CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2. ) +
00295      $                ABS( AIMAG( ZDUM ) / 2. )
00296 *     ..
00297 *     .. Executable Statements ..
00298 *
00299       INFO = 0
00300       UPPER = LSAME( UPLO, 'U' )
00301       NOTRAN = LSAME( TRANS, 'N' )
00302       NOUNIT = LSAME( DIAG, 'N' )
00303 *
00304 *     Test the input parameters.
00305 *
00306       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00307          INFO = -1
00308       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00309      $         LSAME( TRANS, 'C' ) ) THEN
00310          INFO = -2
00311       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00312          INFO = -3
00313       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
00314      $         LSAME( NORMIN, 'N' ) ) THEN
00315          INFO = -4
00316       ELSE IF( N.LT.0 ) THEN
00317          INFO = -5
00318       ELSE IF( KD.LT.0 ) THEN
00319          INFO = -6
00320       ELSE IF( LDAB.LT.KD+1 ) THEN
00321          INFO = -8
00322       END IF
00323       IF( INFO.NE.0 ) THEN
00324          CALL XERBLA( 'CLATBS', -INFO )
00325          RETURN
00326       END IF
00327 *
00328 *     Quick return if possible
00329 *
00330       IF( N.EQ.0 )
00331      $   RETURN
00332 *
00333 *     Determine machine dependent parameters to control overflow.
00334 *
00335       SMLNUM = SLAMCH( 'Safe minimum' )
00336       BIGNUM = ONE / SMLNUM
00337       CALL SLABAD( SMLNUM, BIGNUM )
00338       SMLNUM = SMLNUM / SLAMCH( 'Precision' )
00339       BIGNUM = ONE / SMLNUM
00340       SCALE = ONE
00341 *
00342       IF( LSAME( NORMIN, 'N' ) ) THEN
00343 *
00344 *        Compute the 1-norm of each column, not including the diagonal.
00345 *
00346          IF( UPPER ) THEN
00347 *
00348 *           A is upper triangular.
00349 *
00350             DO 10 J = 1, N
00351                JLEN = MIN( KD, J-1 )
00352                CNORM( J ) = SCASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
00353    10       CONTINUE
00354          ELSE
00355 *
00356 *           A is lower triangular.
00357 *
00358             DO 20 J = 1, N
00359                JLEN = MIN( KD, N-J )
00360                IF( JLEN.GT.0 ) THEN
00361                   CNORM( J ) = SCASUM( JLEN, AB( 2, J ), 1 )
00362                ELSE
00363                   CNORM( J ) = ZERO
00364                END IF
00365    20       CONTINUE
00366          END IF
00367       END IF
00368 *
00369 *     Scale the column norms by TSCAL if the maximum element in CNORM is
00370 *     greater than BIGNUM/2.
00371 *
00372       IMAX = ISAMAX( N, CNORM, 1 )
00373       TMAX = CNORM( IMAX )
00374       IF( TMAX.LE.BIGNUM*HALF ) THEN
00375          TSCAL = ONE
00376       ELSE
00377          TSCAL = HALF / ( SMLNUM*TMAX )
00378          CALL SSCAL( N, TSCAL, CNORM, 1 )
00379       END IF
00380 *
00381 *     Compute a bound on the computed solution vector to see if the
00382 *     Level 2 BLAS routine CTBSV can be used.
00383 *
00384       XMAX = ZERO
00385       DO 30 J = 1, N
00386          XMAX = MAX( XMAX, CABS2( X( J ) ) )
00387    30 CONTINUE
00388       XBND = XMAX
00389       IF( NOTRAN ) THEN
00390 *
00391 *        Compute the growth in A * x = b.
00392 *
00393          IF( UPPER ) THEN
00394             JFIRST = N
00395             JLAST = 1
00396             JINC = -1
00397             MAIND = KD + 1
00398          ELSE
00399             JFIRST = 1
00400             JLAST = N
00401             JINC = 1
00402             MAIND = 1
00403          END IF
00404 *
00405          IF( TSCAL.NE.ONE ) THEN
00406             GROW = ZERO
00407             GO TO 60
00408          END IF
00409 *
00410          IF( NOUNIT ) THEN
00411 *
00412 *           A is non-unit triangular.
00413 *
00414 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
00415 *           Initially, G(0) = max{x(i), i=1,...,n}.
00416 *
00417             GROW = HALF / MAX( XBND, SMLNUM )
00418             XBND = GROW
00419             DO 40 J = JFIRST, JLAST, JINC
00420 *
00421 *              Exit the loop if the growth factor is too small.
00422 *
00423                IF( GROW.LE.SMLNUM )
00424      $            GO TO 60
00425 *
00426                TJJS = AB( MAIND, J )
00427                TJJ = CABS1( TJJS )
00428 *
00429                IF( TJJ.GE.SMLNUM ) THEN
00430 *
00431 *                 M(j) = G(j-1) / abs(A(j,j))
00432 *
00433                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
00434                ELSE
00435 *
00436 *                 M(j) could overflow, set XBND to 0.
00437 *
00438                   XBND = ZERO
00439                END IF
00440 *
00441                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
00442 *
00443 *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
00444 *
00445                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
00446                ELSE
00447 *
00448 *                 G(j) could overflow, set GROW to 0.
00449 *
00450                   GROW = ZERO
00451                END IF
00452    40       CONTINUE
00453             GROW = XBND
00454          ELSE
00455 *
00456 *           A is unit triangular.
00457 *
00458 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
00459 *
00460             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
00461             DO 50 J = JFIRST, JLAST, JINC
00462 *
00463 *              Exit the loop if the growth factor is too small.
00464 *
00465                IF( GROW.LE.SMLNUM )
00466      $            GO TO 60
00467 *
00468 *              G(j) = G(j-1)*( 1 + CNORM(j) )
00469 *
00470                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
00471    50       CONTINUE
00472          END IF
00473    60    CONTINUE
00474 *
00475       ELSE
00476 *
00477 *        Compute the growth in A**T * x = b  or  A**H * x = b.
00478 *
00479          IF( UPPER ) THEN
00480             JFIRST = 1
00481             JLAST = N
00482             JINC = 1
00483             MAIND = KD + 1
00484          ELSE
00485             JFIRST = N
00486             JLAST = 1
00487             JINC = -1
00488             MAIND = 1
00489          END IF
00490 *
00491          IF( TSCAL.NE.ONE ) THEN
00492             GROW = ZERO
00493             GO TO 90
00494          END IF
00495 *
00496          IF( NOUNIT ) THEN
00497 *
00498 *           A is non-unit triangular.
00499 *
00500 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
00501 *           Initially, M(0) = max{x(i), i=1,...,n}.
00502 *
00503             GROW = HALF / MAX( XBND, SMLNUM )
00504             XBND = GROW
00505             DO 70 J = JFIRST, JLAST, JINC
00506 *
00507 *              Exit the loop if the growth factor is too small.
00508 *
00509                IF( GROW.LE.SMLNUM )
00510      $            GO TO 90
00511 *
00512 *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
00513 *
00514                XJ = ONE + CNORM( J )
00515                GROW = MIN( GROW, XBND / XJ )
00516 *
00517                TJJS = AB( MAIND, J )
00518                TJJ = CABS1( TJJS )
00519 *
00520                IF( TJJ.GE.SMLNUM ) THEN
00521 *
00522 *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
00523 *
00524                   IF( XJ.GT.TJJ )
00525      $               XBND = XBND*( TJJ / XJ )
00526                ELSE
00527 *
00528 *                 M(j) could overflow, set XBND to 0.
00529 *
00530                   XBND = ZERO
00531                END IF
00532    70       CONTINUE
00533             GROW = MIN( GROW, XBND )
00534          ELSE
00535 *
00536 *           A is unit triangular.
00537 *
00538 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
00539 *
00540             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
00541             DO 80 J = JFIRST, JLAST, JINC
00542 *
00543 *              Exit the loop if the growth factor is too small.
00544 *
00545                IF( GROW.LE.SMLNUM )
00546      $            GO TO 90
00547 *
00548 *              G(j) = ( 1 + CNORM(j) )*G(j-1)
00549 *
00550                XJ = ONE + CNORM( J )
00551                GROW = GROW / XJ
00552    80       CONTINUE
00553          END IF
00554    90    CONTINUE
00555       END IF
00556 *
00557       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
00558 *
00559 *        Use the Level 2 BLAS solve if the reciprocal of the bound on
00560 *        elements of X is not too small.
00561 *
00562          CALL CTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
00563       ELSE
00564 *
00565 *        Use a Level 1 BLAS solve, scaling intermediate results.
00566 *
00567          IF( XMAX.GT.BIGNUM*HALF ) THEN
00568 *
00569 *           Scale X so that its components are less than or equal to
00570 *           BIGNUM in absolute value.
00571 *
00572             SCALE = ( BIGNUM*HALF ) / XMAX
00573             CALL CSSCAL( N, SCALE, X, 1 )
00574             XMAX = BIGNUM
00575          ELSE
00576             XMAX = XMAX*TWO
00577          END IF
00578 *
00579          IF( NOTRAN ) THEN
00580 *
00581 *           Solve A * x = b
00582 *
00583             DO 110 J = JFIRST, JLAST, JINC
00584 *
00585 *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
00586 *
00587                XJ = CABS1( X( J ) )
00588                IF( NOUNIT ) THEN
00589                   TJJS = AB( MAIND, J )*TSCAL
00590                ELSE
00591                   TJJS = TSCAL
00592                   IF( TSCAL.EQ.ONE )
00593      $               GO TO 105
00594                END IF
00595                   TJJ = CABS1( TJJS )
00596                   IF( TJJ.GT.SMLNUM ) THEN
00597 *
00598 *                    abs(A(j,j)) > SMLNUM:
00599 *
00600                      IF( TJJ.LT.ONE ) THEN
00601                         IF( XJ.GT.TJJ*BIGNUM ) THEN
00602 *
00603 *                          Scale x by 1/b(j).
00604 *
00605                            REC = ONE / XJ
00606                            CALL CSSCAL( N, REC, X, 1 )
00607                            SCALE = SCALE*REC
00608                            XMAX = XMAX*REC
00609                         END IF
00610                      END IF
00611                      X( J ) = CLADIV( X( J ), TJJS )
00612                      XJ = CABS1( X( J ) )
00613                   ELSE IF( TJJ.GT.ZERO ) THEN
00614 *
00615 *                    0 < abs(A(j,j)) <= SMLNUM:
00616 *
00617                      IF( XJ.GT.TJJ*BIGNUM ) THEN
00618 *
00619 *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
00620 *                       to avoid overflow when dividing by A(j,j).
00621 *
00622                         REC = ( TJJ*BIGNUM ) / XJ
00623                         IF( CNORM( J ).GT.ONE ) THEN
00624 *
00625 *                          Scale by 1/CNORM(j) to avoid overflow when
00626 *                          multiplying x(j) times column j.
00627 *
00628                            REC = REC / CNORM( J )
00629                         END IF
00630                         CALL CSSCAL( N, REC, X, 1 )
00631                         SCALE = SCALE*REC
00632                         XMAX = XMAX*REC
00633                      END IF
00634                      X( J ) = CLADIV( X( J ), TJJS )
00635                      XJ = CABS1( X( J ) )
00636                   ELSE
00637 *
00638 *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
00639 *                    scale = 0, and compute a solution to A*x = 0.
00640 *
00641                      DO 100 I = 1, N
00642                         X( I ) = ZERO
00643   100                CONTINUE
00644                      X( J ) = ONE
00645                      XJ = ONE
00646                      SCALE = ZERO
00647                      XMAX = ZERO
00648                   END IF
00649   105          CONTINUE
00650 *
00651 *              Scale x if necessary to avoid overflow when adding a
00652 *              multiple of column j of A.
00653 *
00654                IF( XJ.GT.ONE ) THEN
00655                   REC = ONE / XJ
00656                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
00657 *
00658 *                    Scale x by 1/(2*abs(x(j))).
00659 *
00660                      REC = REC*HALF
00661                      CALL CSSCAL( N, REC, X, 1 )
00662                      SCALE = SCALE*REC
00663                   END IF
00664                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
00665 *
00666 *                 Scale x by 1/2.
00667 *
00668                   CALL CSSCAL( N, HALF, X, 1 )
00669                   SCALE = SCALE*HALF
00670                END IF
00671 *
00672                IF( UPPER ) THEN
00673                   IF( J.GT.1 ) THEN
00674 *
00675 *                    Compute the update
00676 *                       x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
00677 *                                             x(j)* A(max(1,j-kd):j-1,j)
00678 *
00679                      JLEN = MIN( KD, J-1 )
00680                      CALL CAXPY( JLEN, -X( J )*TSCAL,
00681      $                           AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
00682                      I = ICAMAX( J-1, X, 1 )
00683                      XMAX = CABS1( X( I ) )
00684                   END IF
00685                ELSE IF( J.LT.N ) THEN
00686 *
00687 *                 Compute the update
00688 *                    x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
00689 *                                          x(j) * A(j+1:min(j+kd,n),j)
00690 *
00691                   JLEN = MIN( KD, N-J )
00692                   IF( JLEN.GT.0 )
00693      $               CALL CAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
00694      $                           X( J+1 ), 1 )
00695                   I = J + ICAMAX( N-J, X( J+1 ), 1 )
00696                   XMAX = CABS1( X( I ) )
00697                END IF
00698   110       CONTINUE
00699 *
00700          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
00701 *
00702 *           Solve A**T * x = b
00703 *
00704             DO 150 J = JFIRST, JLAST, JINC
00705 *
00706 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
00707 *                                    k<>j
00708 *
00709                XJ = CABS1( X( J ) )
00710                USCAL = TSCAL
00711                REC = ONE / MAX( XMAX, ONE )
00712                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
00713 *
00714 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
00715 *
00716                   REC = REC*HALF
00717                   IF( NOUNIT ) THEN
00718                      TJJS = AB( MAIND, J )*TSCAL
00719                   ELSE
00720                      TJJS = TSCAL
00721                   END IF
00722                      TJJ = CABS1( TJJS )
00723                      IF( TJJ.GT.ONE ) THEN
00724 *
00725 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
00726 *
00727                         REC = MIN( ONE, REC*TJJ )
00728                         USCAL = CLADIV( USCAL, TJJS )
00729                      END IF
00730                   IF( REC.LT.ONE ) THEN
00731                      CALL CSSCAL( N, REC, X, 1 )
00732                      SCALE = SCALE*REC
00733                      XMAX = XMAX*REC
00734                   END IF
00735                END IF
00736 *
00737                CSUMJ = ZERO
00738                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
00739 *
00740 *                 If the scaling needed for A in the dot product is 1,
00741 *                 call CDOTU to perform the dot product.
00742 *
00743                   IF( UPPER ) THEN
00744                      JLEN = MIN( KD, J-1 )
00745                      CSUMJ = CDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
00746      $                       X( J-JLEN ), 1 )
00747                   ELSE
00748                      JLEN = MIN( KD, N-J )
00749                      IF( JLEN.GT.1 )
00750      $                  CSUMJ = CDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
00751      $                          1 )
00752                   END IF
00753                ELSE
00754 *
00755 *                 Otherwise, use in-line code for the dot product.
00756 *
00757                   IF( UPPER ) THEN
00758                      JLEN = MIN( KD, J-1 )
00759                      DO 120 I = 1, JLEN
00760                         CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
00761      $                          X( J-JLEN-1+I )
00762   120                CONTINUE
00763                   ELSE
00764                      JLEN = MIN( KD, N-J )
00765                      DO 130 I = 1, JLEN
00766                         CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
00767   130                CONTINUE
00768                   END IF
00769                END IF
00770 *
00771                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
00772 *
00773 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
00774 *                 was not used to scale the dotproduct.
00775 *
00776                   X( J ) = X( J ) - CSUMJ
00777                   XJ = CABS1( X( J ) )
00778                   IF( NOUNIT ) THEN
00779 *
00780 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
00781 *
00782                      TJJS = AB( MAIND, J )*TSCAL
00783                   ELSE
00784                      TJJS = TSCAL
00785                      IF( TSCAL.EQ.ONE )
00786      $                  GO TO 145
00787                   END IF
00788                      TJJ = CABS1( TJJS )
00789                      IF( TJJ.GT.SMLNUM ) THEN
00790 *
00791 *                       abs(A(j,j)) > SMLNUM:
00792 *
00793                         IF( TJJ.LT.ONE ) THEN
00794                            IF( XJ.GT.TJJ*BIGNUM ) THEN
00795 *
00796 *                             Scale X by 1/abs(x(j)).
00797 *
00798                               REC = ONE / XJ
00799                               CALL CSSCAL( N, REC, X, 1 )
00800                               SCALE = SCALE*REC
00801                               XMAX = XMAX*REC
00802                            END IF
00803                         END IF
00804                         X( J ) = CLADIV( X( J ), TJJS )
00805                      ELSE IF( TJJ.GT.ZERO ) THEN
00806 *
00807 *                       0 < abs(A(j,j)) <= SMLNUM:
00808 *
00809                         IF( XJ.GT.TJJ*BIGNUM ) THEN
00810 *
00811 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
00812 *
00813                            REC = ( TJJ*BIGNUM ) / XJ
00814                            CALL CSSCAL( N, REC, X, 1 )
00815                            SCALE = SCALE*REC
00816                            XMAX = XMAX*REC
00817                         END IF
00818                         X( J ) = CLADIV( X( J ), TJJS )
00819                      ELSE
00820 *
00821 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
00822 *                       scale = 0 and compute a solution to A**T *x = 0.
00823 *
00824                         DO 140 I = 1, N
00825                            X( I ) = ZERO
00826   140                   CONTINUE
00827                         X( J ) = ONE
00828                         SCALE = ZERO
00829                         XMAX = ZERO
00830                      END IF
00831   145             CONTINUE
00832                ELSE
00833 *
00834 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
00835 *                 product has already been divided by 1/A(j,j).
00836 *
00837                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
00838                END IF
00839                XMAX = MAX( XMAX, CABS1( X( J ) ) )
00840   150       CONTINUE
00841 *
00842          ELSE
00843 *
00844 *           Solve A**H * x = b
00845 *
00846             DO 190 J = JFIRST, JLAST, JINC
00847 *
00848 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
00849 *                                    k<>j
00850 *
00851                XJ = CABS1( X( J ) )
00852                USCAL = TSCAL
00853                REC = ONE / MAX( XMAX, ONE )
00854                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
00855 *
00856 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
00857 *
00858                   REC = REC*HALF
00859                   IF( NOUNIT ) THEN
00860                      TJJS = CONJG( AB( MAIND, J ) )*TSCAL
00861                   ELSE
00862                      TJJS = TSCAL
00863                   END IF
00864                      TJJ = CABS1( TJJS )
00865                      IF( TJJ.GT.ONE ) THEN
00866 *
00867 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
00868 *
00869                         REC = MIN( ONE, REC*TJJ )
00870                         USCAL = CLADIV( USCAL, TJJS )
00871                      END IF
00872                   IF( REC.LT.ONE ) THEN
00873                      CALL CSSCAL( N, REC, X, 1 )
00874                      SCALE = SCALE*REC
00875                      XMAX = XMAX*REC
00876                   END IF
00877                END IF
00878 *
00879                CSUMJ = ZERO
00880                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
00881 *
00882 *                 If the scaling needed for A in the dot product is 1,
00883 *                 call CDOTC to perform the dot product.
00884 *
00885                   IF( UPPER ) THEN
00886                      JLEN = MIN( KD, J-1 )
00887                      CSUMJ = CDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
00888      $                       X( J-JLEN ), 1 )
00889                   ELSE
00890                      JLEN = MIN( KD, N-J )
00891                      IF( JLEN.GT.1 )
00892      $                  CSUMJ = CDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
00893      $                          1 )
00894                   END IF
00895                ELSE
00896 *
00897 *                 Otherwise, use in-line code for the dot product.
00898 *
00899                   IF( UPPER ) THEN
00900                      JLEN = MIN( KD, J-1 )
00901                      DO 160 I = 1, JLEN
00902                         CSUMJ = CSUMJ + ( CONJG( AB( KD+I-JLEN, J ) )*
00903      $                          USCAL )*X( J-JLEN-1+I )
00904   160                CONTINUE
00905                   ELSE
00906                      JLEN = MIN( KD, N-J )
00907                      DO 170 I = 1, JLEN
00908                         CSUMJ = CSUMJ + ( CONJG( AB( I+1, J ) )*USCAL )*
00909      $                          X( J+I )
00910   170                CONTINUE
00911                   END IF
00912                END IF
00913 *
00914                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
00915 *
00916 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
00917 *                 was not used to scale the dotproduct.
00918 *
00919                   X( J ) = X( J ) - CSUMJ
00920                   XJ = CABS1( X( J ) )
00921                   IF( NOUNIT ) THEN
00922 *
00923 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
00924 *
00925                      TJJS = CONJG( AB( MAIND, J ) )*TSCAL
00926                   ELSE
00927                      TJJS = TSCAL
00928                      IF( TSCAL.EQ.ONE )
00929      $                  GO TO 185
00930                   END IF
00931                      TJJ = CABS1( TJJS )
00932                      IF( TJJ.GT.SMLNUM ) THEN
00933 *
00934 *                       abs(A(j,j)) > SMLNUM:
00935 *
00936                         IF( TJJ.LT.ONE ) THEN
00937                            IF( XJ.GT.TJJ*BIGNUM ) THEN
00938 *
00939 *                             Scale X by 1/abs(x(j)).
00940 *
00941                               REC = ONE / XJ
00942                               CALL CSSCAL( N, REC, X, 1 )
00943                               SCALE = SCALE*REC
00944                               XMAX = XMAX*REC
00945                            END IF
00946                         END IF
00947                         X( J ) = CLADIV( X( J ), TJJS )
00948                      ELSE IF( TJJ.GT.ZERO ) THEN
00949 *
00950 *                       0 < abs(A(j,j)) <= SMLNUM:
00951 *
00952                         IF( XJ.GT.TJJ*BIGNUM ) THEN
00953 *
00954 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
00955 *
00956                            REC = ( TJJ*BIGNUM ) / XJ
00957                            CALL CSSCAL( N, REC, X, 1 )
00958                            SCALE = SCALE*REC
00959                            XMAX = XMAX*REC
00960                         END IF
00961                         X( J ) = CLADIV( X( J ), TJJS )
00962                      ELSE
00963 *
00964 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
00965 *                       scale = 0 and compute a solution to A**H *x = 0.
00966 *
00967                         DO 180 I = 1, N
00968                            X( I ) = ZERO
00969   180                   CONTINUE
00970                         X( J ) = ONE
00971                         SCALE = ZERO
00972                         XMAX = ZERO
00973                      END IF
00974   185             CONTINUE
00975                ELSE
00976 *
00977 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
00978 *                 product has already been divided by 1/A(j,j).
00979 *
00980                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
00981                END IF
00982                XMAX = MAX( XMAX, CABS1( X( J ) ) )
00983   190       CONTINUE
00984          END IF
00985          SCALE = SCALE / TSCAL
00986       END IF
00987 *
00988 *     Scale the column norms by 1/TSCAL for return.
00989 *
00990       IF( TSCAL.NE.ONE ) THEN
00991          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
00992       END IF
00993 *
00994       RETURN
00995 *
00996 *     End of CLATBS
00997 *
00998       END
 All Files Functions