1 SUBROUTINE ZLAPTM( UPLO, N, NRHS, ALPHA, D, E, X, LDX, BETA, B,
2 $ LDB )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDB, LDX, N, NRHS
11 DOUBLE PRECISION ALPHA, BETA
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION D( * )
15 COMPLEX*16 B( LDB, * ), E( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZLAPTM multiplies an N by NRHS matrix X by a Hermitian tridiagonal
22 * matrix A and stores the result in a matrix B. The operation has the
23 * form
24 *
25 * B := alpha * A * X + beta * B
26 *
27 * where alpha may be either 1. or -1. and beta may be 0., 1., or -1.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER
33 * Specifies whether the superdiagonal or the subdiagonal of the
34 * tridiagonal matrix A is stored.
35 * = 'U': Upper, E is the superdiagonal of A.
36 * = 'L': Lower, E is the subdiagonal of A.
37 *
38 * N (input) INTEGER
39 * The order of the matrix A. N >= 0.
40 *
41 * NRHS (input) INTEGER
42 * The number of right hand sides, i.e., the number of columns
43 * of the matrices X and B.
44 *
45 * ALPHA (input) DOUBLE PRECISION
46 * The scalar alpha. ALPHA must be 1. or -1.; otherwise,
47 * it is assumed to be 0.
48 *
49 * D (input) DOUBLE PRECISION array, dimension (N)
50 * The n diagonal elements of the tridiagonal matrix A.
51 *
52 * E (input) COMPLEX*16 array, dimension (N-1)
53 * The (n-1) subdiagonal or superdiagonal elements of A.
54 *
55 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
56 * The N by NRHS matrix X.
57 *
58 * LDX (input) INTEGER
59 * The leading dimension of the array X. LDX >= max(N,1).
60 *
61 * BETA (input) DOUBLE PRECISION
62 * The scalar beta. BETA must be 0., 1., or -1.; otherwise,
63 * it is assumed to be 1.
64 *
65 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
66 * On entry, the N by NRHS matrix B.
67 * On exit, B is overwritten by the matrix expression
68 * B := alpha * A * X + beta * B.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(N,1).
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ONE, ZERO
77 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, J
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 EXTERNAL LSAME
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC DCONJG
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( N.EQ.0 )
92 $ RETURN
93 *
94 IF( BETA.EQ.ZERO ) THEN
95 DO 20 J = 1, NRHS
96 DO 10 I = 1, N
97 B( I, J ) = ZERO
98 10 CONTINUE
99 20 CONTINUE
100 ELSE IF( BETA.EQ.-ONE ) THEN
101 DO 40 J = 1, NRHS
102 DO 30 I = 1, N
103 B( I, J ) = -B( I, J )
104 30 CONTINUE
105 40 CONTINUE
106 END IF
107 *
108 IF( ALPHA.EQ.ONE ) THEN
109 IF( LSAME( UPLO, 'U' ) ) THEN
110 *
111 * Compute B := B + A*X, where E is the superdiagonal of A.
112 *
113 DO 60 J = 1, NRHS
114 IF( N.EQ.1 ) THEN
115 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
116 ELSE
117 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
118 $ E( 1 )*X( 2, J )
119 B( N, J ) = B( N, J ) + DCONJG( E( N-1 ) )*
120 $ X( N-1, J ) + D( N )*X( N, J )
121 DO 50 I = 2, N - 1
122 B( I, J ) = B( I, J ) + DCONJG( E( I-1 ) )*
123 $ X( I-1, J ) + D( I )*X( I, J ) +
124 $ E( I )*X( I+1, J )
125 50 CONTINUE
126 END IF
127 60 CONTINUE
128 ELSE
129 *
130 * Compute B := B + A*X, where E is the subdiagonal of A.
131 *
132 DO 80 J = 1, NRHS
133 IF( N.EQ.1 ) THEN
134 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
135 ELSE
136 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
137 $ DCONJG( E( 1 ) )*X( 2, J )
138 B( N, J ) = B( N, J ) + E( N-1 )*X( N-1, J ) +
139 $ D( N )*X( N, J )
140 DO 70 I = 2, N - 1
141 B( I, J ) = B( I, J ) + E( I-1 )*X( I-1, J ) +
142 $ D( I )*X( I, J ) +
143 $ DCONJG( E( I ) )*X( I+1, J )
144 70 CONTINUE
145 END IF
146 80 CONTINUE
147 END IF
148 ELSE IF( ALPHA.EQ.-ONE ) THEN
149 IF( LSAME( UPLO, 'U' ) ) THEN
150 *
151 * Compute B := B - A*X, where E is the superdiagonal of A.
152 *
153 DO 100 J = 1, NRHS
154 IF( N.EQ.1 ) THEN
155 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
156 ELSE
157 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
158 $ E( 1 )*X( 2, J )
159 B( N, J ) = B( N, J ) - DCONJG( E( N-1 ) )*
160 $ X( N-1, J ) - D( N )*X( N, J )
161 DO 90 I = 2, N - 1
162 B( I, J ) = B( I, J ) - DCONJG( E( I-1 ) )*
163 $ X( I-1, J ) - D( I )*X( I, J ) -
164 $ E( I )*X( I+1, J )
165 90 CONTINUE
166 END IF
167 100 CONTINUE
168 ELSE
169 *
170 * Compute B := B - A*X, where E is the subdiagonal of A.
171 *
172 DO 120 J = 1, NRHS
173 IF( N.EQ.1 ) THEN
174 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
175 ELSE
176 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
177 $ DCONJG( E( 1 ) )*X( 2, J )
178 B( N, J ) = B( N, J ) - E( N-1 )*X( N-1, J ) -
179 $ D( N )*X( N, J )
180 DO 110 I = 2, N - 1
181 B( I, J ) = B( I, J ) - E( I-1 )*X( I-1, J ) -
182 $ D( I )*X( I, J ) -
183 $ DCONJG( E( I ) )*X( I+1, J )
184 110 CONTINUE
185 END IF
186 120 CONTINUE
187 END IF
188 END IF
189 RETURN
190 *
191 * End of ZLAPTM
192 *
193 END
2 $ LDB )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDB, LDX, N, NRHS
11 DOUBLE PRECISION ALPHA, BETA
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION D( * )
15 COMPLEX*16 B( LDB, * ), E( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZLAPTM multiplies an N by NRHS matrix X by a Hermitian tridiagonal
22 * matrix A and stores the result in a matrix B. The operation has the
23 * form
24 *
25 * B := alpha * A * X + beta * B
26 *
27 * where alpha may be either 1. or -1. and beta may be 0., 1., or -1.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER
33 * Specifies whether the superdiagonal or the subdiagonal of the
34 * tridiagonal matrix A is stored.
35 * = 'U': Upper, E is the superdiagonal of A.
36 * = 'L': Lower, E is the subdiagonal of A.
37 *
38 * N (input) INTEGER
39 * The order of the matrix A. N >= 0.
40 *
41 * NRHS (input) INTEGER
42 * The number of right hand sides, i.e., the number of columns
43 * of the matrices X and B.
44 *
45 * ALPHA (input) DOUBLE PRECISION
46 * The scalar alpha. ALPHA must be 1. or -1.; otherwise,
47 * it is assumed to be 0.
48 *
49 * D (input) DOUBLE PRECISION array, dimension (N)
50 * The n diagonal elements of the tridiagonal matrix A.
51 *
52 * E (input) COMPLEX*16 array, dimension (N-1)
53 * The (n-1) subdiagonal or superdiagonal elements of A.
54 *
55 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
56 * The N by NRHS matrix X.
57 *
58 * LDX (input) INTEGER
59 * The leading dimension of the array X. LDX >= max(N,1).
60 *
61 * BETA (input) DOUBLE PRECISION
62 * The scalar beta. BETA must be 0., 1., or -1.; otherwise,
63 * it is assumed to be 1.
64 *
65 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
66 * On entry, the N by NRHS matrix B.
67 * On exit, B is overwritten by the matrix expression
68 * B := alpha * A * X + beta * B.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(N,1).
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ONE, ZERO
77 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, J
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 EXTERNAL LSAME
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC DCONJG
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( N.EQ.0 )
92 $ RETURN
93 *
94 IF( BETA.EQ.ZERO ) THEN
95 DO 20 J = 1, NRHS
96 DO 10 I = 1, N
97 B( I, J ) = ZERO
98 10 CONTINUE
99 20 CONTINUE
100 ELSE IF( BETA.EQ.-ONE ) THEN
101 DO 40 J = 1, NRHS
102 DO 30 I = 1, N
103 B( I, J ) = -B( I, J )
104 30 CONTINUE
105 40 CONTINUE
106 END IF
107 *
108 IF( ALPHA.EQ.ONE ) THEN
109 IF( LSAME( UPLO, 'U' ) ) THEN
110 *
111 * Compute B := B + A*X, where E is the superdiagonal of A.
112 *
113 DO 60 J = 1, NRHS
114 IF( N.EQ.1 ) THEN
115 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
116 ELSE
117 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
118 $ E( 1 )*X( 2, J )
119 B( N, J ) = B( N, J ) + DCONJG( E( N-1 ) )*
120 $ X( N-1, J ) + D( N )*X( N, J )
121 DO 50 I = 2, N - 1
122 B( I, J ) = B( I, J ) + DCONJG( E( I-1 ) )*
123 $ X( I-1, J ) + D( I )*X( I, J ) +
124 $ E( I )*X( I+1, J )
125 50 CONTINUE
126 END IF
127 60 CONTINUE
128 ELSE
129 *
130 * Compute B := B + A*X, where E is the subdiagonal of A.
131 *
132 DO 80 J = 1, NRHS
133 IF( N.EQ.1 ) THEN
134 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
135 ELSE
136 B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
137 $ DCONJG( E( 1 ) )*X( 2, J )
138 B( N, J ) = B( N, J ) + E( N-1 )*X( N-1, J ) +
139 $ D( N )*X( N, J )
140 DO 70 I = 2, N - 1
141 B( I, J ) = B( I, J ) + E( I-1 )*X( I-1, J ) +
142 $ D( I )*X( I, J ) +
143 $ DCONJG( E( I ) )*X( I+1, J )
144 70 CONTINUE
145 END IF
146 80 CONTINUE
147 END IF
148 ELSE IF( ALPHA.EQ.-ONE ) THEN
149 IF( LSAME( UPLO, 'U' ) ) THEN
150 *
151 * Compute B := B - A*X, where E is the superdiagonal of A.
152 *
153 DO 100 J = 1, NRHS
154 IF( N.EQ.1 ) THEN
155 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
156 ELSE
157 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
158 $ E( 1 )*X( 2, J )
159 B( N, J ) = B( N, J ) - DCONJG( E( N-1 ) )*
160 $ X( N-1, J ) - D( N )*X( N, J )
161 DO 90 I = 2, N - 1
162 B( I, J ) = B( I, J ) - DCONJG( E( I-1 ) )*
163 $ X( I-1, J ) - D( I )*X( I, J ) -
164 $ E( I )*X( I+1, J )
165 90 CONTINUE
166 END IF
167 100 CONTINUE
168 ELSE
169 *
170 * Compute B := B - A*X, where E is the subdiagonal of A.
171 *
172 DO 120 J = 1, NRHS
173 IF( N.EQ.1 ) THEN
174 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
175 ELSE
176 B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
177 $ DCONJG( E( 1 ) )*X( 2, J )
178 B( N, J ) = B( N, J ) - E( N-1 )*X( N-1, J ) -
179 $ D( N )*X( N, J )
180 DO 110 I = 2, N - 1
181 B( I, J ) = B( I, J ) - E( I-1 )*X( I-1, J ) -
182 $ D( I )*X( I, J ) -
183 $ DCONJG( E( I ) )*X( I+1, J )
184 110 CONTINUE
185 END IF
186 120 CONTINUE
187 END IF
188 END IF
189 RETURN
190 *
191 * End of ZLAPTM
192 *
193 END