LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clahef_aa.f
Go to the documentation of this file.
1*> \brief \b CLAHEF_AA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CLAHEF_AA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahef_aa.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahef_aa.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahef_aa.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
22* H, LDH, WORK )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER J1, M, NB, LDA, LDH
27* ..
28* .. Array Arguments ..
29* INTEGER IPIV( * )
30* COMPLEX A( LDA, * ), H( LDH, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CLAHEF_AA factorizes a panel of a complex hermitian matrix A using
40*> the Aasen's algorithm. The panel consists of a set of NB rows of A
41*> when UPLO is U, or a set of NB columns when UPLO is L.
42*>
43*> In order to factorize the panel, the Aasen's algorithm requires the
44*> last row, or column, of the previous panel. The first row, or column,
45*> of A is set to be the first row, or column, of an identity matrix,
46*> which is used to factorize the first panel.
47*>
48*> The resulting J-th row of U, or J-th column of L, is stored in the
49*> (J-1)-th row, or column, of A (without the unit diagonals), while
50*> the diagonal and subdiagonal of A are overwritten by those of T.
51*>
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] UPLO
58*> \verbatim
59*> UPLO is CHARACTER*1
60*> = 'U': Upper triangle of A is stored;
61*> = 'L': Lower triangle of A is stored.
62*> \endverbatim
63*>
64*> \param[in] J1
65*> \verbatim
66*> J1 is INTEGER
67*> The location of the first row, or column, of the panel
68*> within the submatrix of A, passed to this routine, e.g.,
69*> when called by CHETRF_AA, for the first panel, J1 is 1,
70*> while for the remaining panels, J1 is 2.
71*> \endverbatim
72*>
73*> \param[in] M
74*> \verbatim
75*> M is INTEGER
76*> The dimension of the submatrix. M >= 0.
77*> \endverbatim
78*>
79*> \param[in] NB
80*> \verbatim
81*> NB is INTEGER
82*> The dimension of the panel to be facotorized.
83*> \endverbatim
84*>
85*> \param[in,out] A
86*> \verbatim
87*> A is COMPLEX array, dimension (LDA,M) for
88*> the first panel, while dimension (LDA,M+1) for the
89*> remaining panels.
90*>
91*> On entry, A contains the last row, or column, of
92*> the previous panel, and the trailing submatrix of A
93*> to be factorized, except for the first panel, only
94*> the panel is passed.
95*>
96*> On exit, the leading panel is factorized.
97*> \endverbatim
98*>
99*> \param[in] LDA
100*> \verbatim
101*> LDA is INTEGER
102*> The leading dimension of the array A. LDA >= max(1,N).
103*> \endverbatim
104*>
105*> \param[out] IPIV
106*> \verbatim
107*> IPIV is INTEGER array, dimension (N)
108*> Details of the row and column interchanges,
109*> the row and column k were interchanged with the row and
110*> column IPIV(k).
111*> \endverbatim
112*>
113*> \param[in,out] H
114*> \verbatim
115*> H is COMPLEX workspace, dimension (LDH,NB).
116*>
117*> \endverbatim
118*>
119*> \param[in] LDH
120*> \verbatim
121*> LDH is INTEGER
122*> The leading dimension of the workspace H. LDH >= max(1,M).
123*> \endverbatim
124*>
125*> \param[out] WORK
126*> \verbatim
127*> WORK is COMPLEX workspace, dimension (M).
128*> \endverbatim
129*>
130*
131* Authors:
132* ========
133*
134*> \author Univ. of Tennessee
135*> \author Univ. of California Berkeley
136*> \author Univ. of Colorado Denver
137*> \author NAG Ltd.
138*
139*> \ingroup lahef_aa
140*
141* =====================================================================
142 SUBROUTINE clahef_aa( UPLO, J1, M, NB, A, LDA, IPIV,
143 $ H, LDH, WORK )
144*
145* -- LAPACK computational routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149 IMPLICIT NONE
150*
151* .. Scalar Arguments ..
152 CHARACTER UPLO
153 INTEGER M, NB, J1, LDA, LDH
154* ..
155* .. Array Arguments ..
156 INTEGER IPIV( * )
157 COMPLEX A( LDA, * ), H( LDH, * ), WORK( * )
158* ..
159*
160* =====================================================================
161* .. Parameters ..
162 COMPLEX ZERO, ONE
163 parameter( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
164*
165* .. Local Scalars ..
166 INTEGER J, K, K1, I1, I2, MJ
167 COMPLEX PIV, ALPHA
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER ICAMAX, ILAENV
172 EXTERNAL lsame, ilaenv, icamax
173* ..
174* .. External Subroutines ..
175 EXTERNAL clacgv, cgemv, cscal, caxpy, ccopy, cswap, claset,
176 $ xerbla
177* ..
178* .. Intrinsic Functions ..
179 INTRINSIC real, conjg, max
180* ..
181* .. Executable Statements ..
182*
183 j = 1
184*
185* K1 is the first column of the panel to be factorized
186* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
187*
188 k1 = (2-j1)+1
189*
190 IF( lsame( uplo, 'U' ) ) THEN
191*
192* .....................................................
193* Factorize A as U**T*D*U using the upper triangle of A
194* .....................................................
195*
196 10 CONTINUE
197 IF ( j.GT.min(m, nb) )
198 $ GO TO 20
199*
200* K is the column to be factorized
201* when being called from CHETRF_AA,
202* > for the first block column, J1 is 1, hence J1+J-1 is J,
203* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
204*
205 k = j1+j-1
206 IF( j.EQ.m ) THEN
207*
208* Only need to compute T(J, J)
209*
210 mj = 1
211 ELSE
212 mj = m-j+1
213 END IF
214*
215* H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
216* where H(J:N, J) has been initialized to be A(J, J:N)
217*
218 IF( k.GT.2 ) THEN
219*
220* K is the column to be factorized
221* > for the first block column, K is J, skipping the first two
222* columns
223* > for the rest of the columns, K is J+1, skipping only the
224* first column
225*
226 CALL clacgv( j-k1, a( 1, j ), 1 )
227 CALL cgemv( 'No transpose', mj, j-k1,
228 $ -one, h( j, k1 ), ldh,
229 $ a( 1, j ), 1,
230 $ one, h( j, j ), 1 )
231 CALL clacgv( j-k1, a( 1, j ), 1 )
232 END IF
233*
234* Copy H(i:n, i) into WORK
235*
236 CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
237*
238 IF( j.GT.k1 ) THEN
239*
240* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
241* where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
242*
243 alpha = -conjg( a( k-1, j ) )
244 CALL caxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
245 END IF
246*
247* Set A(J, J) = T(J, J)
248*
249 a( k, j ) = real( work( 1 ) )
250*
251 IF( j.LT.m ) THEN
252*
253* Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
254* where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
255*
256 IF( k.GT.1 ) THEN
257 alpha = -a( k, j )
258 CALL caxpy( m-j, alpha, a( k-1, j+1 ), lda,
259 $ work( 2 ), 1 )
260 ENDIF
261*
262* Find max(|WORK(2:n)|)
263*
264 i2 = icamax( m-j, work( 2 ), 1 ) + 1
265 piv = work( i2 )
266*
267* Apply hermitian pivot
268*
269 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
270*
271* Swap WORK(I1) and WORK(I2)
272*
273 i1 = 2
274 work( i2 ) = work( i1 )
275 work( i1 ) = piv
276*
277* Swap A(I1, I1+1:N) with A(I1+1:N, I2)
278*
279 i1 = i1+j-1
280 i2 = i2+j-1
281 CALL cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
282 $ a( j1+i1, i2 ), 1 )
283 CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
284 CALL clacgv( i2-i1-1, a( j1+i1, i2 ), 1 )
285*
286* Swap A(I1, I2+1:N) with A(I2, I2+1:N)
287*
288 IF( i2.LT.m )
289 $ CALL cswap( m-i2, a( j1+i1-1, i2+1 ), lda,
290 $ a( j1+i2-1, i2+1 ), lda )
291*
292* Swap A(I1, I1) with A(I2,I2)
293*
294 piv = a( i1+j1-1, i1 )
295 a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
296 a( j1+i2-1, i2 ) = piv
297*
298* Swap H(I1, 1:J1) with H(I2, 1:J1)
299*
300 CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
301 ipiv( i1 ) = i2
302*
303 IF( i1.GT.(k1-1) ) THEN
304*
305* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
306* skipping the first column
307*
308 CALL cswap( i1-k1+1, a( 1, i1 ), 1,
309 $ a( 1, i2 ), 1 )
310 END IF
311 ELSE
312 ipiv( j+1 ) = j+1
313 ENDIF
314*
315* Set A(J, J+1) = T(J, J+1)
316*
317 a( k, j+1 ) = work( 2 )
318*
319 IF( j.LT.nb ) THEN
320*
321* Copy A(J+1:N, J+1) into H(J:N, J),
322*
323 CALL ccopy( m-j, a( k+1, j+1 ), lda,
324 $ h( j+1, j+1 ), 1 )
325 END IF
326*
327* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
328* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
329*
330 IF( j.LT.(m-1) ) THEN
331 IF( a( k, j+1 ).NE.zero ) THEN
332 alpha = one / a( k, j+1 )
333 CALL ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
334 CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
335 ELSE
336 CALL claset( 'Full', 1, m-j-1, zero, zero,
337 $ a( k, j+2 ), lda)
338 END IF
339 END IF
340 END IF
341 j = j + 1
342 GO TO 10
343 20 CONTINUE
344*
345 ELSE
346*
347* .....................................................
348* Factorize A as L*D*L**T using the lower triangle of A
349* .....................................................
350*
351 30 CONTINUE
352 IF( j.GT.min( m, nb ) )
353 $ GO TO 40
354*
355* K is the column to be factorized
356* when being called from CHETRF_AA,
357* > for the first block column, J1 is 1, hence J1+J-1 is J,
358* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
359*
360 k = j1+j-1
361 IF( j.EQ.m ) THEN
362*
363* Only need to compute T(J, J)
364*
365 mj = 1
366 ELSE
367 mj = m-j+1
368 END IF
369*
370* H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
371* where H(J:N, J) has been initialized to be A(J:N, J)
372*
373 IF( k.GT.2 ) THEN
374*
375* K is the column to be factorized
376* > for the first block column, K is J, skipping the first two
377* columns
378* > for the rest of the columns, K is J+1, skipping only the
379* first column
380*
381 CALL clacgv( j-k1, a( j, 1 ), lda )
382 CALL cgemv( 'No transpose', mj, j-k1,
383 $ -one, h( j, k1 ), ldh,
384 $ a( j, 1 ), lda,
385 $ one, h( j, j ), 1 )
386 CALL clacgv( j-k1, a( j, 1 ), lda )
387 END IF
388*
389* Copy H(J:N, J) into WORK
390*
391 CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
392*
393 IF( j.GT.k1 ) THEN
394*
395* Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
396* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
397*
398 alpha = -conjg( a( j, k-1 ) )
399 CALL caxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
400 END IF
401*
402* Set A(J, J) = T(J, J)
403*
404 a( j, k ) = real( work( 1 ) )
405*
406 IF( j.LT.m ) THEN
407*
408* Compute WORK(2:N) = T(J, J) L((J+1):N, J)
409* where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
410*
411 IF( k.GT.1 ) THEN
412 alpha = -a( j, k )
413 CALL caxpy( m-j, alpha, a( j+1, k-1 ), 1,
414 $ work( 2 ), 1 )
415 ENDIF
416*
417* Find max(|WORK(2:n)|)
418*
419 i2 = icamax( m-j, work( 2 ), 1 ) + 1
420 piv = work( i2 )
421*
422* Apply hermitian pivot
423*
424 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
425*
426* Swap WORK(I1) and WORK(I2)
427*
428 i1 = 2
429 work( i2 ) = work( i1 )
430 work( i1 ) = piv
431*
432* Swap A(I1+1:N, I1) with A(I2, I1+1:N)
433*
434 i1 = i1+j-1
435 i2 = i2+j-1
436 CALL cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
437 $ a( i2, j1+i1 ), lda )
438 CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
439 CALL clacgv( i2-i1-1, a( i2, j1+i1 ), lda )
440*
441* Swap A(I2+1:N, I1) with A(I2+1:N, I2)
442*
443 IF( i2.LT.m )
444 $ CALL cswap( m-i2, a( i2+1, j1+i1-1 ), 1,
445 $ a( i2+1, j1+i2-1 ), 1 )
446*
447* Swap A(I1, I1) with A(I2, I2)
448*
449 piv = a( i1, j1+i1-1 )
450 a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
451 a( i2, j1+i2-1 ) = piv
452*
453* Swap H(I1, I1:J1) with H(I2, I2:J1)
454*
455 CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
456 ipiv( i1 ) = i2
457*
458 IF( i1.GT.(k1-1) ) THEN
459*
460* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
461* skipping the first column
462*
463 CALL cswap( i1-k1+1, a( i1, 1 ), lda,
464 $ a( i2, 1 ), lda )
465 END IF
466 ELSE
467 ipiv( j+1 ) = j+1
468 ENDIF
469*
470* Set A(J+1, J) = T(J+1, J)
471*
472 a( j+1, k ) = work( 2 )
473*
474 IF( j.LT.nb ) THEN
475*
476* Copy A(J+1:N, J+1) into H(J+1:N, J),
477*
478 CALL ccopy( m-j, a( j+1, k+1 ), 1,
479 $ h( j+1, j+1 ), 1 )
480 END IF
481*
482* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
483* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
484*
485 IF( j.LT.(m-1) ) THEN
486 IF( a( j+1, k ).NE.zero ) THEN
487 alpha = one / a( j+1, k )
488 CALL ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
489 CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
490 ELSE
491 CALL claset( 'Full', m-j-1, 1, zero, zero,
492 $ a( j+2, k ), lda )
493 END IF
494 END IF
495 END IF
496 j = j + 1
497 GO TO 30
498 40 CONTINUE
499 END IF
500 RETURN
501*
502* End of CLAHEF_AA
503*
504 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
subroutine clahef_aa(uplo, j1, m, nb, a, lda, ipiv, h, ldh, work)
CLAHEF_AA
Definition clahef_aa.f:144
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81