LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ stftri()

subroutine stftri ( character  transr,
character  uplo,
character  diag,
integer  n,
real, dimension( 0: * )  a,
integer  info 
)

STFTRI

Download STFTRI + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STFTRI computes the inverse of a triangular matrix A stored in RFP
 format.

 This is a Level 3 BLAS version of the algorithm.
Parameters
[in]TRANSR
          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.
[in]DIAG
          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension (NT);
          NT=N*(N+1)/2. On entry, the triangular factor of a Hermitian
          Positive Definite matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the nt elements of
          upper packed A; If UPLO = 'L' the RFP A contains the nt
          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
          TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
          even and N is odd. See the Note below for more details.

          On exit, the (triangular) inverse of the original matrix, in
          the same storage format.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
               matrix is singular and its inverse can not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

Definition at line 200 of file stftri.f.

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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
strtri
subroutine strtri(uplo, diag, n, a, lda, info)
STRTRI
Definition strtri.f:109
Here is the call graph for this function:
Here is the caller graph for this function: