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