crscl.f

Go to the documentation of this file.

*> \brief \b CRSCL multiplies a vector by the reciprocal of a real scalar.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download CRSCL + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/crscl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/crscl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/crscl.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE CRSCL( N, A, X, INCX )
*
*       .. Scalar Arguments ..
*       INTEGER            INCX, N
*       COMPLEX            A
*       ..
*       .. Array Arguments ..
*       COMPLEX            X( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CRSCL multiplies an n-element complex vector x by the complex scalar
*> 1/a.  This is done without overflow or underflow as long as
*> the final result x/a does not overflow or underflow.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of components of the vector x.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is COMPLEX
*>          The scalar a which is used to divide each component of x.
*>          A must not be 0, or the subroutine will divide by zero.
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*>          X is COMPLEX array, dimension
*>                         (1+(N-1)*abs(INCX))
*>          The n-element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*>          INCX is INTEGER
*>          The increment between successive values of the vector X.
*>          > 0:  X(1) = X(1) and X(1+(i-1)*INCX) = x(i),     1< i<= n
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complexOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE crscl( N, A, X, INCX )
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INCX, N
      COMPLEX            A
*     ..
*     .. Array Arguments ..
      COMPLEX            X( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      REAL               SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR
     %                   , UI
*     ..
*     .. External Functions ..
      REAL               SLAMCH
      COMPLEX            CLADIV
      EXTERNAL           slamch, cladiv
*     ..
*     .. External Subroutines ..
      EXTERNAL           cscal, csscal, csrscl
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( n.LE.0 )
     $   RETURN
*
*     Get machine parameters
*
      safmin = slamch( 'S' )
      safmax = one / safmin
      ov   = slamch( 'O' )
*
*     Initialize constants related to A.
*
      ar = real( a )
      ai = aimag( a )
      absr = abs( ar )
      absi = abs( ai )
*
      IF( ai.EQ.zero ) THEN
*        If alpha is real, then we can use csrscl
         CALL csrscl( n, ar, x, incx )
*
      ELSE IF( ar.EQ.zero ) THEN
*        If alpha has a zero real part, then we follow the same rules as if
*        alpha were real.
         IF( absi.GT.safmax ) THEN
            CALL csscal( n, safmin, x, incx )
            CALL cscal( n, cmplx( zero, -safmax / ai ), x, incx )
         ELSE IF( absi.LT.safmin ) THEN
            CALL cscal( n, cmplx( zero, -safmin / ai ), x, incx )
            CALL csscal( n, safmax, x, incx )
         ELSE
            CALL cscal( n, cmplx( zero, -one / ai ), x, incx )
         END IF
*
      ELSE
*        The following numbers can be computed.
*        They are the inverse of the real and imaginary parts of 1/alpha.
*        Note that a and b are always different from zero.
*        NaNs are only possible if either:
*        1. alphaR or alphaI is NaN.
*        2. alphaR and alphaI are both infinite, in which case it makes sense
*        to propagate a NaN.
         ur = ar + ai * ( ai / ar )
         ui = ai + ar * ( ar / ai )
*
         IF( (abs( ur ).LT.safmin).OR.(abs( ui ).LT.safmin) ) THEN
*           This means that both alphaR and alphaI are very small.
            CALL cscal( n, cmplx( safmin / ur, -safmin / ui ), x,
     $                  incx )
            CALL csscal( n, safmax, x, incx )
         ELSE IF( (abs( ur ).GT.safmax).OR.(abs( ui ).GT.safmax) ) THEN
            IF( (absr.GT.ov).OR.(absi.GT.ov) ) THEN
*              This means that a and b are both Inf. No need for scaling.
               CALL cscal( n, cmplx( one / ur, -one / ui ), x, incx )
            ELSE
               CALL csscal( n, safmin, x, incx )
               IF( (abs( ur ).GT.ov).OR.(abs( ui ).GT.ov) ) THEN
*                 Infs were generated. We do proper scaling to avoid them.
                  IF( absr.GE.absi ) THEN
*                    ABS( UR ) <= ABS( UI )
                     ur = (safmin * ar) + safmin * (ai * ( ai / ar ))
                     ui = (safmin * ai) + ar * ( (safmin * ar) / ai )
                  ELSE
*                    ABS( UR ) > ABS( UI )
                     ur = (safmin * ar) + ai * ( (safmin * ai) / ar )
                     ui = (safmin * ai) + safmin * (ar * ( ar / ai ))
                  END IF
                  CALL cscal( n, cmplx( one / ur, -one / ui ), x,
     $                        incx )
               ELSE
                  CALL cscal( n, cmplx( safmax / ur, -safmax / ui ),
     $                        x, incx )
               END IF
            END IF
         ELSE
            CALL cscal( n, cmplx( one / ur, -one / ui ), x, incx )
         END IF
      END IF
*
      RETURN
*
*     End of CRSCL
*
      END