LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 complexSYcomputational
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)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
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 clahef_aa(UPLO, J1, M, NB, A, LDA, IPIV, H, LDH, WORK)
CLAHEF_AA
Definition: clahef_aa.f:144