Skip to content

Commit 718b92b

Browse files
iyamazakiiyamazaki
committed
two-stage Aasens in C, D, and Z precisions
1 parent 813984b commit 718b92b

9 files changed

Lines changed: 3647 additions & 0 deletions

‎SRC/chesv_aa_2stage.f‎

Lines changed: 285 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,285 @@
1+
*> \brief <b> CHESV_AA_2STAGE computes the solution to system of linear equations A * X = B for HE matrices</b>
2+
*
3+
* @precisions fortran c -> s d
4+
*
5+
* =========== DOCUMENTATION ===========
6+
*
7+
* Online html documentation available at
8+
* http://www.netlib.org/lapack/explore-html/
9+
*
10+
*> \htmlonly
11+
*> Download CHESV_AA_2STAGE + dependencies
12+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chesv_aa_2stage.f">
13+
*> [TGZ]</a>
14+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chesv_aa_2stage.f">
15+
*> [ZIP]</a>
16+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chesv_aa_2stage.f">
17+
*> [TXT]</a>
18+
*> \endhtmlonly
19+
*
20+
* Definition:
21+
* ===========
22+
*
23+
* SUBROUTINE CHESV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
24+
* IPIV, IPIV2, B, LDB, WORK, LWORK,
25+
* INFO )
26+
*
27+
* .. Scalar Arguments ..
28+
* CHARACTER UPLO
29+
* INTEGER N, NRHS, LDA, LTB, LDB, LWORK, INFO
30+
* ..
31+
* .. Array Arguments ..
32+
* INTEGER IPIV( * ), IPIV2( * )
33+
* COMPLEX A( LDA, * ), TB( * ), B( LDB, *), WORK( * )
34+
* ..
35+
*
36+
*> \par Purpose:
37+
* =============
38+
*>
39+
*> \verbatim
40+
*>
41+
*> CHESV_AA_2STAGE computes the solution to a complex system of
42+
*> linear equations
43+
*> A * X = B,
44+
*> where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS
45+
*> matrices.
46+
*>
47+
*> Aasen's 2-stage algorithm is used to factor A as
48+
*> A = U * T * U**H, if UPLO = 'U', or
49+
*> A = L * T * L**H, if UPLO = 'L',
50+
*> where U (or L) is a product of permutation and unit upper (lower)
51+
*> triangular matrices, and T is Hermitian and band. The matrix T is
52+
*> then LU-factored with partial pivoting. The factored form of A
53+
*> is then used to solve the system of equations A * X = B.
54+
*>
55+
*> This is the blocked version of the algorithm, calling Level 3 BLAS.
56+
*> \endverbatim
57+
*
58+
* Arguments:
59+
* ==========
60+
*
61+
*> \param[in] UPLO
62+
*> \verbatim
63+
*> UPLO is CHARACTER*1
64+
*> = 'U': Upper triangle of A is stored;
65+
*> = 'L': Lower triangle of A is stored.
66+
*> \endverbatim
67+
*>
68+
*> \param[in] N
69+
*> \verbatim
70+
*> N is INTEGER
71+
*> The order of the matrix A. N >= 0.
72+
*> \endverbatim
73+
*>
74+
*> \param[in] NRHS
75+
*> \verbatim
76+
*> NRHS is INTEGER
77+
*> The number of right hand sides, i.e., the number of columns
78+
*> of the matrix B. NRHS >= 0.
79+
*> \endverbatim
80+
*>
81+
*> \param[in,out] A
82+
*> \verbatim
83+
*> A is COMPLEX array, dimension (LDA,N)
84+
*> On entry, the hermitian matrix A. If UPLO = 'U', the leading
85+
*> N-by-N upper triangular part of A contains the upper
86+
*> triangular part of the matrix A, and the strictly lower
87+
*> triangular part of A is not referenced. If UPLO = 'L', the
88+
*> leading N-by-N lower triangular part of A contains the lower
89+
*> triangular part of the matrix A, and the strictly upper
90+
*> triangular part of A is not referenced.
91+
*>
92+
*> On exit, L is stored below (or above) the subdiaonal blocks,
93+
*> when UPLO is 'L' (or 'U').
94+
*> \endverbatim
95+
*>
96+
*> \param[in] LDA
97+
*> \verbatim
98+
*> LDA is INTEGER
99+
*> The leading dimension of the array A. LDA >= max(1,N).
100+
*> \endverbatim
101+
*>
102+
*> \param[out] TB
103+
*> \verbatim
104+
*> TB is COMPLEX array, dimension (LTB)
105+
*> On exit, details of the LU factorization of the band matrix.
106+
*> \endverbatim
107+
*>
108+
*> \param[in] LTB
109+
*> \verbatim
110+
*> The size of the array TB. LTB >= 4*N, internally
111+
*> used to select NB such that LTB >= (3*NB+1)*N.
112+
*>
113+
*> If LTB = -1, then a workspace query is assumed; the
114+
*> routine only calculates the optimal size of LTB,
115+
*> returns this value as the first entry of TB, and
116+
*> no error message related to LTB is issued by XERBLA.
117+
*> \endverbatim
118+
*>
119+
*> \param[out] IPIV
120+
*> \verbatim
121+
*> IPIV is INTEGER array, dimension (N)
122+
*> On exit, it contains the details of the interchanges, i.e.,
123+
*> the row and column k of A were interchanged with the
124+
*> row and column IPIV(k).
125+
*> \endverbatim
126+
*>
127+
*> \param[out] IPIV2
128+
*> \verbatim
129+
*> IPIV is INTEGER array, dimension (N)
130+
*> On exit, it contains the details of the interchanges, i.e.,
131+
*> the row and column k of T were interchanged with the
132+
*> row and column IPIV(k).
133+
*> \endverbatim
134+
*>
135+
*> \param[in,out] B
136+
*> \verbatim
137+
*> B is COMPLEX array, dimension (LDB,NRHS)
138+
*> On entry, the right hand side matrix B.
139+
*> On exit, the solution matrix X.
140+
*> \endverbatim
141+
*>
142+
*> \param[in] LDB
143+
*> \verbatim
144+
*> LDB is INTEGER
145+
*> The leading dimension of the array B. LDB >= max(1,N).
146+
*> \endverbatim
147+
*>
148+
*> \param[out] WORK
149+
*> \verbatim
150+
*> WORK is COMPLEX workspace of size LWORK
151+
*> \endverbatim
152+
*>
153+
*> \param[in] LWORK
154+
*> \verbatim
155+
*> The size of WORK. LWORK >= N, internally used to select NB
156+
*> such that LWORK >= N*NB.
157+
*>
158+
*> If LWORK = -1, then a workspace query is assumed; the
159+
*> routine only calculates the optimal size of the WORK array,
160+
*> returns this value as the first entry of the WORK array, and
161+
*> no error message related to LWORK is issued by XERBLA.
162+
*> \endverbatim
163+
*>
164+
*> \param[out] INFO
165+
*> \verbatim
166+
*> INFO is INTEGER
167+
*> = 0: successful exit
168+
*> < 0: if INFO = -i, the i-th argument had an illegal value.
169+
*> > 0: if INFO = i, band LU factorization failed on i-th column
170+
*> \endverbatim
171+
*
172+
* Authors:
173+
* ========
174+
*
175+
*> \author Univ. of Tennessee
176+
*> \author Univ. of California Berkeley
177+
*> \author Univ. of Colorado Denver
178+
*> \author NAG Ltd.
179+
*
180+
*> \date December 2016
181+
*
182+
*> \ingroup complexSYcomputational
183+
*
184+
* =====================================================================
185+
SUBROUTINE CHESV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
186+
$ IPIV, IPIV2, B, LDB, WORK, LWORK,
187+
$ INFO )
188+
*
189+
* -- LAPACK computational routine (version 3.7.0) --
190+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
191+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192+
* December 2016
193+
*
194+
IMPLICIT NONE
195+
*
196+
* .. Scalar Arguments ..
197+
CHARACTER UPLO
198+
INTEGER N, NRHS, LDA, LDB, LTB, LWORK, INFO
199+
* ..
200+
* .. Array Arguments ..
201+
INTEGER IPIV( * ), IPIV2( * )
202+
COMPLEX A( LDA, * ), B( LDB, * ), TB( * ), WORK( * )
203+
* ..
204+
*
205+
* =====================================================================
206+
* .. Parameters ..
207+
COMPLEX ZERO, ONE
208+
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
209+
$ ONE = ( 1.0E+0, 0.0E+0 ) )
210+
*
211+
* .. Local Scalars ..
212+
LOGICAL UPPER, TQUERY, WQUERY
213+
INTEGER I, J, K, I1, I2, TD
214+
INTEGER LDTB, LWKOPT, NB, KB, NT, IINFO
215+
COMPLEX PIV
216+
* ..
217+
* .. External Functions ..
218+
LOGICAL LSAME
219+
INTEGER ILAENV
220+
EXTERNAL LSAME, ILAENV
221+
* ..
222+
* .. External Subroutines ..
223+
EXTERNAL XERBLA
224+
* ..
225+
* .. Intrinsic Functions ..
226+
INTRINSIC CONJG, MIN, MAX
227+
* ..
228+
* .. Executable Statements ..
229+
*
230+
* Test the input parameters.
231+
*
232+
INFO = 0
233+
UPPER = LSAME( UPLO, 'U' )
234+
WQUERY = ( LWORK.EQ.-1 )
235+
TQUERY = ( LTB.EQ.-1 )
236+
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
237+
INFO = -1
238+
ELSE IF( N.LT.0 ) THEN
239+
INFO = -2
240+
ELSE IF( NRHS.LT.0 ) THEN
241+
INFO = -3
242+
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
243+
INFO = -5
244+
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
245+
INFO = -11
246+
END IF
247+
*
248+
IF( INFO.EQ.0 ) THEN
249+
CALL CHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, -1, IPIV,
250+
$ IPIV2, WORK, -1, INFO )
251+
LWKOPT = INT( WORK(1) )
252+
IF( LWORK.LT.LWKOPT .AND. .NOT.WQUERY ) THEN
253+
INFO = -13
254+
END IF
255+
IF( LTB.LT.INT( TB(1) ) .AND. .NOT.TQUERY ) THEN
256+
INFO = -7
257+
END IF
258+
END IF
259+
*
260+
IF( INFO.NE.0 ) THEN
261+
CALL XERBLA( 'CHESV_AA_2STAGE', -INFO )
262+
RETURN
263+
ELSE IF( WQUERY .OR. TQUERY ) THEN
264+
RETURN
265+
END IF
266+
*
267+
*
268+
* Compute the factorization A = U*T*U**H or A = L*T*L**H.
269+
*
270+
CALL CHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
271+
$ WORK, LWORK, INFO )
272+
IF( INFO.EQ.0 ) THEN
273+
*
274+
* Solve the system A*X = B, overwriting B with X.
275+
*
276+
CALL CHETRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB, IPIV,
277+
$ IPIV2, B, LDB, INFO )
278+
*
279+
END IF
280+
*
281+
WORK( 1 ) = LWKOPT
282+
*
283+
* End of CHESV_AA_2STAGE
284+
*
285+
END

0 commit comments

Comments
 (0)