|
5340
|
1 SUBROUTINE ZLAUUM( UPLO, N, A, LDA, INFO ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
5340
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 CHARACTER UPLO |
|
|
9 INTEGER INFO, LDA, N |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 COMPLEX*16 A( LDA, * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * ZLAUUM computes the product U * U' or L' * L, where the triangular |
|
|
19 * factor U or L is stored in the upper or lower triangular part of |
|
|
20 * the array A. |
|
|
21 * |
|
|
22 * If UPLO = 'U' or 'u' then the upper triangle of the result is stored, |
|
|
23 * overwriting the factor U in A. |
|
|
24 * If UPLO = 'L' or 'l' then the lower triangle of the result is stored, |
|
|
25 * overwriting the factor L in A. |
|
|
26 * |
|
|
27 * This is the blocked form of the algorithm, calling Level 3 BLAS. |
|
|
28 * |
|
|
29 * Arguments |
|
|
30 * ========= |
|
|
31 * |
|
|
32 * UPLO (input) CHARACTER*1 |
|
|
33 * Specifies whether the triangular factor stored in the array A |
|
|
34 * is upper or lower triangular: |
|
|
35 * = 'U': Upper triangular |
|
|
36 * = 'L': Lower triangular |
|
|
37 * |
|
|
38 * N (input) INTEGER |
|
|
39 * The order of the triangular factor U or L. N >= 0. |
|
|
40 * |
|
|
41 * A (input/output) COMPLEX*16 array, dimension (LDA,N) |
|
|
42 * On entry, the triangular factor U or L. |
|
|
43 * On exit, if UPLO = 'U', the upper triangle of A is |
|
|
44 * overwritten with the upper triangle of the product U * U'; |
|
|
45 * if UPLO = 'L', the lower triangle of A is overwritten with |
|
|
46 * the lower triangle of the product L' * L. |
|
|
47 * |
|
|
48 * LDA (input) INTEGER |
|
|
49 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
50 * |
|
|
51 * INFO (output) INTEGER |
|
|
52 * = 0: successful exit |
|
|
53 * < 0: if INFO = -k, the k-th argument had an illegal value |
|
|
54 * |
|
|
55 * ===================================================================== |
|
|
56 * |
|
|
57 * .. Parameters .. |
|
|
58 DOUBLE PRECISION ONE |
|
|
59 PARAMETER ( ONE = 1.0D+0 ) |
|
|
60 COMPLEX*16 CONE |
|
|
61 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) |
|
|
62 * .. |
|
|
63 * .. Local Scalars .. |
|
|
64 LOGICAL UPPER |
|
|
65 INTEGER I, IB, NB |
|
|
66 * .. |
|
|
67 * .. External Functions .. |
|
|
68 LOGICAL LSAME |
|
|
69 INTEGER ILAENV |
|
|
70 EXTERNAL LSAME, ILAENV |
|
|
71 * .. |
|
|
72 * .. External Subroutines .. |
|
|
73 EXTERNAL XERBLA, ZGEMM, ZHERK, ZLAUU2, ZTRMM |
|
|
74 * .. |
|
|
75 * .. Intrinsic Functions .. |
|
|
76 INTRINSIC MAX, MIN |
|
|
77 * .. |
|
|
78 * .. Executable Statements .. |
|
|
79 * |
|
|
80 * Test the input parameters. |
|
|
81 * |
|
|
82 INFO = 0 |
|
|
83 UPPER = LSAME( UPLO, 'U' ) |
|
|
84 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN |
|
|
85 INFO = -1 |
|
|
86 ELSE IF( N.LT.0 ) THEN |
|
|
87 INFO = -2 |
|
|
88 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
89 INFO = -4 |
|
|
90 END IF |
|
|
91 IF( INFO.NE.0 ) THEN |
|
|
92 CALL XERBLA( 'ZLAUUM', -INFO ) |
|
|
93 RETURN |
|
|
94 END IF |
|
|
95 * |
|
|
96 * Quick return if possible |
|
|
97 * |
|
|
98 IF( N.EQ.0 ) |
|
|
99 $ RETURN |
|
|
100 * |
|
|
101 * Determine the block size for this environment. |
|
|
102 * |
|
|
103 NB = ILAENV( 1, 'ZLAUUM', UPLO, N, -1, -1, -1 ) |
|
|
104 * |
|
|
105 IF( NB.LE.1 .OR. NB.GE.N ) THEN |
|
|
106 * |
|
|
107 * Use unblocked code |
|
|
108 * |
|
|
109 CALL ZLAUU2( UPLO, N, A, LDA, INFO ) |
|
|
110 ELSE |
|
|
111 * |
|
|
112 * Use blocked code |
|
|
113 * |
|
|
114 IF( UPPER ) THEN |
|
|
115 * |
|
|
116 * Compute the product U * U'. |
|
|
117 * |
|
|
118 DO 10 I = 1, N, NB |
|
|
119 IB = MIN( NB, N-I+1 ) |
|
|
120 CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', |
|
|
121 $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, |
|
|
122 $ A( 1, I ), LDA ) |
|
|
123 CALL ZLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) |
|
|
124 IF( I+IB.LE.N ) THEN |
|
|
125 CALL ZGEMM( 'No transpose', 'Conjugate transpose', |
|
|
126 $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), |
|
|
127 $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), |
|
|
128 $ LDA ) |
|
|
129 CALL ZHERK( 'Upper', 'No transpose', IB, N-I-IB+1, |
|
|
130 $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), |
|
|
131 $ LDA ) |
|
|
132 END IF |
|
|
133 10 CONTINUE |
|
|
134 ELSE |
|
|
135 * |
|
|
136 * Compute the product L' * L. |
|
|
137 * |
|
|
138 DO 20 I = 1, N, NB |
|
|
139 IB = MIN( NB, N-I+1 ) |
|
|
140 CALL ZTRMM( 'Left', 'Lower', 'Conjugate transpose', |
|
|
141 $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, |
|
|
142 $ A( I, 1 ), LDA ) |
|
|
143 CALL ZLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) |
|
|
144 IF( I+IB.LE.N ) THEN |
|
|
145 CALL ZGEMM( 'Conjugate transpose', 'No transpose', IB, |
|
|
146 $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, |
|
|
147 $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) |
|
|
148 CALL ZHERK( 'Lower', 'Conjugate transpose', IB, |
|
|
149 $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, |
|
|
150 $ A( I, I ), LDA ) |
|
|
151 END IF |
|
|
152 20 CONTINUE |
|
|
153 END IF |
|
|
154 END IF |
|
|
155 * |
|
|
156 RETURN |
|
|
157 * |
|
|
158 * End of ZLAUUM |
|
|
159 * |
|
|
160 END |
