OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex16
LAPACK » General tridiagonal Matrix » Computational routines

Functions

subroutine zgtcon (norm, n, dl, d, du, du2, ipiv, anorm, rcond, work, info)
 ZGTCON
subroutine zgtrfs (trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 ZGTRFS
subroutine zgttrf (n, dl, d, du, du2, ipiv, info)
 ZGTTRF
subroutine zgttrs (trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
 ZGTTRS
subroutine zgtts2 (itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb)
 ZGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.

Detailed Description

This is the group of complex16 computational functions for GT matrices

Function Documentation

◆ zgtcon()

subroutine zgtcon ( character norm,
integer n,
complex*16, dimension( * ) dl,
complex*16, dimension( * ) d,
complex*16, dimension( * ) du,
complex*16, dimension( * ) du2,
integer, dimension( * ) ipiv,
double precision anorm,
double precision rcond,
complex*16, dimension( * ) work,
integer info )

ZGTCON

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

Purpose:
!>
!> ZGTCON estimates the reciprocal of the condition number of a complex
!> tridiagonal matrix A using the LU factorization as computed by
!> ZGTTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
!>
Parameters
[in]NORM
!>          NORM is CHARACTER*1
!>          Specifies whether the 1-norm condition number or the
!>          infinity-norm condition number is required:
!>          = '1' or 'O':  1-norm;
!>          = 'I':         Infinity-norm.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]DL
!>          DL is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) multipliers that define the matrix L from the
!>          LU factorization of A as computed by ZGTTRF.
!>
[in]D
!>          D is COMPLEX*16 array, dimension (N)
!>          The n diagonal elements of the upper triangular matrix U from
!>          the LU factorization of A.
!>
[in]DU
!>          DU is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) elements of the first superdiagonal of U.
!>
[in]DU2
!>          DU2 is COMPLEX*16 array, dimension (N-2)
!>          The (n-2) elements of the second superdiagonal of U.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
!>          If NORM = 'I', the infinity-norm of the original matrix A.
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*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.

Definition at line 139 of file zgtcon.f.

