LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cpftrf.f
Go to the documentation of this file.
1 *> \brief \b CPFTRF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CPFTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpftrf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpftrf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpftrf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CPFTRF( TRANSR, UPLO, N, A, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER TRANSR, UPLO
25 * INTEGER N, INFO
26 * ..
27 * .. Array Arguments ..
28 * COMPLEX A( 0: * )
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> CPFTRF 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 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 *> \date December 2016
208 *
209 *> \ingroup complexOTHERcomputational
210 *
211 * =====================================================================
212  SUBROUTINE cpftrf( TRANSR, UPLO, N, A, INFO )
213 *
214 * -- LAPACK computational routine (version 3.7.0) --
215 * -- LAPACK is a software package provided by Univ. of Tennessee, --
216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * December 2016
218 *
219 * .. Scalar Arguments ..
220  CHARACTER TRANSR, UPLO
221  INTEGER N, INFO
222 * ..
223 * .. Array Arguments ..
224  COMPLEX A( 0: * )
225 *
226 * =====================================================================
227 *
228 * .. Parameters ..
229  REAL ONE
230  COMPLEX CONE
231  parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+0 ) )
232 * ..
233 * .. Local Scalars ..
234  LOGICAL LOWER, NISODD, NORMALTRANSR
235  INTEGER N1, N2, K
236 * ..
237 * .. External Functions ..
238  LOGICAL LSAME
239  EXTERNAL lsame
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL xerbla, cherk, cpotrf, ctrsm
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC mod
246 * ..
247 * .. Executable Statements ..
248 *
249 * Test the input parameters.
250 *
251  info = 0
252  normaltransr = lsame( transr, 'N' )
253  lower = lsame( uplo, 'L' )
254  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
255  info = -1
256  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
257  info = -2
258  ELSE IF( n.LT.0 ) THEN
259  info = -3
260  END IF
261  IF( info.NE.0 ) THEN
262  CALL xerbla( 'CPFTRF', -info )
263  RETURN
264  END IF
265 *
266 * Quick return if possible
267 *
268  IF( n.EQ.0 )
269  $ RETURN
270 *
271 * If N is odd, set NISODD = .TRUE.
272 * If N is even, set K = N/2 and NISODD = .FALSE.
273 *
274  IF( mod( n, 2 ).EQ.0 ) THEN
275  k = n / 2
276  nisodd = .false.
277  ELSE
278  nisodd = .true.
279  END IF
280 *
281 * Set N1 and N2 depending on LOWER
282 *
283  IF( lower ) THEN
284  n2 = n / 2
285  n1 = n - n2
286  ELSE
287  n1 = n / 2
288  n2 = n - n1
289  END IF
290 *
291 * start execution: there are eight cases
292 *
293  IF( nisodd ) THEN
294 *
295 * N is odd
296 *
297  IF( normaltransr ) THEN
298 *
299 * N is odd and TRANSR = 'N'
300 *
301  IF( lower ) THEN
302 *
303 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
304 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
305 * T1 -> a(0), T2 -> a(n), S -> a(n1)
306 *
307  CALL cpotrf( 'L', n1, a( 0 ), n, info )
308  IF( info.GT.0 )
309  $ RETURN
310  CALL ctrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ), n,
311  $ a( n1 ), n )
312  CALL cherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
313  $ a( n ), n )
314  CALL cpotrf( 'U', n2, a( n ), n, info )
315  IF( info.GT.0 )
316  $ info = info + n1
317 *
318  ELSE
319 *
320 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
321 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
322 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
323 *
324  CALL cpotrf( 'L', n1, a( n2 ), n, info )
325  IF( info.GT.0 )
326  $ RETURN
327  CALL ctrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ), n,
328  $ a( 0 ), n )
329  CALL cherk( 'U', 'C', n2, n1, -one, a( 0 ), n, one,
330  $ a( n1 ), n )
331  CALL cpotrf( 'U', n2, a( n1 ), n, info )
332  IF( info.GT.0 )
333  $ info = info + n1
334 *
335  END IF
336 *
337  ELSE
338 *
339 * N is odd and TRANSR = 'C'
340 *
341  IF( lower ) THEN
342 *
343 * SRPA for LOWER, TRANSPOSE and N is odd
344 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
345 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
346 *
347  CALL cpotrf( 'U', n1, a( 0 ), n1, info )
348  IF( info.GT.0 )
349  $ RETURN
350  CALL ctrsm( 'L', 'U', 'C', 'N', n1, n2, cone, a( 0 ), n1,
351  $ a( n1*n1 ), n1 )
352  CALL cherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1, one,
353  $ a( 1 ), n1 )
354  CALL cpotrf( 'L', n2, a( 1 ), n1, info )
355  IF( info.GT.0 )
356  $ info = info + n1
357 *
358  ELSE
359 *
360 * SRPA for UPPER, TRANSPOSE and N is odd
361 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
362 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
363 *
364  CALL cpotrf( 'U', n1, a( n2*n2 ), n2, info )
365  IF( info.GT.0 )
366  $ RETURN
367  CALL ctrsm( 'R', 'U', 'N', 'N', n2, n1, cone, a( n2*n2 ),
368  $ n2, a( 0 ), n2 )
369  CALL cherk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
370  $ a( n1*n2 ), n2 )
371  CALL cpotrf( 'L', n2, a( n1*n2 ), n2, info )
372  IF( info.GT.0 )
373  $ info = info + n1
374 *
375  END IF
376 *
377  END IF
378 *
379  ELSE
380 *
381 * N is even
382 *
383  IF( normaltransr ) THEN
384 *
385 * N is even and TRANSR = 'N'
386 *
387  IF( lower ) THEN
388 *
389 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
390 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
391 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
392 *
393  CALL cpotrf( 'L', k, a( 1 ), n+1, info )
394  IF( info.GT.0 )
395  $ RETURN
396  CALL ctrsm( 'R', 'L', 'C', 'N', k, k, cone, a( 1 ), n+1,
397  $ a( k+1 ), n+1 )
398  CALL cherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
399  $ a( 0 ), n+1 )
400  CALL cpotrf( 'U', k, a( 0 ), n+1, info )
401  IF( info.GT.0 )
402  $ info = info + k
403 *
404  ELSE
405 *
406 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
407 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
408 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
409 *
410  CALL cpotrf( 'L', k, a( k+1 ), n+1, info )
411  IF( info.GT.0 )
412  $ RETURN
413  CALL ctrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
414  $ n+1, a( 0 ), n+1 )
415  CALL cherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
416  $ a( k ), n+1 )
417  CALL cpotrf( 'U', k, a( k ), n+1, info )
418  IF( info.GT.0 )
419  $ info = info + k
420 *
421  END IF
422 *
423  ELSE
424 *
425 * N is even and TRANSR = 'C'
426 *
427  IF( lower ) THEN
428 *
429 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
430 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
431 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
432 *
433  CALL cpotrf( 'U', k, a( 0+k ), k, info )
434  IF( info.GT.0 )
435  $ RETURN
436  CALL ctrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ), n1,
437  $ a( k*( k+1 ) ), k )
438  CALL cherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k, one,
439  $ a( 0 ), k )
440  CALL cpotrf( 'L', k, a( 0 ), k, info )
441  IF( info.GT.0 )
442  $ info = info + k
443 *
444  ELSE
445 *
446 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
447 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
448 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
449 *
450  CALL cpotrf( 'U', k, a( k*( k+1 ) ), k, info )
451  IF( info.GT.0 )
452  $ RETURN
453  CALL ctrsm( 'R', 'U', 'N', 'N', k, k, cone,
454  $ a( k*( k+1 ) ), k, a( 0 ), k )
455  CALL cherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
456  $ a( k*k ), k )
457  CALL cpotrf( 'L', k, a( k*k ), k, info )
458  IF( info.GT.0 )
459  $ info = info + k
460 *
461  END IF
462 *
463  END IF
464 *
465  END IF
466 *
467  RETURN
468 *
469 * End of CPFTRF
470 *
471  END
cpftrf
subroutine cpftrf(TRANSR, UPLO, N, A, INFO)
CPFTRF
Definition: cpftrf.f:213
ctrsm
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
cpotrf
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:109
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62