|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZLAPTM( 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 DOUBLE PRECISION ALPHA, BETA 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION D( * ) 00015 COMPLEX*16 B( LDB, * ), E( * ), X( LDX, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZLAPTM 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) DOUBLE PRECISION 00046 * The scalar alpha. ALPHA must be 1. or -1.; otherwise, 00047 * it is assumed to be 0. 00048 * 00049 * D (input) DOUBLE PRECISION array, dimension (N) 00050 * The n diagonal elements of the tridiagonal matrix A. 00051 * 00052 * E (input) COMPLEX*16 array, dimension (N-1) 00053 * The (n-1) subdiagonal or superdiagonal elements of A. 00054 * 00055 * X (input) COMPLEX*16 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) DOUBLE PRECISION 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*16 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 DOUBLE PRECISION ONE, ZERO 00077 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+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 DCONJG 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 ) + DCONJG( 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 ) + DCONJG( 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 $ DCONJG( 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 $ DCONJG( 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 ) - DCONJG( 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 ) - DCONJG( 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 $ DCONJG( 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 $ DCONJG( 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 ZLAPTM 00192 * 00193 END
1.7.3 