141*
142* -- LAPACK computational routine --
143* -- LAPACK is a software package provided by Univ. of Tennessee, --
144* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145*
146* .. Scalar Arguments ..
147 CHARACTER NORM
148 INTEGER INFO, N
149 DOUBLE PRECISION ANORM, RCOND
150* ..
151* .. Array Arguments ..
152 INTEGER IPIV( * )
153 COMPLEX*16 D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 DOUBLE PRECISION ONE, ZERO
160 parameter( one = 1.0d+0, zero = 0.0d+0 )
161* ..
162* .. Local Scalars ..
163 LOGICAL ONENRM
164 INTEGER I, KASE, KASE1
165 DOUBLE PRECISION AINVNM
166* ..
167* .. Local Arrays ..
168 INTEGER ISAVE( 3 )
169* ..
170* .. External Functions ..
171 LOGICAL LSAME
172 EXTERNAL lsame
173* ..
174* .. External Subroutines ..
175 EXTERNAL xerbla, zgttrs, zlacn2
176* ..
177* .. Intrinsic Functions ..
178 INTRINSIC dcmplx
179* ..
180* .. Executable Statements ..
181*
182* Test the input arguments.
183*
184 info = 0
185 onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
186 IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
187 info = -1
188 ELSE IF( n.LT.0 ) THEN
189 info = -2
190 ELSE IF( anorm.LT.zero ) THEN
191 info = -8
192 END IF
193 IF( info.NE.0 ) THEN
194 CALL xerbla( 'ZGTCON', -info )
195 RETURN
196 END IF
197*
198* Quick return if possible
199*
200 rcond = zero
201 IF( n.EQ.0 ) THEN
202 rcond = one
203 RETURN
204 ELSE IF( anorm.EQ.zero ) THEN
205 RETURN
206 END IF
207*
208* Check that D(1:N) is non-zero.
209*
210 DO 10 i = 1, n
211 IF( d( i ).EQ.dcmplx( zero ) )
212 $ RETURN
213 10 CONTINUE
214*
215 ainvnm = zero
216 IF( onenrm ) THEN
217 kase1 = 1
218 ELSE
219 kase1 = 2
220 END IF
221 kase = 0
222 20 CONTINUE
223 CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
224 IF( kase.NE.0 ) THEN
225 IF( kase.EQ.kase1 ) THEN
226*
227* Multiply by inv(U)*inv(L).
228*
229 CALL zgttrs( 'No transpose', n, 1, dl, d, du, du2, ipiv,
230 $ work, n, info )
231 ELSE
232*
233* Multiply by inv(L**H)*inv(U**H).
234*
235 CALL zgttrs( 'Conjugate transpose', n, 1, dl, d, du, du2,
236 $ ipiv, work, n, info )
237 END IF
238 GO TO 20
239 END IF
240*
241* Compute the estimate of the reciprocal condition number.
242*
243 IF( ainvnm.NE.zero )
244 $ rcond = ( one / ainvnm ) / anorm
245*
246 RETURN
247*
248* End of ZGTCON
249*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
zgttrs
subroutine zgttrs(trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
ZGTTRS
Definition zgttrs.f:138
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133

◆ zgtrfs()

subroutine zgtrfs ( character trans,
integer n,
integer nrhs,
complex*16, dimension( * ) dl,
complex*16, dimension( * ) d,
complex*16, dimension( * ) du,
complex*16, dimension( * ) dlf,
complex*16, dimension( * ) df,
complex*16, dimension( * ) duf,
complex*16, dimension( * ) du2,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldx, * ) x,
integer ldx,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZGTRFS

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

Purpose:
!>
!> ZGTRFS improves the computed solution to a system of linear
!> equations when the coefficient matrix is tridiagonal, and provides
!> error bounds and backward error estimates for the solution.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations:
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Conjugate transpose)
!>
[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]DL
!>          DL is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) subdiagonal elements of A.
!>
[in]D
!>          D is COMPLEX*16 array, dimension (N)
!>          The diagonal elements of A.
!>
[in]DU
!>          DU is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) superdiagonal elements of A.
!>
[in]DLF
!>          DLF is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) multipliers that define the matrix L from the
!>          LU factorization of A as computed by ZGTTRF.
!>
[in]DF
!>          DF is COMPLEX*16 array, dimension (N)
!>          The n diagonal elements of the upper triangular matrix U from
!>          the LU factorization of A.
!>
[in]DUF
!>          DUF is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) elements of the first superdiagonal of U.
!>
[in]DU2
!>          DU2 is COMPLEX*16 array, dimension (N-2)
!>          The (n-2) elements of the second superdiagonal of U.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!>
[in]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          The right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]X
!>          X is COMPLEX*16 array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by ZGTTRS.
!>          On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Internal Parameters:
!>  ITMAX is the maximum number of steps of iterative refinement.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 207 of file zgtrfs.f.

