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