LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zpftrf.f
Go to the documentation of this file.
1 *> \brief \b ZPFTRF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZPFTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpftrf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpftrf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpftrf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZPFTRF( TRANSR, UPLO, N, A, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER TRANSR, UPLO
25 * INTEGER N, INFO
26 * ..
27 * .. Array Arguments ..
28 * COMPLEX*16 A( 0: * )
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> ZPFTRF computes the Cholesky factorization of a complex Hermitian
37 *> positive definite matrix A.
38 *>
39 *> The factorization has the form
40 *> A = U**H * U, if UPLO = 'U', or
41 *> A = L * L**H, if UPLO = 'L',
42 *> where U is an upper triangular matrix and L is lower triangular.
43 *>
44 *> This is the block version of the algorithm, calling Level 3 BLAS.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] TRANSR
51 *> \verbatim
52 *> TRANSR is CHARACTER*1
53 *> = 'N': The Normal TRANSR of RFP A is stored;
54 *> = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
55 *> \endverbatim
56 *>
57 *> \param[in] UPLO
58 *> \verbatim
59 *> UPLO is CHARACTER*1
60 *> = 'U': Upper triangle of RFP A is stored;
61 *> = 'L': Lower triangle of RFP A is stored.
62 *> \endverbatim
63 *>
64 *> \param[in] N
65 *> \verbatim
66 *> N is INTEGER
67 *> The order of the matrix A. N >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in,out] A
71 *> \verbatim
72 *> A is COMPLEX*16 array, dimension ( N*(N+1)/2 );
73 *> On entry, the Hermitian matrix A in RFP format. RFP format is
74 *> described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
75 *> then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
76 *> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
77 *> the Conjugate-transpose of RFP A as defined when
78 *> TRANSR = 'N'. The contents of RFP A are defined by UPLO as
79 *> follows: If UPLO = 'U' the RFP A contains the nt elements of
80 *> upper packed A. If UPLO = 'L' the RFP A contains the elements
81 *> of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
82 *> 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
83 *> is odd. See the Note below for more details.
84 *>
85 *> On exit, if INFO = 0, the factor U or L from the Cholesky
86 *> factorization RFP A = U**H*U or RFP A = L*L**H.
87 *> \endverbatim
88 *>
89 *> \param[out] INFO
90 *> \verbatim
91 *> INFO is INTEGER
92 *> = 0: successful exit
93 *> < 0: if INFO = -i, the i-th argument had an illegal value
94 *> > 0: if INFO = i, the leading minor of order i is not
95 *> positive definite, and the factorization could not be
96 *> completed.
97 *>
98 *> Further Notes on RFP Format:
99 *> ============================
100 *>
101 *> We first consider Standard Packed Format when N is even.
102 *> We give an example where N = 6.
103 *>
104 *> AP is Upper AP is Lower
105 *>
106 *> 00 01 02 03 04 05 00
107 *> 11 12 13 14 15 10 11
108 *> 22 23 24 25 20 21 22
109 *> 33 34 35 30 31 32 33
110 *> 44 45 40 41 42 43 44
111 *> 55 50 51 52 53 54 55
112 *>
113 *> Let TRANSR = 'N'. RFP holds AP as follows:
114 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
115 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
116 *> conjugate-transpose of the first three columns of AP upper.
117 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
118 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
119 *> conjugate-transpose of the last three columns of AP lower.
120 *> To denote conjugate we place -- above the element. This covers the
121 *> case N even and TRANSR = 'N'.
122 *>
123 *> RFP A RFP A
124 *>
125 *> -- -- --
126 *> 03 04 05 33 43 53
127 *> -- --
128 *> 13 14 15 00 44 54
129 *> --
130 *> 23 24 25 10 11 55
131 *>
132 *> 33 34 35 20 21 22
133 *> --
134 *> 00 44 45 30 31 32
135 *> -- --
136 *> 01 11 55 40 41 42
137 *> -- -- --
138 *> 02 12 22 50 51 52
139 *>
140 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
141 *> transpose of RFP A above. One therefore gets:
142 *>
143 *> RFP A RFP A
144 *>
145 *> -- -- -- -- -- -- -- -- -- --
146 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
147 *> -- -- -- -- -- -- -- -- -- --
148 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
149 *> -- -- -- -- -- -- -- -- -- --
150 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
151 *>
152 *> We next consider Standard Packed Format when N is odd.
153 *> We give an example where N = 5.
154 *>
155 *> AP is Upper AP is Lower
156 *>
157 *> 00 01 02 03 04 00
158 *> 11 12 13 14 10 11
159 *> 22 23 24 20 21 22
160 *> 33 34 30 31 32 33
161 *> 44 40 41 42 43 44
162 *>
163 *> Let TRANSR = 'N'. RFP holds AP as follows:
164 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
165 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
166 *> conjugate-transpose of the first two columns of AP upper.
167 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
168 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
169 *> conjugate-transpose of the last two columns of AP lower.
170 *> To denote conjugate we place -- above the element. This covers the
171 *> case N odd and TRANSR = 'N'.
172 *>
173 *> RFP A RFP A
174 *>
175 *> -- --
176 *> 02 03 04 00 33 43
177 *> --
178 *> 12 13 14 10 11 44
179 *>
180 *> 22 23 24 20 21 22
181 *> --
182 *> 00 33 34 30 31 32
183 *> -- --
184 *> 01 11 44 40 41 42
185 *>
186 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
187 *> transpose of RFP A above. One therefore gets:
188 *>
189 *> RFP A RFP A
190 *>
191 *> -- -- -- -- -- -- -- -- --
192 *> 02 12 22 00 01 00 10 20 30 40 50
193 *> -- -- -- -- -- -- -- -- --
194 *> 03 13 23 33 11 33 11 21 31 41 51
195 *> -- -- -- -- -- -- -- -- --
196 *> 04 14 24 34 44 43 44 22 32 42 52
197 *> \endverbatim
198 *
199 * Authors:
200 * ========
201 *
202 *> \author Univ. of Tennessee
203 *> \author Univ. of California Berkeley
204 *> \author Univ. of Colorado Denver
205 *> \author NAG Ltd.
206 *
207 *> \ingroup complex16OTHERcomputational
208 *
209 * =====================================================================
210  SUBROUTINE zpftrf( TRANSR, UPLO, N, A, INFO )
211 *
212 * -- LAPACK computational routine --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 *
216 * .. Scalar Arguments ..
217  CHARACTER TRANSR, UPLO
218  INTEGER N, INFO
219 * ..
220 * .. Array Arguments ..
221  COMPLEX*16 A( 0: * )
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  DOUBLE PRECISION ONE
227  COMPLEX*16 CONE
228  parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ) )
229 * ..
230 * .. Local Scalars ..
231  LOGICAL LOWER, NISODD, NORMALTRANSR
232  INTEGER N1, N2, K
233 * ..
234 * .. External Functions ..
235  LOGICAL LSAME
236  EXTERNAL lsame
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL xerbla, zherk, zpotrf, ztrsm
240 * ..
241 * .. Intrinsic Functions ..
242  INTRINSIC mod
243 * ..
244 * .. Executable Statements ..
245 *
246 * Test the input parameters.
247 *
248  info = 0
249  normaltransr = lsame( transr, 'N' )
250  lower = lsame( uplo, 'L' )
251  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
252  info = -1
253  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
254  info = -2
255  ELSE IF( n.LT.0 ) THEN
256  info = -3
257  END IF
258  IF( info.NE.0 ) THEN
259  CALL xerbla( 'ZPFTRF', -info )
260  RETURN
261  END IF
262 *
263 * Quick return if possible
264 *
265  IF( n.EQ.0 )
266  $ RETURN
267 *
268 * If N is odd, set NISODD = .TRUE.
269 * If N is even, set K = N/2 and NISODD = .FALSE.
270 *
271  IF( mod( n, 2 ).EQ.0 ) THEN
272  k = n / 2
273  nisodd = .false.
274  ELSE
275  nisodd = .true.
276  END IF
277 *
278 * Set N1 and N2 depending on LOWER
279 *
280  IF( lower ) THEN
281  n2 = n / 2
282  n1 = n - n2
283  ELSE
284  n1 = n / 2
285  n2 = n - n1
286  END IF
287 *
288 * start execution: there are eight cases
289 *
290  IF( nisodd ) THEN
291 *
292 * N is odd
293 *
294  IF( normaltransr ) THEN
295 *
296 * N is odd and TRANSR = 'N'
297 *
298  IF( lower ) THEN
299 *
300 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
301 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
302 * T1 -> a(0), T2 -> a(n), S -> a(n1)
303 *
304  CALL zpotrf( 'L', n1, a( 0 ), n, info )
305  IF( info.GT.0 )
306  $ RETURN
307  CALL ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ), n,
308  $ a( n1 ), n )
309  CALL zherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
310  $ a( n ), n )
311  CALL zpotrf( 'U', n2, a( n ), n, info )
312  IF( info.GT.0 )
313  $ info = info + n1
314 *
315  ELSE
316 *
317 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
318 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
319 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
320 *
321  CALL zpotrf( 'L', n1, a( n2 ), n, info )
322  IF( info.GT.0 )
323  $ RETURN
324  CALL ztrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ), n,
325  $ a( 0 ), n )
326  CALL zherk( 'U', 'C', n2, n1, -one, a( 0 ), n, one,
327  $ a( n1 ), n )
328  CALL zpotrf( 'U', n2, a( n1 ), n, info )
329  IF( info.GT.0 )
330  $ info = info + n1
331 *
332  END IF
333 *
334  ELSE
335 *
336 * N is odd and TRANSR = 'C'
337 *
338  IF( lower ) THEN
339 *
340 * SRPA for LOWER, TRANSPOSE and N is odd
341 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
342 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
343 *
344  CALL zpotrf( 'U', n1, a( 0 ), n1, info )
345  IF( info.GT.0 )
346  $ RETURN
347  CALL ztrsm( 'L', 'U', 'C', 'N', n1, n2, cone, a( 0 ), n1,
348  $ a( n1*n1 ), n1 )
349  CALL zherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1, one,
350  $ a( 1 ), n1 )
351  CALL zpotrf( 'L', n2, a( 1 ), n1, info )
352  IF( info.GT.0 )
353  $ info = info + n1
354 *
355  ELSE
356 *
357 * SRPA for UPPER, TRANSPOSE and N is odd
358 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
359 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
360 *
361  CALL zpotrf( 'U', n1, a( n2*n2 ), n2, info )
362  IF( info.GT.0 )
363  $ RETURN
364  CALL ztrsm( 'R', 'U', 'N', 'N', n2, n1, cone, a( n2*n2 ),
365  $ n2, a( 0 ), n2 )
366  CALL zherk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
367  $ a( n1*n2 ), n2 )
368  CALL zpotrf( 'L', n2, a( n1*n2 ), n2, info )
369  IF( info.GT.0 )
370  $ info = info + n1
371 *
372  END IF
373 *
374  END IF
375 *
376  ELSE
377 *
378 * N is even
379 *
380  IF( normaltransr ) THEN
381 *
382 * N is even and TRANSR = 'N'
383 *
384  IF( lower ) THEN
385 *
386 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
387 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
388 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
389 *
390  CALL zpotrf( 'L', k, a( 1 ), n+1, info )
391  IF( info.GT.0 )
392  $ RETURN
393  CALL ztrsm( 'R', 'L', 'C', 'N', k, k, cone, a( 1 ), n+1,
394  $ a( k+1 ), n+1 )
395  CALL zherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
396  $ a( 0 ), n+1 )
397  CALL zpotrf( 'U', k, a( 0 ), n+1, info )
398  IF( info.GT.0 )
399  $ info = info + k
400 *
401  ELSE
402 *
403 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
404 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
405 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
406 *
407  CALL zpotrf( 'L', k, a( k+1 ), n+1, info )
408  IF( info.GT.0 )
409  $ RETURN
410  CALL ztrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
411  $ n+1, a( 0 ), n+1 )
412  CALL zherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
413  $ a( k ), n+1 )
414  CALL zpotrf( 'U', k, a( k ), n+1, info )
415  IF( info.GT.0 )
416  $ info = info + k
417 *
418  END IF
419 *
420  ELSE
421 *
422 * N is even and TRANSR = 'C'
423 *
424  IF( lower ) THEN
425 *
426 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
427 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
428 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
429 *
430  CALL zpotrf( 'U', k, a( 0+k ), k, info )
431  IF( info.GT.0 )
432  $ RETURN
433  CALL ztrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ), n1,
434  $ a( k*( k+1 ) ), k )
435  CALL zherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k, one,
436  $ a( 0 ), k )
437  CALL zpotrf( 'L', k, a( 0 ), k, info )
438  IF( info.GT.0 )
439  $ info = info + k
440 *
441  ELSE
442 *
443 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
444 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
445 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
446 *
447  CALL zpotrf( 'U', k, a( k*( k+1 ) ), k, info )
448  IF( info.GT.0 )
449  $ RETURN
450  CALL ztrsm( 'R', 'U', 'N', 'N', k, k, cone,
451  $ a( k*( k+1 ) ), k, a( 0 ), k )
452  CALL zherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
453  $ a( k*k ), k )
454  CALL zpotrf( 'L', k, a( k*k ), k, info )
455  IF( info.GT.0 )
456  $ info = info + k
457 *
458  END IF
459 *
460  END IF
461 *
462  END IF
463 *
464  RETURN
465 *
466 * End of ZPFTRF
467 *
468  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:173
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
subroutine zpftrf(TRANSR, UPLO, N, A, INFO)
ZPFTRF
Definition: zpftrf.f:211
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102