210*
211* -- LAPACK computational routine --
212* -- LAPACK is a software package provided by Univ. of Tennessee, --
213* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214*
215* .. Scalar Arguments ..
216 CHARACTER TRANS
217 INTEGER INFO, LDB, LDX, N, NRHS
218* ..
219* .. Array Arguments ..
220 INTEGER IPIV( * )
221 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
222 COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ),
223 $ DLF( * ), DU( * ), DU2( * ), DUF( * ),
224 $ WORK( * ), X( LDX, * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 INTEGER ITMAX
231 parameter( itmax = 5 )
232 DOUBLE PRECISION ZERO, ONE
233 parameter( zero = 0.0d+0, one = 1.0d+0 )
234 DOUBLE PRECISION TWO
235 parameter( two = 2.0d+0 )
236 DOUBLE PRECISION THREE
237 parameter( three = 3.0d+0 )
238* ..
239* .. Local Scalars ..
240 LOGICAL NOTRAN
241 CHARACTER TRANSN, TRANST
242 INTEGER COUNT, I, J, KASE, NZ
243 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
244 COMPLEX*16 ZDUM
245* ..
246* .. Local Arrays ..
247 INTEGER ISAVE( 3 )
248* ..
249* .. External Subroutines ..
250 EXTERNAL xerbla, zaxpy, zcopy, zgttrs, zlacn2, zlagtm
251* ..
252* .. Intrinsic Functions ..
253 INTRINSIC abs, dble, dcmplx, dimag, max
254* ..
255* .. External Functions ..
256 LOGICAL LSAME
257 DOUBLE PRECISION DLAMCH
258 EXTERNAL lsame, dlamch
259* ..
260* .. Statement Functions ..
261 DOUBLE PRECISION CABS1
262* ..
263* .. Statement Function definitions ..
264 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
265* ..
266* .. Executable Statements ..
267*
268* Test the input parameters.
269*
270 info = 0
271 notran = lsame( trans, 'N' )
272 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
273 $ lsame( trans, 'C' ) ) THEN
274 info = -1
275 ELSE IF( n.LT.0 ) THEN
276 info = -2
277 ELSE IF( nrhs.LT.0 ) THEN
278 info = -3
279 ELSE IF( ldb.LT.max( 1, n ) ) THEN
280 info = -13
281 ELSE IF( ldx.LT.max( 1, n ) ) THEN
282 info = -15
283 END IF
284 IF( info.NE.0 ) THEN
285 CALL xerbla( 'ZGTRFS', -info )
286 RETURN
287 END IF
288*
289* Quick return if possible
290*
291 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
292 DO 10 j = 1, nrhs
293 ferr( j ) = zero
294 berr( j ) = zero
295 10 CONTINUE
296 RETURN
297 END IF
298*
299 IF( notran ) THEN
300 transn = 'N'
301 transt = 'C'
302 ELSE
303 transn = 'C'
304 transt = 'N'
305 END IF
306*
307* NZ = maximum number of nonzero elements in each row of A, plus 1
308*
309 nz = 4
310 eps = dlamch( 'Epsilon' )
311 safmin = dlamch( 'Safe minimum' )
312 safe1 = nz*safmin
313 safe2 = safe1 / eps
314*
315* Do for each right hand side
316*
317 DO 110 j = 1, nrhs
318*
319 count = 1
320 lstres = three
321 20 CONTINUE
322*
323* Loop until stopping criterion is satisfied.
324*
325* Compute residual R = B - op(A) * X,
326* where op(A) = A, A**T, or A**H, depending on TRANS.
327*
328 CALL zcopy( n, b( 1, j ), 1, work, 1 )
329 CALL zlagtm( trans, n, 1, -one, dl, d, du, x( 1, j ), ldx, one,
330 $ work, n )
331*
332* Compute abs(op(A))*abs(x) + abs(b) for use in the backward
333* error bound.
334*
335 IF( notran ) THEN
336 IF( n.EQ.1 ) THEN
337 rwork( 1 ) = cabs1( b( 1, j ) ) +
338 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
339 ELSE
340 rwork( 1 ) = cabs1( b( 1, j ) ) +
341 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
342 $ cabs1( du( 1 ) )*cabs1( x( 2, j ) )
343 DO 30 i = 2, n - 1
344 rwork( i ) = cabs1( b( i, j ) ) +
345 $ cabs1( dl( i-1 ) )*cabs1( x( i-1, j ) ) +
346 $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
347 $ cabs1( du( i ) )*cabs1( x( i+1, j ) )
348 30 CONTINUE
349 rwork( n ) = cabs1( b( n, j ) ) +
350 $ cabs1( dl( n-1 ) )*cabs1( x( n-1, j ) ) +
351 $ cabs1( d( n ) )*cabs1( x( n, j ) )
352 END IF
353 ELSE
354 IF( n.EQ.1 ) THEN
355 rwork( 1 ) = cabs1( b( 1, j ) ) +
356 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
357 ELSE
358 rwork( 1 ) = cabs1( b( 1, j ) ) +
359 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
360 $ cabs1( dl( 1 ) )*cabs1( x( 2, j ) )
361 DO 40 i = 2, n - 1
362 rwork( i ) = cabs1( b( i, j ) ) +
363 $ cabs1( du( i-1 ) )*cabs1( x( i-1, j ) ) +
364 $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
365 $ cabs1( dl( i ) )*cabs1( x( i+1, j ) )
366 40 CONTINUE
367 rwork( n ) = cabs1( b( n, j ) ) +
368 $ cabs1( du( n-1 ) )*cabs1( x( n-1, j ) ) +
369 $ cabs1( d( n ) )*cabs1( x( n, j ) )
370 END IF
371 END IF
372*
373* Compute componentwise relative backward error from formula
374*
375* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
376*
377* where abs(Z) is the componentwise absolute value of the matrix
378* or vector Z. If the i-th component of the denominator is less
379* than SAFE2, then SAFE1 is added to the i-th components of the
380* numerator and denominator before dividing.
381*
382 s = zero
383 DO 50 i = 1, n
384 IF( rwork( i ).GT.safe2 ) THEN
385 s = max( s, cabs1( work( i ) ) / rwork( i ) )
386 ELSE
387 s = max( s, ( cabs1( work( i ) )+safe1 ) /
388 $ ( rwork( i )+safe1 ) )
389 END IF
390 50 CONTINUE
391 berr( j ) = s
392*
393* Test stopping criterion. Continue iterating if
394* 1) The residual BERR(J) is larger than machine epsilon, and
395* 2) BERR(J) decreased by at least a factor of 2 during the
396* last iteration, and
397* 3) At most ITMAX iterations tried.
398*
399 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
400 $ count.LE.itmax ) THEN
401*
402* Update solution and try again.
403*
404 CALL zgttrs( trans, n, 1, dlf, df, duf, du2, ipiv, work, n,
405 $ info )
406 CALL zaxpy( n, dcmplx( one ), work, 1, x( 1, j ), 1 )
407 lstres = berr( j )
408 count = count + 1
409 GO TO 20
410 END IF
411*
412* Bound error from formula
413*
414* norm(X - XTRUE) / norm(X) .le. FERR =
415* norm( abs(inv(op(A)))*
416* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
417*
418* where
419* norm(Z) is the magnitude of the largest component of Z
420* inv(op(A)) is the inverse of op(A)
421* abs(Z) is the componentwise absolute value of the matrix or
422* vector Z
423* NZ is the maximum number of nonzeros in any row of A, plus 1
424* EPS is machine epsilon
425*
426* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
427* is incremented by SAFE1 if the i-th component of
428* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
429*
430* Use ZLACN2 to estimate the infinity-norm of the matrix
431* inv(op(A)) * diag(W),
432* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
433*
434 DO 60 i = 1, n
435 IF( rwork( i ).GT.safe2 ) THEN
436 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
437 ELSE
438 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
439 $ safe1
440 END IF
441 60 CONTINUE
442*
443 kase = 0
444 70 CONTINUE
445 CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
446 IF( kase.NE.0 ) THEN
447 IF( kase.EQ.1 ) THEN
448*
449* Multiply by diag(W)*inv(op(A)**H).
450*
451 CALL zgttrs( transt, n, 1, dlf, df, duf, du2, ipiv, work,
452 $ n, info )
453 DO 80 i = 1, n
454 work( i ) = rwork( i )*work( i )
455 80 CONTINUE
456 ELSE
457*
458* Multiply by inv(op(A))*diag(W).
459*
460 DO 90 i = 1, n
461 work( i ) = rwork( i )*work( i )
462 90 CONTINUE
463 CALL zgttrs( transn, n, 1, dlf, df, duf, du2, ipiv, work,
464 $ n, info )
465 END IF
466 GO TO 70
467 END IF
468*
469* Normalize error.
470*
471 lstres = zero
472 DO 100 i = 1, n
473 lstres = max( lstres, cabs1( x( i, j ) ) )
474 100 CONTINUE
475 IF( lstres.NE.zero )
476 $ ferr( j ) = ferr( j ) / lstres
477*
478 110 CONTINUE
479*
480 RETURN
481*
482* End of ZGTRFS
483*
zlagtm
subroutine zlagtm(trans, n, nrhs, alpha, dl, d, du, x, ldx, beta, b, ldb)
ZLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix,...
Definition zlagtm.f:145
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21

◆ zgttrf()

subroutine zgttrf ( integer n,
complex*16, dimension( * ) dl,
complex*16, dimension( * ) d,
complex*16, dimension( * ) du,
complex*16, dimension( * ) du2,
integer, dimension( * ) ipiv,
integer info )

ZGTTRF

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

Purpose:
!>
!> ZGTTRF computes an LU factorization of a complex tridiagonal matrix A
!> using elimination with partial pivoting and row interchanges.
!>
!> The factorization has the form
!>    A = L * U
!> where L is a product of permutation and unit lower bidiagonal
!> matrices and U is upper triangular with nonzeros in only the main
!> diagonal and first two superdiagonals.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.
!>
[in,out]DL
!>          DL is COMPLEX*16 array, dimension (N-1)
!>          On entry, DL must contain the (n-1) sub-diagonal elements of
!>          A.
!>
!>          On exit, DL is overwritten by the (n-1) multipliers that
!>          define the matrix L from the LU factorization of A.
!>
[in,out]D
!>          D is COMPLEX*16 array, dimension (N)
!>          On entry, D must contain the diagonal elements of A.
!>
!>          On exit, D is overwritten by the n diagonal elements of the
!>          upper triangular matrix U from the LU factorization of A.
!>
[in,out]DU
!>          DU is COMPLEX*16 array, dimension (N-1)
!>          On entry, DU must contain the (n-1) super-diagonal elements
!>          of A.
!>
!>          On exit, DU is overwritten by the (n-1) elements of the first
!>          super-diagonal of U.
!>
[out]DU2
!>          DU2 is COMPLEX*16 array, dimension (N-2)
!>          On exit, DU2 is overwritten by the (n-2) elements of the
!>          second super-diagonal of U.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -k, the k-th argument had an illegal value
!>          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
!>                has been completed, but the factor U is exactly
!>                singular, and division by zero will occur if it is used
!>                to solve a system of equations.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 123 of file zgttrf.f.

124*
125* -- LAPACK computational routine --
126* -- LAPACK is a software package provided by Univ. of Tennessee, --
127* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
128*
129* .. Scalar Arguments ..
130 INTEGER INFO, N
131* ..
132* .. Array Arguments ..
133 INTEGER IPIV( * )
134 COMPLEX*16 D( * ), DL( * ), DU( * ), DU2( * )
135* ..
136*
137* =====================================================================
138*
139* .. Parameters ..
140 DOUBLE PRECISION ZERO
141 parameter( zero = 0.0d+0 )
142* ..
143* .. Local Scalars ..
144 INTEGER I
145 COMPLEX*16 FACT, TEMP, ZDUM
146* ..
147* .. External Subroutines ..
148 EXTERNAL xerbla
149* ..
150* .. Intrinsic Functions ..
151 INTRINSIC abs, dble, dimag
152* ..
153* .. Statement Functions ..
154 DOUBLE PRECISION CABS1
155* ..
156* .. Statement Function definitions ..
157 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
158* ..
159* .. Executable Statements ..
160*
161 info = 0
162 IF( n.LT.0 ) THEN
163 info = -1
164 CALL xerbla( 'ZGTTRF', -info )
165 RETURN
166 END IF
167*
168* Quick return if possible
169*
170 IF( n.EQ.0 )
171 $ RETURN
172*
173* Initialize IPIV(i) = i and DU2(i) = 0
174*
175 DO 10 i = 1, n
176 ipiv( i ) = i
177 10 CONTINUE
178 DO 20 i = 1, n - 2
179 du2( i ) = zero
180 20 CONTINUE
181*
182 DO 30 i = 1, n - 2
183 IF( cabs1( d( i ) ).GE.cabs1( dl( i ) ) ) THEN
184*
185* No row interchange required, eliminate DL(I)
186*
187 IF( cabs1( d( i ) ).NE.zero ) THEN
188 fact = dl( i ) / d( i )
189 dl( i ) = fact
190 d( i+1 ) = d( i+1 ) - fact*du( i )
191 END IF
192 ELSE
193*
194* Interchange rows I and I+1, eliminate DL(I)
195*
196 fact = d( i ) / dl( i )
197 d( i ) = dl( i )
198 dl( i ) = fact
199 temp = du( i )
200 du( i ) = d( i+1 )
201 d( i+1 ) = temp - fact*d( i+1 )
202 du2( i ) = du( i+1 )
203 du( i+1 ) = -fact*du( i+1 )
204 ipiv( i ) = i + 1
205 END IF
206 30 CONTINUE
207 IF( n.GT.1 ) THEN
208 i = n - 1
209 IF( cabs1( d( i ) ).GE.cabs1( dl( i ) ) ) THEN
210 IF( cabs1( d( i ) ).NE.zero ) THEN
211 fact = dl( i ) / d( i )
212 dl( i ) = fact
213 d( i+1 ) = d( i+1 ) - fact*du( i )
214 END IF
215 ELSE
216 fact = d( i ) / dl( i )
217 d( i ) = dl( i )
218 dl( i ) = fact
219 temp = du( i )
220 du( i ) = d( i+1 )
221 d( i+1 ) = temp - fact*d( i+1 )
222 ipiv( i ) = i + 1
223 END IF
224 END IF
225*
226* Check for a zero on the diagonal of U.
227*
228 DO 40 i = 1, n
229 IF( cabs1( d( i ) ).EQ.zero ) THEN
230 info = i
231 GO TO 50
232 END IF
233 40 CONTINUE
234 50 CONTINUE
235*
236 RETURN
237*
238* End of ZGTTRF
239*

◆ zgttrs()

subroutine zgttrs ( character trans,
integer n,
integer nrhs,
complex*16, dimension( * ) dl,
complex*16, dimension( * ) d,
complex*16, dimension( * ) du,
complex*16, dimension( * ) du2,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZGTTRS

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

Purpose:
!>
!> ZGTTRS solves one of the systems of equations
!>    A * X = B,  A**T * X = B,  or  A**H * X = B,
!> with a tridiagonal matrix A using the LU factorization computed
!> by ZGTTRF.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations.
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Conjugate transpose)
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in]DL
!>          DL is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) multipliers that define the matrix L from the
!>          LU factorization of A.
!>
[in]D
!>          D is COMPLEX*16 array, dimension (N)
!>          The n diagonal elements of the upper triangular matrix U from
!>          the LU factorization of A.
!>
[in]DU
!>          DU is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) elements of the first super-diagonal of U.
!>
[in]DU2
!>          DU2 is COMPLEX*16 array, dimension (N-2)
!>          The (n-2) elements of the second super-diagonal of U.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the matrix of right hand side vectors B.
!>          On exit, B is overwritten by the solution vectors 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 = -k, the k-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 136 of file zgttrs.f.

138*
139* -- LAPACK computational routine --
140* -- LAPACK is a software package provided by Univ. of Tennessee, --
141* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142*
143* .. Scalar Arguments ..
144 CHARACTER TRANS
145 INTEGER INFO, LDB, N, NRHS
146* ..
147* .. Array Arguments ..
148 INTEGER IPIV( * )
149 COMPLEX*16 B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
150* ..
151*
152* =====================================================================
153*
154* .. Local Scalars ..
155 LOGICAL NOTRAN
156 INTEGER ITRANS, J, JB, NB
157* ..
158* .. External Functions ..
159 INTEGER ILAENV
160 EXTERNAL ilaenv
161* ..
162* .. External Subroutines ..
163 EXTERNAL xerbla, zgtts2
164* ..
165* .. Intrinsic Functions ..
166 INTRINSIC max, min
167* ..
168* .. Executable Statements ..
169*
170 info = 0
171 notran = ( trans.EQ.'N' .OR. trans.EQ.'n' )
172 IF( .NOT.notran .AND. .NOT.( trans.EQ.'T' .OR. trans.EQ.
173 $ 't' ) .AND. .NOT.( trans.EQ.'C' .OR. trans.EQ.'c' ) ) THEN
174 info = -1
175 ELSE IF( n.LT.0 ) THEN
176 info = -2
177 ELSE IF( nrhs.LT.0 ) THEN
178 info = -3
179 ELSE IF( ldb.LT.max( n, 1 ) ) THEN
180 info = -10
181 END IF
182 IF( info.NE.0 ) THEN
183 CALL xerbla( 'ZGTTRS', -info )
184 RETURN
185 END IF
186*
187* Quick return if possible
188*
189 IF( n.EQ.0 .OR. nrhs.EQ.0 )
190 $ RETURN
191*
192* Decode TRANS
193*
194 IF( notran ) THEN
195 itrans = 0
196 ELSE IF( trans.EQ.'T' .OR. trans.EQ.'t' ) THEN
197 itrans = 1
198 ELSE
199 itrans = 2
200 END IF
201*
202* Determine the number of right-hand sides to solve at a time.
203*
204 IF( nrhs.EQ.1 ) THEN
205 nb = 1
206 ELSE
207 nb = max( 1, ilaenv( 1, 'ZGTTRS', trans, n, nrhs, -1, -1 ) )
208 END IF
209*
210 IF( nb.GE.nrhs ) THEN
211 CALL zgtts2( itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb )
212 ELSE
213 DO 10 j = 1, nrhs, nb
214 jb = min( nrhs-j+1, nb )
215 CALL zgtts2( itrans, n, jb, dl, d, du, du2, ipiv, b( 1, j ),
216 $ ldb )
217 10 CONTINUE
218 END IF
219*
220* End of ZGTTRS
221*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
zgtts2
subroutine zgtts2(itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb)
ZGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization compu...
Definition zgtts2.f:128
min
#define min(a, b)
Definition macros.h:20

◆ zgtts2()

subroutine zgtts2 ( integer itrans,
integer n,
integer nrhs,
complex*16, dimension( * ) dl,
complex*16, dimension( * ) d,
complex*16, dimension( * ) du,
complex*16, dimension( * ) du2,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb )

ZGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.

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

Purpose:
!>
!> ZGTTS2 solves one of the systems of equations
!>    A * X = B,  A**T * X = B,  or  A**H * X = B,
!> with a tridiagonal matrix A using the LU factorization computed
!> by ZGTTRF.
!>
Parameters
[in]ITRANS
!>          ITRANS is INTEGER
!>          Specifies the form of the system of equations.
!>          = 0:  A * X = B     (No transpose)
!>          = 1:  A**T * X = B  (Transpose)
!>          = 2:  A**H * X = B  (Conjugate transpose)
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in]DL
!>          DL is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) multipliers that define the matrix L from the
!>          LU factorization of A.
!>
[in]D
!>          D is COMPLEX*16 array, dimension (N)
!>          The n diagonal elements of the upper triangular matrix U from
!>          the LU factorization of A.
!>
[in]DU
!>          DU is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) elements of the first super-diagonal of U.
!>
[in]DU2
!>          DU2 is COMPLEX*16 array, dimension (N-2)
!>          The (n-2) elements of the second super-diagonal of U.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the matrix of right hand side vectors B.
!>          On exit, B is overwritten by the solution vectors X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 127 of file zgtts2.f.

