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