LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine chetrs_3 ( character  UPLO,
integer  N,
integer  NRHS,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( * )  E,
integer, dimension( * )  IPIV,
complex, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

CHETRS_3

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

Purpose: 
 CHETRS_3 solves a system of linear equations A * X = B with a complex
 Hermitian matrix A using the factorization computed
 by CHETRF_RK or CHETRF_BK:

    A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),

 where U (or L) is unit upper (or lower) triangular matrix,
 U**H (or L**H) is the conjugate of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is Hermitian and block
 diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This algorithm is using Level 3 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are
          stored as an upper or lower triangular matrix:
          = 'U':  Upper triangular, form is A = P*U*D*(U**H)*(P**T);
          = 'L':  Lower triangular, form is A = P*L*D*(L**H)*(P**T).
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.
[in]A
          A is COMPLEX array, dimension (LDA,N)
          Diagonal of the block diagonal matrix D and factors U or L
          as computed by CHETRF_RK and CHETRF_BK:
            a) ONLY diagonal elements of the Hermitian block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                should be provided on entry in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]E
          E is COMPLEX array, dimension (N)
          On entry, contains the superdiagonal (or subdiagonal)
          elements of the Hermitian block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not refernced;
          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is not referenced in both
          UPLO = 'U' or UPLO = 'L' cases.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by CHETRF_RK or CHETRF_BK.
[in,out]B
          B is COMPLEX array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors: 
  December 2016,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 167 of file chetrs_3.f.