128*
129* -- LAPACK computational routine --
130* -- LAPACK is a software package provided by Univ. of Tennessee, --
131* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132*
133* .. Scalar Arguments ..
134 INTEGER ITRANS, LDB, N, NRHS
135* ..
136* .. Array Arguments ..
137 INTEGER IPIV( * )
138 COMPLEX*16 B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
139* ..
140*
141* =====================================================================
142*
143* .. Local Scalars ..
144 INTEGER I, J
145 COMPLEX*16 TEMP
146* ..
147* .. Intrinsic Functions ..
148 INTRINSIC dconjg
149* ..
150* .. Executable Statements ..
151*
152* Quick return if possible
153*
154 IF( n.EQ.0 .OR. nrhs.EQ.0 )
155 $ RETURN
156*
157 IF( itrans.EQ.0 ) THEN
158*
159* Solve A*X = B using the LU factorization of A,
160* overwriting each right hand side vector with its solution.
161*
162 IF( nrhs.LE.1 ) THEN
163 j = 1
164 10 CONTINUE
165*
166* Solve L*x = b.
167*
168 DO 20 i = 1, n - 1
169 IF( ipiv( i ).EQ.i ) THEN
170 b( i+1, j ) = b( i+1, j ) - dl( i )*b( i, j )
171 ELSE
172 temp = b( i, j )
173 b( i, j ) = b( i+1, j )
174 b( i+1, j ) = temp - dl( i )*b( i, j )
175 END IF
176 20 CONTINUE
177*
178* Solve U*x = b.
179*
180 b( n, j ) = b( n, j ) / d( n )
181 IF( n.GT.1 )
182 $ b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /
183 $ d( n-1 )
184 DO 30 i = n - 2, 1, -1
185 b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-du2( i )*
186 $ b( i+2, j ) ) / d( i )
187 30 CONTINUE
188 IF( j.LT.nrhs ) THEN
189 j = j + 1
190 GO TO 10
191 END IF
192 ELSE
193 DO 60 j = 1, nrhs
194*
195* Solve L*x = b.
196*
197 DO 40 i = 1, n - 1
198 IF( ipiv( i ).EQ.i ) THEN
199 b( i+1, j ) = b( i+1, j ) - dl( i )*b( i, j )
200 ELSE
201 temp = b( i, j )
202 b( i, j ) = b( i+1, j )
203 b( i+1, j ) = temp - dl( i )*b( i, j )
204 END IF
205 40 CONTINUE
206*
207* Solve U*x = b.
208*
209 b( n, j ) = b( n, j ) / d( n )
210 IF( n.GT.1 )
211 $ b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /
212 $ d( n-1 )
213 DO 50 i = n - 2, 1, -1
214 b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-du2( i )*
215 $ b( i+2, j ) ) / d( i )
216 50 CONTINUE
217 60 CONTINUE
218 END IF
219 ELSE IF( itrans.EQ.1 ) THEN
220*
221* Solve A**T * X = B.
222*
223 IF( nrhs.LE.1 ) THEN
224 j = 1
225 70 CONTINUE
226*
227* Solve U**T * x = b.
228*
229 b( 1, j ) = b( 1, j ) / d( 1 )
230 IF( n.GT.1 )
231 $ b( 2, j ) = ( b( 2, j )-du( 1 )*b( 1, j ) ) / d( 2 )
232 DO 80 i = 3, n
233 b( i, j ) = ( b( i, j )-du( i-1 )*b( i-1, j )-du2( i-2 )*
234 $ b( i-2, j ) ) / d( i )
235 80 CONTINUE
236*
237* Solve L**T * x = b.
238*
239 DO 90 i = n - 1, 1, -1
240 IF( ipiv( i ).EQ.i ) THEN
241 b( i, j ) = b( i, j ) - dl( i )*b( i+1, j )
242 ELSE
243 temp = b( i+1, j )
244 b( i+1, j ) = b( i, j ) - dl( i )*temp
245 b( i, j ) = temp
246 END IF
247 90 CONTINUE
248 IF( j.LT.nrhs ) THEN
249 j = j + 1
250 GO TO 70
251 END IF
252 ELSE
253 DO 120 j = 1, nrhs
254*
255* Solve U**T * x = b.
256*
257 b( 1, j ) = b( 1, j ) / d( 1 )
258 IF( n.GT.1 )
259 $ b( 2, j ) = ( b( 2, j )-du( 1 )*b( 1, j ) ) / d( 2 )
260 DO 100 i = 3, n
261 b( i, j ) = ( b( i, j )-du( i-1 )*b( i-1, j )-
262 $ du2( i-2 )*b( i-2, j ) ) / d( i )
263 100 CONTINUE
264*
265* Solve L**T * x = b.
266*
267 DO 110 i = n - 1, 1, -1
268 IF( ipiv( i ).EQ.i ) THEN
269 b( i, j ) = b( i, j ) - dl( i )*b( i+1, j )
270 ELSE
271 temp = b( i+1, j )
272 b( i+1, j ) = b( i, j ) - dl( i )*temp
273 b( i, j ) = temp
274 END IF
275 110 CONTINUE
276 120 CONTINUE
277 END IF
278 ELSE
279*
280* Solve A**H * X = B.
281*
282 IF( nrhs.LE.1 ) THEN
283 j = 1
284 130 CONTINUE
285*
286* Solve U**H * x = b.
287*
288 b( 1, j ) = b( 1, j ) / dconjg( d( 1 ) )
289 IF( n.GT.1 )
290 $ b( 2, j ) = ( b( 2, j )-dconjg( du( 1 ) )*b( 1, j ) ) /
291 $ dconjg( d( 2 ) )
292 DO 140 i = 3, n
293 b( i, j ) = ( b( i, j )-dconjg( du( i-1 ) )*b( i-1, j )-
294 $ dconjg( du2( i-2 ) )*b( i-2, j ) ) /
295 $ dconjg( d( i ) )
296 140 CONTINUE
297*
298* Solve L**H * x = b.
299*
300 DO 150 i = n - 1, 1, -1
301 IF( ipiv( i ).EQ.i ) THEN
302 b( i, j ) = b( i, j ) - dconjg( dl( i ) )*b( i+1, j )
303 ELSE
304 temp = b( i+1, j )
305 b( i+1, j ) = b( i, j ) - dconjg( dl( i ) )*temp
306 b( i, j ) = temp
307 END IF
308 150 CONTINUE
309 IF( j.LT.nrhs ) THEN
310 j = j + 1
311 GO TO 130
312 END IF
313 ELSE
314 DO 180 j = 1, nrhs
315*
316* Solve U**H * x = b.
317*
318 b( 1, j ) = b( 1, j ) / dconjg( d( 1 ) )
319 IF( n.GT.1 )
320 $ b( 2, j ) = ( b( 2, j )-dconjg( du( 1 ) )*b( 1, j ) )
321 $ / dconjg( d( 2 ) )
322 DO 160 i = 3, n
323 b( i, j ) = ( b( i, j )-dconjg( du( i-1 ) )*
324 $ b( i-1, j )-dconjg( du2( i-2 ) )*
325 $ b( i-2, j ) ) / dconjg( d( i ) )
326 160 CONTINUE
327*
328* Solve L**H * x = b.
329*
330 DO 170 i = n - 1, 1, -1
331 IF( ipiv( i ).EQ.i ) THEN
332 b( i, j ) = b( i, j ) - dconjg( dl( i ) )*
333 $ b( i+1, j )
334 ELSE
335 temp = b( i+1, j )
336 b( i+1, j ) = b( i, j ) - dconjg( dl( i ) )*temp
337 b( i, j ) = temp
338 END IF
339 170 CONTINUE
340 180 CONTINUE
341 END IF
342 END IF
343*
344* End of ZGTTS2
345*