![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpotrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zpotrf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZPOTRF( UPLO, N, A, LDA, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
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: * ZPOTRF computes the Cholesky factorization of a complex Hermitian
20: * positive definite matrix A.
21: *
22: * The factorization has the form
23: * A = U**H * U, if UPLO = 'U', or
24: * A = L * L**H, if UPLO = 'L',
25: * where U is an upper triangular matrix and L is lower triangular.
26: *
27: * This is the block version of the algorithm, calling Level 3 BLAS.
28: *
29: * Arguments
30: * =========
31: *
32: * UPLO (input) CHARACTER*1
33: * = 'U': Upper triangle of A is stored;
34: * = 'L': Lower triangle of A is stored.
35: *
36: * N (input) INTEGER
37: * The order of the matrix A. N >= 0.
38: *
39: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
40: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
41: * N-by-N upper triangular part of A contains the upper
42: * triangular part of the matrix A, and the strictly lower
43: * triangular part of A is not referenced. If UPLO = 'L', the
44: * leading N-by-N lower triangular part of A contains the lower
45: * triangular part of the matrix A, and the strictly upper
46: * triangular part of A is not referenced.
47: *
48: * On exit, if INFO = 0, the factor U or L from the Cholesky
49: * factorization A = U**H*U or A = L*L**H.
50: *
51: * LDA (input) INTEGER
52: * The leading dimension of the array A. LDA >= max(1,N).
53: *
54: * INFO (output) INTEGER
55: * = 0: successful exit
56: * < 0: if INFO = -i, the i-th argument had an illegal value
57: * > 0: if INFO = i, the leading minor of order i is not
58: * positive definite, and the factorization could not be
59: * completed.
60: *
61: * =====================================================================
62: *
63: * .. Parameters ..
64: DOUBLE PRECISION ONE
65: COMPLEX*16 CONE
66: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ) )
67: * ..
68: * .. Local Scalars ..
69: LOGICAL UPPER
70: INTEGER J, JB, NB
71: * ..
72: * .. External Functions ..
73: LOGICAL LSAME
74: INTEGER ILAENV
75: EXTERNAL LSAME, ILAENV
76: * ..
77: * .. External Subroutines ..
78: EXTERNAL XERBLA, ZGEMM, ZHERK, ZPOTF2, ZTRSM
79: * ..
80: * .. Intrinsic Functions ..
81: INTRINSIC MAX, MIN
82: * ..
83: * .. Executable Statements ..
84: *
85: * Test the input parameters.
86: *
87: INFO = 0
88: UPPER = LSAME( UPLO, 'U' )
89: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
90: INFO = -1
91: ELSE IF( N.LT.0 ) THEN
92: INFO = -2
93: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
94: INFO = -4
95: END IF
96: IF( INFO.NE.0 ) THEN
97: CALL XERBLA( 'ZPOTRF', -INFO )
98: RETURN
99: END IF
100: *
101: * Quick return if possible
102: *
103: IF( N.EQ.0 )
104: $ RETURN
105: *
106: * Determine the block size for this environment.
107: *
108: NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
109: IF( NB.LE.1 .OR. NB.GE.N ) THEN
110: *
111: * Use unblocked code.
112: *
113: CALL ZPOTF2( UPLO, N, A, LDA, INFO )
114: ELSE
115: *
116: * Use blocked code.
117: *
118: IF( UPPER ) THEN
119: *
120: * Compute the Cholesky factorization A = U'*U.
121: *
122: DO 10 J = 1, N, NB
123: *
124: * Update and factorize the current diagonal block and test
125: * for non-positive-definiteness.
126: *
127: JB = MIN( NB, N-J+1 )
128: CALL ZHERK( 'Upper', 'Conjugate transpose', JB, J-1,
129: $ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA )
130: CALL ZPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
131: IF( INFO.NE.0 )
132: $ GO TO 30
133: IF( J+JB.LE.N ) THEN
134: *
135: * Compute the current block row.
136: *
137: CALL ZGEMM( 'Conjugate transpose', 'No transpose', JB,
138: $ N-J-JB+1, J-1, -CONE, A( 1, J ), LDA,
139: $ A( 1, J+JB ), LDA, CONE, A( J, J+JB ),
140: $ LDA )
141: CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
142: $ 'Non-unit', JB, N-J-JB+1, CONE, A( J, J ),
143: $ LDA, A( J, J+JB ), LDA )
144: END IF
145: 10 CONTINUE
146: *
147: ELSE
148: *
149: * Compute the Cholesky factorization A = L*L'.
150: *
151: DO 20 J = 1, N, NB
152: *
153: * Update and factorize the current diagonal block and test
154: * for non-positive-definiteness.
155: *
156: JB = MIN( NB, N-J+1 )
157: CALL ZHERK( 'Lower', 'No transpose', JB, J-1, -ONE,
158: $ A( J, 1 ), LDA, ONE, A( J, J ), LDA )
159: CALL ZPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
160: IF( INFO.NE.0 )
161: $ GO TO 30
162: IF( J+JB.LE.N ) THEN
163: *
164: * Compute the current block column.
165: *
166: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
167: $ N-J-JB+1, JB, J-1, -CONE, A( J+JB, 1 ),
168: $ LDA, A( J, 1 ), LDA, CONE, A( J+JB, J ),
169: $ LDA )
170: CALL ZTRSM( 'Right', 'Lower', 'Conjugate transpose',
171: $ 'Non-unit', N-J-JB+1, JB, CONE, A( J, J ),
172: $ LDA, A( J+JB, J ), LDA )
173: END IF
174: 20 CONTINUE
175: END IF
176: END IF
177: GO TO 40
178: *
179: 30 CONTINUE
180: INFO = INFO + J - 1
181: *
182: 40 CONTINUE
183: RETURN
184: *
185: * End of ZPOTRF
186: *
187: END