167 *
168 * -- LAPACK computational routine (version 3.7.0) --
169 * -- LAPACK is a software package provided by Univ. of Tennessee, --
170 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171 * December 2016
172 *
173 * .. Scalar Arguments ..
174  CHARACTER uplo
175  INTEGER info, lda, ldb, n, nrhs
176 * ..
177 * .. Array Arguments ..
178  INTEGER ipiv( * )
179  COMPLEX a( lda, * ), b( ldb, * ), e( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  COMPLEX one
186  parameter ( one = ( 1.0e+0,0.0e+0 ) )
187 * ..
188 * .. Local Scalars ..
189  LOGICAL upper
190  INTEGER i, j, k, kp
191  REAL s
192  COMPLEX ak, akm1, akm1k, bk, bkm1, denom
193 * ..
194 * .. External Functions ..
195  LOGICAL lsame
196  EXTERNAL lsame
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL csscal, cswap, ctrsm, xerbla
200 * ..
201 * .. Intrinsic Functions ..
202  INTRINSIC abs, conjg, max, real
203 * ..
204 * .. Executable Statements ..
205 *
206  info = 0
207  upper = lsame( uplo, 'U' )
208  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
209  info = -1
210  ELSE IF( n.LT.0 ) THEN
211  info = -2
212  ELSE IF( nrhs.LT.0 ) THEN
213  info = -3
214  ELSE IF( lda.LT.max( 1, n ) ) THEN
215  info = -5
216  ELSE IF( ldb.LT.max( 1, n ) ) THEN
217  info = -9
218  END IF
219  IF( info.NE.0 ) THEN
220  CALL xerbla( 'CHETRS_3', -info )
221  RETURN
222  END IF
223 *
224 * Quick return if possible
225 *
226  IF( n.EQ.0 .OR. nrhs.EQ.0 )
227  $ RETURN
228 *
229  IF( upper ) THEN
230 *
231 * Begin Upper
232 *
233 * Solve A*X = B, where A = U*D*U**H.
234 *
235 * P**T * B
236 *
237 * Interchange rows K and IPIV(K) of matrix B in the same order
238 * that the formation order of IPIV(I) vector for Upper case.
239 *
240 * (We can do the simple loop over IPIV with decrement -1,
241 * since the ABS value of IPIV(I) represents the row index
242 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
243 *
244  DO k = n, 1, -1
245  kp = abs( ipiv( k ) )
246  IF( kp.NE.k ) THEN
247  CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
248  END IF
249  END DO
250 *
251 * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
252 *
253  CALL ctrsm( 'L', 'U', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
254 *
255 * Compute D \ B -> B [ D \ (U \P**T * B) ]
256 *
257  i = n
258  DO WHILE ( i.GE.1 )
259  IF( ipiv( i ).GT.0 ) THEN
260  s = REAL( ONE ) / REAL( A( I, I ) )
261  CALL csscal( nrhs, s, b( i, 1 ), ldb )
262  ELSE IF ( i.GT.1 ) THEN
263  akm1k = e( i )
264  akm1 = a( i-1, i-1 ) / akm1k
265  ak = a( i, i ) / conjg( akm1k )
266  denom = akm1*ak - one
267  DO j = 1, nrhs
268  bkm1 = b( i-1, j ) / akm1k
269  bk = b( i, j ) / conjg( akm1k )
270  b( i-1, j ) = ( ak*bkm1-bk ) / denom
271  b( i, j ) = ( akm1*bk-bkm1 ) / denom
272  END DO
273  i = i - 1
274  END IF
275  i = i - 1
276  END DO
277 *
278 * Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
279 *
280  CALL ctrsm( 'L', 'U', 'C', 'U', n, nrhs, one, a, lda, b, ldb )
281 *
282 * P * B [ P * (U**H \ (D \ (U \P**T * B) )) ]
283 *
284 * Interchange rows K and IPIV(K) of matrix B in reverse order
285 * from the formation order of IPIV(I) vector for Upper case.
286 *
287 * (We can do the simple loop over IPIV with increment 1,
288 * since the ABS value of IPIV(I) represents the row index
289 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
290 *
291  DO k = 1, n, 1
292  kp = abs( ipiv( k ) )
293  IF( kp.NE.k ) THEN
294  CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
295  END IF
296  END DO
297 *
298  ELSE
299 *
300 * Begin Lower
301 *
302 * Solve A*X = B, where A = L*D*L**H.
303 *
304 * P**T * B
305 * Interchange rows K and IPIV(K) of matrix B in the same order
306 * that the formation order of IPIV(I) vector for Lower case.
307 *
308 * (We can do the simple loop over IPIV with increment 1,
309 * since the ABS value of IPIV(I) represents the row index
310 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
311 *
312  DO k = 1, n, 1
313  kp = abs( ipiv( k ) )
314  IF( kp.NE.k ) THEN
315  CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
316  END IF
317  END DO
318 *
319 * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
320 *
321  CALL ctrsm( 'L', 'L', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
322 *
323 * Compute D \ B -> B [ D \ (L \P**T * B) ]
324 *
325  i = 1
326  DO WHILE ( i.LE.n )
327  IF( ipiv( i ).GT.0 ) THEN
328  s = REAL( ONE ) / REAL( A( I, I ) )
329  CALL csscal( nrhs, s, b( i, 1 ), ldb )
330  ELSE IF( i.LT.n ) THEN
331  akm1k = e( i )
332  akm1 = a( i, i ) / conjg( akm1k )
333  ak = a( i+1, i+1 ) / akm1k
334  denom = akm1*ak - one
335  DO j = 1, nrhs
336  bkm1 = b( i, j ) / conjg( akm1k )
337  bk = b( i+1, j ) / akm1k
338  b( i, j ) = ( ak*bkm1-bk ) / denom
339  b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
340  END DO
341  i = i + 1
342  END IF
343  i = i + 1
344  END DO
345 *
346 * Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
347 *
348  CALL ctrsm('L', 'L', 'C', 'U', n, nrhs, one, a, lda, b, ldb )
349 *
350 * P * B [ P * (L**H \ (D \ (L \P**T * B) )) ]
351 *
352 * Interchange rows K and IPIV(K) of matrix B in reverse order
353 * from the formation order of IPIV(I) vector for Lower case.
354 *
355 * (We can do the simple loop over IPIV with decrement -1,
356 * since the ABS value of IPIV(I) represents the row index
357 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
358 *
359  DO k = n, 1, -1
360  kp = abs( ipiv( k ) )
361  IF( kp.NE.k ) THEN
362  CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
363  END IF
364  END DO
365 *
366 * END Lower
367 *
368  END IF
369 *
370  RETURN
371 *
372 * End of CHETRS_3
373 *
ctrsm
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: