LAPACK: TESTING/LIN/claptm.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAPTM( UPLO, N, NRHS, ALPHA, D, E, X, LDX, BETA, B,
00002      $                   LDB )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            LDB, LDX, N, NRHS
00011       REAL               ALPHA, BETA
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               D( * )
00015       COMPLEX            B( LDB, * ), E( * ), X( LDX, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CLAPTM multiplies an N by NRHS matrix X by a Hermitian tridiagonal
00022 *  matrix A and stores the result in a matrix B.  The operation has the
00023 *  form
00024 *
00025 *     B := alpha * A * X + beta * B
00026 *
00027 *  where alpha may be either 1. or -1. and beta may be 0., 1., or -1.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  UPLO    (input) CHARACTER
00033 *          Specifies whether the superdiagonal or the subdiagonal of the
00034 *          tridiagonal matrix A is stored.
00035 *          = 'U':  Upper, E is the superdiagonal of A.
00036 *          = 'L':  Lower, E is the subdiagonal of A.
00037 *
00038 *  N       (input) INTEGER
00039 *          The order of the matrix A.  N >= 0.
00040 *
00041 *  NRHS    (input) INTEGER
00042 *          The number of right hand sides, i.e., the number of columns
00043 *          of the matrices X and B.
00044 *
00045 *  ALPHA   (input) REAL
00046 *          The scalar alpha.  ALPHA must be 1. or -1.; otherwise,
00047 *          it is assumed to be 0.
00048 *
00049 *  D       (input) REAL array, dimension (N)
00050 *          The n diagonal elements of the tridiagonal matrix A.
00051 *
00052 *  E       (input) COMPLEX array, dimension (N-1)
00053 *          The (n-1) subdiagonal or superdiagonal elements of A.
00054 *
00055 *  X       (input) COMPLEX array, dimension (LDX,NRHS)
00056 *          The N by NRHS matrix X.
00057 *
00058 *  LDX     (input) INTEGER
00059 *          The leading dimension of the array X.  LDX >= max(N,1).
00060 *
00061 *  BETA    (input) REAL
00062 *          The scalar beta.  BETA must be 0., 1., or -1.; otherwise,
00063 *          it is assumed to be 1.
00064 *
00065 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00066 *          On entry, the N by NRHS matrix B.
00067 *          On exit, B is overwritten by the matrix expression
00068 *          B := alpha * A * X + beta * B.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of the array B.  LDB >= max(N,1).
00072 *
00073 *  =====================================================================
00074 *
00075 *     .. Parameters ..
00076       REAL               ONE, ZERO
00077       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00078 *     ..
00079 *     .. Local Scalars ..
00080       INTEGER            I, J
00081 *     ..
00082 *     .. External Functions ..
00083       LOGICAL            LSAME
00084       EXTERNAL           LSAME
00085 *     ..
00086 *     .. Intrinsic Functions ..
00087       INTRINSIC          CONJG
00088 *     ..
00089 *     .. Executable Statements ..
00090 *
00091       IF( N.EQ.0 )
00092      $   RETURN
00093 *
00094       IF( BETA.EQ.ZERO ) THEN
00095          DO 20 J = 1, NRHS
00096             DO 10 I = 1, N
00097                B( I, J ) = ZERO
00098    10       CONTINUE
00099    20    CONTINUE
00100       ELSE IF( BETA.EQ.-ONE ) THEN
00101          DO 40 J = 1, NRHS
00102             DO 30 I = 1, N
00103                B( I, J ) = -B( I, J )
00104    30       CONTINUE
00105    40    CONTINUE
00106       END IF
00107 *
00108       IF( ALPHA.EQ.ONE ) THEN
00109          IF( LSAME( UPLO, 'U' ) ) THEN
00110 *
00111 *           Compute B := B + A*X, where E is the superdiagonal of A.
00112 *
00113             DO 60 J = 1, NRHS
00114                IF( N.EQ.1 ) THEN
00115                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
00116                ELSE
00117                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
00118      $                        E( 1 )*X( 2, J )
00119                   B( N, J ) = B( N, J ) + CONJG( E( N-1 ) )*
00120      $                        X( N-1, J ) + D( N )*X( N, J )
00121                   DO 50 I = 2, N - 1
00122                      B( I, J ) = B( I, J ) + CONJG( E( I-1 ) )*
00123      $                           X( I-1, J ) + D( I )*X( I, J ) +
00124      $                           E( I )*X( I+1, J )
00125    50             CONTINUE
00126                END IF
00127    60       CONTINUE
00128          ELSE
00129 *
00130 *           Compute B := B + A*X, where E is the subdiagonal of A.
00131 *
00132             DO 80 J = 1, NRHS
00133                IF( N.EQ.1 ) THEN
00134                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
00135                ELSE
00136                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
00137      $                        CONJG( E( 1 ) )*X( 2, J )
00138                   B( N, J ) = B( N, J ) + E( N-1 )*X( N-1, J ) +
00139      $                        D( N )*X( N, J )
00140                   DO 70 I = 2, N - 1
00141                      B( I, J ) = B( I, J ) + E( I-1 )*X( I-1, J ) +
00142      $                           D( I )*X( I, J ) +
00143      $                           CONJG( E( I ) )*X( I+1, J )
00144    70             CONTINUE
00145                END IF
00146    80       CONTINUE
00147          END IF
00148       ELSE IF( ALPHA.EQ.-ONE ) THEN
00149          IF( LSAME( UPLO, 'U' ) ) THEN
00150 *
00151 *           Compute B := B - A*X, where E is the superdiagonal of A.
00152 *
00153             DO 100 J = 1, NRHS
00154                IF( N.EQ.1 ) THEN
00155                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
00156                ELSE
00157                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
00158      $                        E( 1 )*X( 2, J )
00159                   B( N, J ) = B( N, J ) - CONJG( E( N-1 ) )*
00160      $                        X( N-1, J ) - D( N )*X( N, J )
00161                   DO 90 I = 2, N - 1
00162                      B( I, J ) = B( I, J ) - CONJG( E( I-1 ) )*
00163      $                           X( I-1, J ) - D( I )*X( I, J ) -
00164      $                           E( I )*X( I+1, J )
00165    90             CONTINUE
00166                END IF
00167   100       CONTINUE
00168          ELSE
00169 *
00170 *           Compute B := B - A*X, where E is the subdiagonal of A.
00171 *
00172             DO 120 J = 1, NRHS
00173                IF( N.EQ.1 ) THEN
00174                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
00175                ELSE
00176                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
00177      $                        CONJG( E( 1 ) )*X( 2, J )
00178                   B( N, J ) = B( N, J ) - E( N-1 )*X( N-1, J ) -
00179      $                        D( N )*X( N, J )
00180                   DO 110 I = 2, N - 1
00181                      B( I, J ) = B( I, J ) - E( I-1 )*X( I-1, J ) -
00182      $                           D( I )*X( I, J ) -
00183      $                           CONJG( E( I ) )*X( I+1, J )
00184   110             CONTINUE
00185                END IF
00186   120       CONTINUE
00187          END IF
00188       END IF
00189       RETURN
00190 *
00191 *     End of CLAPTM
00192 *
00193       END