LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup complexOTHERcomputational
110 *
111 *> \par Further Details:
112 * =====================
113 *>
114 *> \verbatim
115 *>
116 *> We first consider Standard Packed Format when N is even.
117 *> We give an example where N = 6.
118 *>
119 *> AP is Upper AP is Lower
120 *>
121 *> 00 01 02 03 04 05 00
122 *> 11 12 13 14 15 10 11
123 *> 22 23 24 25 20 21 22
124 *> 33 34 35 30 31 32 33
125 *> 44 45 40 41 42 43 44
126 *> 55 50 51 52 53 54 55
127 *>
128 *>
129 *> Let TRANSR = 'N'. RFP holds AP as follows:
130 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
131 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
132 *> conjugate-transpose of the first three columns of AP upper.
133 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
134 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
135 *> conjugate-transpose of the last three columns of AP lower.
136 *> To denote conjugate we place -- above the element. This covers the
137 *> case N even and TRANSR = 'N'.
138 *>
139 *> RFP A RFP A
140 *>
141 *> -- -- --
142 *> 03 04 05 33 43 53
143 *> -- --
144 *> 13 14 15 00 44 54
145 *> --
146 *> 23 24 25 10 11 55
147 *>
148 *> 33 34 35 20 21 22
149 *> --
150 *> 00 44 45 30 31 32
151 *> -- --
152 *> 01 11 55 40 41 42
153 *> -- -- --
154 *> 02 12 22 50 51 52
155 *>
156 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
157 *> transpose of RFP A above. One therefore gets:
158 *>
159 *>
160 *> RFP A RFP A
161 *>
162 *> -- -- -- -- -- -- -- -- -- --
163 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
164 *> -- -- -- -- -- -- -- -- -- --
165 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
166 *> -- -- -- -- -- -- -- -- -- --
167 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
168 *>
169 *>
170 *> We next consider Standard Packed Format when N is odd.
171 *> We give an example where N = 5.
172 *>
173 *> AP is Upper AP is Lower
174 *>
175 *> 00 01 02 03 04 00
176 *> 11 12 13 14 10 11
177 *> 22 23 24 20 21 22
178 *> 33 34 30 31 32 33
179 *> 44 40 41 42 43 44
180 *>
181 *>
182 *> Let TRANSR = 'N'. RFP holds AP as follows:
183 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
184 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
185 *> conjugate-transpose of the first two columns of AP upper.
186 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
187 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
188 *> conjugate-transpose of the last two columns of AP lower.
189 *> To denote conjugate we place -- above the element. This covers the
190 *> case N odd and TRANSR = 'N'.
191 *>
192 *> RFP A RFP A
193 *>
194 *> -- --
195 *> 02 03 04 00 33 43
196 *> --
197 *> 12 13 14 10 11 44
198 *>
199 *> 22 23 24 20 21 22
200 *> --
201 *> 00 33 34 30 31 32
202 *> -- --
203 *> 01 11 44 40 41 42
204 *>
205 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
206 *> transpose of RFP A above. One therefore gets:
207 *>
208 *>
209 *> RFP A RFP A
210 *>
211 *> -- -- -- -- -- -- -- -- --
212 *> 02 12 22 00 01 00 10 20 30 40 50
213 *> -- -- -- -- -- -- -- -- --
214 *> 03 13 23 33 11 33 11 21 31 41 51
215 *> -- -- -- -- -- -- -- -- --
216 *> 04 14 24 34 44 43 44 22 32 42 52
217 *> \endverbatim
218 *>
219 * =====================================================================
220  SUBROUTINE ctftri( TRANSR, UPLO, DIAG, N, A, INFO )
221 *
222 * -- LAPACK computational routine --
223 * -- LAPACK is a software package provided by Univ. of Tennessee, --
224 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225 *
226 * .. Scalar Arguments ..
227  CHARACTER TRANSR, UPLO, DIAG
228  INTEGER INFO, N
229 * ..
230 * .. Array Arguments ..
231  COMPLEX A( 0: * )
232 * ..
233 *
234 * =====================================================================
235 *
236 * .. Parameters ..
237  COMPLEX CONE
238  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
239 * ..
240 * .. Local Scalars ..
241  LOGICAL LOWER, NISODD, NORMALTRANSR
242  INTEGER N1, N2, K
243 * ..
244 * .. External Functions ..
245  LOGICAL LSAME
246  EXTERNAL lsame
247 * ..
248 * .. External Subroutines ..
249  EXTERNAL xerbla, ctrmm, ctrtri
250 * ..
251 * .. Intrinsic Functions ..
252  INTRINSIC mod
253 * ..
254 * .. Executable Statements ..
255 *
256 * Test the input parameters.
257 *
258  info = 0
259  normaltransr = lsame( transr, 'N' )
260  lower = lsame( uplo, 'L' )
261  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
262  info = -1
263  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
264  info = -2
265  ELSE IF( .NOT.lsame( diag, 'N' ) .AND. .NOT.lsame( diag, 'U' ) )
266  $ THEN
267  info = -3
268  ELSE IF( n.LT.0 ) THEN
269  info = -4
270  END IF
271  IF( info.NE.0 ) THEN
272  CALL xerbla( 'CTFTRI', -info )
273  RETURN
274  END IF
275 *
276 * Quick return if possible
277 *
278  IF( n.EQ.0 )
279  $ RETURN
280 *
281 * If N is odd, set NISODD = .TRUE.
282 * If N is even, set K = N/2 and NISODD = .FALSE.
283 *
284  IF( mod( n, 2 ).EQ.0 ) THEN
285  k = n / 2
286  nisodd = .false.
287  ELSE
288  nisodd = .true.
289  END IF
290 *
291 * Set N1 and N2 depending on LOWER
292 *
293  IF( lower ) THEN
294  n2 = n / 2
295  n1 = n - n2
296  ELSE
297  n1 = n / 2
298  n2 = n - n1
299  END IF
300 *
301 *
302 * start execution: there are eight cases
303 *
304  IF( nisodd ) THEN
305 *
306 * N is odd
307 *
308  IF( normaltransr ) THEN
309 *
310 * N is odd and TRANSR = 'N'
311 *
312  IF( lower ) THEN
313 *
314 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
315 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
316 * T1 -> a(0), T2 -> a(n), S -> a(n1)
317 *
318  CALL ctrtri( 'L', diag, n1, a( 0 ), n, info )
319  IF( info.GT.0 )
320  $ RETURN
321  CALL ctrmm( 'R', 'L', 'N', diag, n2, n1, -cone, a( 0 ),
322  $ n, a( n1 ), n )
323  CALL ctrtri( 'U', diag, n2, a( n ), n, info )
324  IF( info.GT.0 )
325  $ info = info + n1
326  IF( info.GT.0 )
327  $ RETURN
328  CALL ctrmm( 'L', 'U', 'C', diag, n2, n1, cone, a( n ), n,
329  $ a( n1 ), n )
330 *
331  ELSE
332 *
333 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
334 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
335 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
336 *
337  CALL ctrtri( 'L', diag, n1, a( n2 ), n, info )
338  IF( info.GT.0 )
339  $ RETURN
340  CALL ctrmm( 'L', 'L', 'C', diag, n1, n2, -cone, a( n2 ),
341  $ n, a( 0 ), n )
342  CALL ctrtri( 'U', diag, n2, a( n1 ), n, info )
343  IF( info.GT.0 )
344  $ info = info + n1
345  IF( info.GT.0 )
346  $ RETURN
347  CALL ctrmm( 'R', 'U', 'N', diag, n1, n2, cone, a( n1 ),
348  $ n, a( 0 ), n )
349 *
350  END IF
351 *
352  ELSE
353 *
354 * N is odd and TRANSR = 'C'
355 *
356  IF( lower ) THEN
357 *
358 * SRPA for LOWER, TRANSPOSE and N is odd
359 * T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
360 *
361  CALL ctrtri( 'U', diag, n1, a( 0 ), n1, info )
362  IF( info.GT.0 )
363  $ RETURN
364  CALL ctrmm( 'L', 'U', 'N', diag, n1, n2, -cone, a( 0 ),
365  $ n1, a( n1*n1 ), n1 )
366  CALL ctrtri( 'L', diag, n2, a( 1 ), n1, info )
367  IF( info.GT.0 )
368  $ info = info + n1
369  IF( info.GT.0 )
370  $ RETURN
371  CALL ctrmm( 'R', 'L', 'C', diag, n1, n2, cone, a( 1 ),
372  $ n1, a( n1*n1 ), n1 )
373 *
374  ELSE
375 *
376 * SRPA for UPPER, TRANSPOSE and N is odd
377 * T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
378 *
379  CALL ctrtri( 'U', diag, n1, a( n2*n2 ), n2, info )
380  IF( info.GT.0 )
381  $ RETURN
382  CALL ctrmm( 'R', 'U', 'C', diag, n2, n1, -cone,
383  $ a( n2*n2 ), n2, a( 0 ), n2 )
384  CALL ctrtri( 'L', diag, n2, a( n1*n2 ), n2, info )
385  IF( info.GT.0 )
386  $ info = info + n1
387  IF( info.GT.0 )
388  $ RETURN
389  CALL ctrmm( 'L', 'L', 'N', diag, n2, n1, cone,
390  $ a( n1*n2 ), n2, a( 0 ), n2 )
391  END IF
392 *
393  END IF
394 *
395  ELSE
396 *
397 * N is even
398 *
399  IF( normaltransr ) THEN
400 *
401 * N is even and TRANSR = 'N'
402 *
403  IF( lower ) THEN
404 *
405 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
406 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
407 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
408 *
409  CALL ctrtri( 'L', diag, k, a( 1 ), n+1, info )
410  IF( info.GT.0 )
411  $ RETURN
412  CALL ctrmm( 'R', 'L', 'N', diag, k, k, -cone, a( 1 ),
413  $ n+1, a( k+1 ), n+1 )
414  CALL ctrtri( 'U', diag, k, a( 0 ), n+1, info )
415  IF( info.GT.0 )
416  $ info = info + k
417  IF( info.GT.0 )
418  $ RETURN
419  CALL ctrmm( 'L', 'U', 'C', diag, k, k, cone, a( 0 ), n+1,
420  $ a( k+1 ), n+1 )
421 *
422  ELSE
423 *
424 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
425 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
426 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
427 *
428  CALL ctrtri( 'L', diag, k, a( k+1 ), n+1, info )
429  IF( info.GT.0 )
430  $ RETURN
431  CALL ctrmm( 'L', 'L', 'C', diag, k, k, -cone, a( k+1 ),
432  $ n+1, a( 0 ), n+1 )
433  CALL ctrtri( 'U', diag, k, a( k ), n+1, info )
434  IF( info.GT.0 )
435  $ info = info + k
436  IF( info.GT.0 )
437  $ RETURN
438  CALL ctrmm( 'R', 'U', 'N', diag, k, k, cone, a( k ), n+1,
439  $ a( 0 ), n+1 )
440  END IF
441  ELSE
442 *
443 * N is even and TRANSR = 'C'
444 *
445  IF( lower ) THEN
446 *
447 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
448 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
449 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
450 *
451  CALL ctrtri( 'U', diag, k, a( k ), k, info )
452  IF( info.GT.0 )
453  $ RETURN
454  CALL ctrmm( 'L', 'U', 'N', diag, k, k, -cone, a( k ), k,
455  $ a( k*( k+1 ) ), k )
456  CALL ctrtri( 'L', diag, k, a( 0 ), k, info )
457  IF( info.GT.0 )
458  $ info = info + k
459  IF( info.GT.0 )
460  $ RETURN
461  CALL ctrmm( 'R', 'L', 'C', diag, k, k, cone, a( 0 ), k,
462  $ a( k*( k+1 ) ), k )
463  ELSE
464 *
465 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
466 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
467 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
468 *
469  CALL ctrtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
470  IF( info.GT.0 )
471  $ RETURN
472  CALL ctrmm( 'R', 'U', 'C', diag, k, k, -cone,
473  $ a( k*( k+1 ) ), k, a( 0 ), k )
474  CALL ctrtri( 'L', diag, k, a( k*k ), k, info )
475  IF( info.GT.0 )
476  $ info = info + k
477  IF( info.GT.0 )
478  $ RETURN
479  CALL ctrmm( 'L', 'L', 'N', diag, k, k, cone, a( k*k ), k,
480  $ a( 0 ), k )
481  END IF
482  END IF
483  END IF
484 *
485  RETURN
486 *
487 * End of CTFTRI
488 *
489  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:177
subroutine ctftri(TRANSR, UPLO, DIAG, N, A, INFO)
CTFTRI
Definition: ctftri.f:221
subroutine ctrtri(UPLO, DIAG, N, A, LDA, INFO)
CTRTRI
Definition: ctrtri.f:109