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

Functions

subroutine zptcon (n, d, e, anorm, rcond, rwork, info)
 ZPTCON
subroutine zpteqr (compz, n, d, e, z, ldz, work, info)
 ZPTEQR
subroutine zptrfs (uplo, n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 ZPTRFS
subroutine zpttrf (n, d, e, info)
 ZPTTRF
subroutine zpttrs (uplo, n, nrhs, d, e, b, ldb, info)
 ZPTTRS
subroutine zptts2 (iuplo, n, nrhs, d, e, b, ldb)
 ZPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf.

Detailed Description

This is the group of complex16 computational functions for PT matrices

Function Documentation

◆ zptcon()

subroutine zptcon ( integer n,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
double precision anorm,
double precision rcond,
double precision, dimension( * ) rwork,
integer info )

ZPTCON

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

Purpose:
!>
!> ZPTCON computes the reciprocal of the condition number (in the
!> 1-norm) of a complex Hermitian positive definite tridiagonal matrix
!> using the factorization A = L*D*L**H or A = U**H*D*U computed by
!> ZPTTRF.
!>
!> Norm(inv(A)) is computed by a direct method, and the reciprocal of
!> the condition number is computed as
!>                  RCOND = 1 / (ANORM * norm(inv(A))).
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from the
!>          factorization of A, as computed by ZPTTRF.
!>
[in]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the unit bidiagonal factor
!>          U or L from the factorization of A, as computed by ZPTTRF.
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>          The 1-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 the
!>          1-norm of inv(A) computed in this routine.
!>
[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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The method used is described in Nicholas J. Higham, , SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
!>

Definition at line 118 of file zptcon.f.

119*
120* -- LAPACK computational routine --
121* -- LAPACK is a software package provided by Univ. of Tennessee, --
122* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123*
124* .. Scalar Arguments ..
125 INTEGER INFO, N
126 DOUBLE PRECISION ANORM, RCOND
127* ..
128* .. Array Arguments ..
129 DOUBLE PRECISION D( * ), RWORK( * )
130 COMPLEX*16 E( * )
131* ..
132*
133* =====================================================================
134*
135* .. Parameters ..
136 DOUBLE PRECISION ONE, ZERO
137 parameter( one = 1.0d+0, zero = 0.0d+0 )
138* ..
139* .. Local Scalars ..
140 INTEGER I, IX
141 DOUBLE PRECISION AINVNM
142* ..
143* .. External Functions ..
144 INTEGER IDAMAX
145 EXTERNAL idamax
146* ..
147* .. External Subroutines ..
148 EXTERNAL xerbla
149* ..
150* .. Intrinsic Functions ..
151 INTRINSIC abs
152* ..
153* .. Executable Statements ..
154*
155* Test the input arguments.
156*
157 info = 0
158 IF( n.LT.0 ) THEN
159 info = -1
160 ELSE IF( anorm.LT.zero ) THEN
161 info = -4
162 END IF
163 IF( info.NE.0 ) THEN
164 CALL xerbla( 'ZPTCON', -info )
165 RETURN
166 END IF
167*
168* Quick return if possible
169*
170 rcond = zero
171 IF( n.EQ.0 ) THEN
172 rcond = one
173 RETURN
174 ELSE IF( anorm.EQ.zero ) THEN
175 RETURN
176 END IF
177*
178* Check that D(1:N) is positive.
179*
180 DO 10 i = 1, n
181 IF( d( i ).LE.zero )
182 $ RETURN
183 10 CONTINUE
184*
185* Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
186*
187* m(i,j) = abs(A(i,j)), i = j,
188* m(i,j) = -abs(A(i,j)), i .ne. j,
189*
190* and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
191*
192* Solve M(L) * x = e.
193*
194 rwork( 1 ) = one
195 DO 20 i = 2, n
196 rwork( i ) = one + rwork( i-1 )*abs( e( i-1 ) )
197 20 CONTINUE
198*
199* Solve D * M(L)**H * x = b.
200*
201 rwork( n ) = rwork( n ) / d( n )
202 DO 30 i = n - 1, 1, -1
203 rwork( i ) = rwork( i ) / d( i ) + rwork( i+1 )*abs( e( i ) )
204 30 CONTINUE
205*
206* Compute AINVNM = max(x(i)), 1<=i<=n.
207*
208 ix = idamax( n, rwork, 1 )
209 ainvnm = abs( rwork( ix ) )
210*
211* Compute the reciprocal condition number.
212*
213 IF( ainvnm.NE.zero )
214 $ rcond = ( one / ainvnm ) / anorm
215*
216 RETURN
217*
218* End of ZPTCON
219*
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71

◆ zpteqr()

subroutine zpteqr ( character compz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
complex*16, dimension( ldz, * ) z,
integer ldz,
double precision, dimension( * ) work,
integer info )

ZPTEQR

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

Purpose:
!>
!> ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a
!> symmetric positive definite tridiagonal matrix by first factoring the
!> matrix using DPTTRF and then calling ZBDSQR to compute the singular
!> values of the bidiagonal factor.
!>
!> This routine computes the eigenvalues of the positive definite
!> tridiagonal matrix to high relative accuracy.  This means that if the
!> eigenvalues range over many orders of magnitude in size, then the
!> small eigenvalues and corresponding eigenvectors will be computed
!> more accurately than, for example, with the standard QR method.
!>
!> The eigenvectors of a full or band positive definite Hermitian matrix
!> can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
!> reduce this matrix to tridiagonal form.  (The reduction to
!> tridiagonal form, however, may preclude the possibility of obtaining
!> high relative accuracy in the small eigenvalues of the original
!> matrix, if these eigenvalues range over many orders of magnitude.)
!>
Parameters
[in]COMPZ
!>          COMPZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only.
!>          = 'V':  Compute eigenvectors of original Hermitian
!>                  matrix also.  Array Z contains the unitary matrix
!>                  used to reduce the original matrix to tridiagonal
!>                  form.
!>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix.
!>          On normal exit, D contains the eigenvalues, in descending
!>          order.
!>
[in,out]E
!>          E is DOUBLE PRECISION array, dimension (N-1)
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix.
!>          On exit, E has been destroyed.
!>
[in,out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>          On entry, if COMPZ = 'V', the unitary matrix used in the
!>          reduction to tridiagonal form.
!>          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
!>          original Hermitian matrix;
!>          if COMPZ = 'I', the orthonormal eigenvectors of the
!>          tridiagonal matrix.
!>          If INFO > 0 on exit, Z contains the eigenvectors associated
!>          with only the stored eigenvalues.
!>          If  COMPZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          COMPZ = 'V' or 'I', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if INFO = i, and i is:
!>                <= N  the Cholesky factorization of the matrix could
!>                      not be performed because the i-th principal minor
!>                      was not positive definite.
!>                > N   the SVD algorithm failed to converge;
!>                      if INFO = N+i, i off-diagonal elements of the
!>                      bidiagonal factor did not converge to zero.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 144 of file zpteqr.f.

145*
146* -- LAPACK computational routine --
147* -- LAPACK is a software package provided by Univ. of Tennessee, --
148* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149*
150* .. Scalar Arguments ..
151 CHARACTER COMPZ
152 INTEGER INFO, LDZ, N
153* ..
154* .. Array Arguments ..
155 DOUBLE PRECISION D( * ), E( * ), WORK( * )
156 COMPLEX*16 Z( LDZ, * )
157* ..
158*
159* ====================================================================
160*
161* .. Parameters ..
162 COMPLEX*16 CZERO, CONE
163 parameter( czero = ( 0.0d+0, 0.0d+0 ),
164 $ cone = ( 1.0d+0, 0.0d+0 ) )
165* ..
166* .. External Functions ..
167 LOGICAL LSAME
168 EXTERNAL lsame
169* ..
170* .. External Subroutines ..
171 EXTERNAL dpttrf, xerbla, zbdsqr, zlaset
172* ..
173* .. Local Arrays ..
174 COMPLEX*16 C( 1, 1 ), VT( 1, 1 )
175* ..
176* .. Local Scalars ..
177 INTEGER I, ICOMPZ, NRU
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC max, sqrt
181* ..
182* .. Executable Statements ..
183*
184* Test the input parameters.
185*
186 info = 0
187*
188 IF( lsame( compz, 'N' ) ) THEN
189 icompz = 0
190 ELSE IF( lsame( compz, 'V' ) ) THEN
191 icompz = 1
192 ELSE IF( lsame( compz, 'I' ) ) THEN
193 icompz = 2
194 ELSE
195 icompz = -1
196 END IF
197 IF( icompz.LT.0 ) THEN
198 info = -1
199 ELSE IF( n.LT.0 ) THEN
200 info = -2
201 ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
202 $ n ) ) ) THEN
203 info = -6
204 END IF
205 IF( info.NE.0 ) THEN
206 CALL xerbla( 'ZPTEQR', -info )
207 RETURN
208 END IF
209*
210* Quick return if possible
211*
212 IF( n.EQ.0 )
213 $ RETURN
214*
215 IF( n.EQ.1 ) THEN
216 IF( icompz.GT.0 )
217 $ z( 1, 1 ) = cone
218 RETURN
219 END IF
220 IF( icompz.EQ.2 )
221 $ CALL zlaset( 'Full', n, n, czero, cone, z, ldz )
222*
223* Call DPTTRF to factor the matrix.
224*
225 CALL dpttrf( n, d, e, info )
226 IF( info.NE.0 )
227 $ RETURN
228 DO 10 i = 1, n
229 d( i ) = sqrt( d( i ) )
230 10 CONTINUE
231 DO 20 i = 1, n - 1
232 e( i ) = e( i )*d( i )
233 20 CONTINUE
234*
235* Call ZBDSQR to compute the singular values/vectors of the
236* bidiagonal factor.
237*
238 IF( icompz.GT.0 ) THEN
239 nru = n
240 ELSE
241 nru = 0
242 END IF
243 CALL zbdsqr( 'Lower', n, 0, nru, 0, d, e, vt, 1, z, ldz, c, 1,
244 $ work, info )
245*
246* Square the singular values.
247*
248 IF( info.EQ.0 ) THEN
249 DO 30 i = 1, n
250 d( i ) = d( i )*d( i )
251 30 CONTINUE
252 ELSE
253 info = n + info
254 END IF
255*
256 RETURN
257*
258* End of ZPTEQR
259*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
zbdsqr
subroutine zbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
ZBDSQR
Definition zbdsqr.f:222
dpttrf
subroutine dpttrf(n, d, e, info)
DPTTRF
Definition dpttrf.f:91
max
#define max(a, b)
Definition macros.h:21

◆ zptrfs()

subroutine zptrfs ( character uplo,
integer n,
integer nrhs,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
double precision, dimension( * ) df,
complex*16, dimension( * ) ef,
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 )

ZPTRFS

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

Purpose:
!>
!> ZPTRFS improves the computed solution to a system of linear
!> equations when the coefficient matrix is Hermitian positive definite
!> and tridiagonal, and provides error bounds and backward error
!> estimates for the solution.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          Specifies whether the superdiagonal or the subdiagonal of the
!>          tridiagonal matrix A is stored and the form of the
!>          factorization:
!>          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
!>          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
!>          (The two forms are equivalent if A is real.)
!>
[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]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n real diagonal elements of the tridiagonal matrix A.
!>
[in]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the tridiagonal matrix A
!>          (see UPLO).
!>
[in]DF
!>          DF is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from
!>          the factorization computed by ZPTTRF.
!>
[in]EF
!>          EF is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the unit bidiagonal
!>          factor U or L from the factorization computed by ZPTTRF
!>          (see UPLO).
!>
[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 ZPTTRS.
!>          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 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).
!>
[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 (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 181 of file zptrfs.f.

183*
184* -- LAPACK computational routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER UPLO
190 INTEGER INFO, LDB, LDX, N, NRHS
191* ..
192* .. Array Arguments ..
193 DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
194 $ RWORK( * )
195 COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
196 $ X( LDX, * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 INTEGER ITMAX
203 parameter( itmax = 5 )
204 DOUBLE PRECISION ZERO
205 parameter( zero = 0.0d+0 )
206 DOUBLE PRECISION ONE
207 parameter( one = 1.0d+0 )
208 DOUBLE PRECISION TWO
209 parameter( two = 2.0d+0 )
210 DOUBLE PRECISION THREE
211 parameter( three = 3.0d+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL UPPER
215 INTEGER COUNT, I, IX, J, NZ
216 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
217 COMPLEX*16 BI, CX, DX, EX, ZDUM
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 INTEGER IDAMAX
222 DOUBLE PRECISION DLAMCH
223 EXTERNAL lsame, idamax, dlamch
224* ..
225* .. External Subroutines ..
226 EXTERNAL xerbla, zaxpy, zpttrs
227* ..
228* .. Intrinsic Functions ..
229 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max
230* ..
231* .. Statement Functions ..
232 DOUBLE PRECISION CABS1
233* ..
234* .. Statement Function definitions ..
235 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
236* ..
237* .. Executable Statements ..
238*
239* Test the input parameters.
240*
241 info = 0
242 upper = lsame( uplo, 'U' )
243 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
244 info = -1
245 ELSE IF( n.LT.0 ) THEN
246 info = -2
247 ELSE IF( nrhs.LT.0 ) THEN
248 info = -3
249 ELSE IF( ldb.LT.max( 1, n ) ) THEN
250 info = -9
251 ELSE IF( ldx.LT.max( 1, n ) ) THEN
252 info = -11
253 END IF
254 IF( info.NE.0 ) THEN
255 CALL xerbla( 'ZPTRFS', -info )
256 RETURN
257 END IF
258*
259* Quick return if possible
260*
261 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
262 DO 10 j = 1, nrhs
263 ferr( j ) = zero
264 berr( j ) = zero
265 10 CONTINUE
266 RETURN
267 END IF
268*
269* NZ = maximum number of nonzero elements in each row of A, plus 1
270*
271 nz = 4
272 eps = dlamch( 'Epsilon' )
273 safmin = dlamch( 'Safe minimum' )
274 safe1 = nz*safmin
275 safe2 = safe1 / eps
276*
277* Do for each right hand side
278*
279 DO 100 j = 1, nrhs
280*
281 count = 1
282 lstres = three
283 20 CONTINUE
284*
285* Loop until stopping criterion is satisfied.
286*
287* Compute residual R = B - A * X. Also compute
288* abs(A)*abs(x) + abs(b) for use in the backward error bound.
289*
290 IF( upper ) THEN
291 IF( n.EQ.1 ) THEN
292 bi = b( 1, j )
293 dx = d( 1 )*x( 1, j )
294 work( 1 ) = bi - dx
295 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
296 ELSE
297 bi = b( 1, j )
298 dx = d( 1 )*x( 1, j )
299 ex = e( 1 )*x( 2, j )
300 work( 1 ) = bi - dx - ex
301 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
302 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
303 DO 30 i = 2, n - 1
304 bi = b( i, j )
305 cx = dconjg( e( i-1 ) )*x( i-1, j )
306 dx = d( i )*x( i, j )
307 ex = e( i )*x( i+1, j )
308 work( i ) = bi - cx - dx - ex
309 rwork( i ) = cabs1( bi ) +
310 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
311 $ cabs1( dx ) + cabs1( e( i ) )*
312 $ cabs1( x( i+1, j ) )
313 30 CONTINUE
314 bi = b( n, j )
315 cx = dconjg( e( n-1 ) )*x( n-1, j )
316 dx = d( n )*x( n, j )
317 work( n ) = bi - cx - dx
318 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
319 $ cabs1( x( n-1, j ) ) + cabs1( dx )
320 END IF
321 ELSE
322 IF( n.EQ.1 ) THEN
323 bi = b( 1, j )
324 dx = d( 1 )*x( 1, j )
325 work( 1 ) = bi - dx
326 rwork( 1 ) = cabs1( bi ) + cabs1( dx )
327 ELSE
328 bi = b( 1, j )
329 dx = d( 1 )*x( 1, j )
330 ex = dconjg( e( 1 ) )*x( 2, j )
331 work( 1 ) = bi - dx - ex
332 rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
333 $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
334 DO 40 i = 2, n - 1
335 bi = b( i, j )
336 cx = e( i-1 )*x( i-1, j )
337 dx = d( i )*x( i, j )
338 ex = dconjg( e( i ) )*x( i+1, j )
339 work( i ) = bi - cx - dx - ex
340 rwork( i ) = cabs1( bi ) +
341 $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
342 $ cabs1( dx ) + cabs1( e( i ) )*
343 $ cabs1( x( i+1, j ) )
344 40 CONTINUE
345 bi = b( n, j )
346 cx = e( n-1 )*x( n-1, j )
347 dx = d( n )*x( n, j )
348 work( n ) = bi - cx - dx
349 rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
350 $ cabs1( x( n-1, j ) ) + cabs1( dx )
351 END IF
352 END IF
353*
354* Compute componentwise relative backward error from formula
355*
356* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
357*
358* where abs(Z) is the componentwise absolute value of the matrix
359* or vector Z. If the i-th component of the denominator is less
360* than SAFE2, then SAFE1 is added to the i-th components of the
361* numerator and denominator before dividing.
362*
363 s = zero
364 DO 50 i = 1, n
365 IF( rwork( i ).GT.safe2 ) THEN
366 s = max( s, cabs1( work( i ) ) / rwork( i ) )
367 ELSE
368 s = max( s, ( cabs1( work( i ) )+safe1 ) /
369 $ ( rwork( i )+safe1 ) )
370 END IF
371 50 CONTINUE
372 berr( j ) = s
373*
374* Test stopping criterion. Continue iterating if
375* 1) The residual BERR(J) is larger than machine epsilon, and
376* 2) BERR(J) decreased by at least a factor of 2 during the
377* last iteration, and
378* 3) At most ITMAX iterations tried.
379*
380 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
381 $ count.LE.itmax ) THEN
382*
383* Update solution and try again.
384*
385 CALL zpttrs( uplo, n, 1, df, ef, work, n, info )
386 CALL zaxpy( n, dcmplx( one ), work, 1, x( 1, j ), 1 )
387 lstres = berr( j )
388 count = count + 1
389 GO TO 20
390 END IF
391*
392* Bound error from formula
393*
394* norm(X - XTRUE) / norm(X) .le. FERR =
395* norm( abs(inv(A))*
396* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
397*
398* where
399* norm(Z) is the magnitude of the largest component of Z
400* inv(A) is the inverse of A
401* abs(Z) is the componentwise absolute value of the matrix or
402* vector Z
403* NZ is the maximum number of nonzeros in any row of A, plus 1
404* EPS is machine epsilon
405*
406* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
407* is incremented by SAFE1 if the i-th component of
408* abs(A)*abs(X) + abs(B) is less than SAFE2.
409*
410 DO 60 i = 1, n
411 IF( rwork( i ).GT.safe2 ) THEN
412 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
413 ELSE
414 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
415 $ safe1
416 END IF
417 60 CONTINUE
418 ix = idamax( n, rwork, 1 )
419 ferr( j ) = rwork( ix )
420*
421* Estimate the norm of inv(A).
422*
423* Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
424*
425* m(i,j) = abs(A(i,j)), i = j,
426* m(i,j) = -abs(A(i,j)), i .ne. j,
427*
428* and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
429*
430* Solve M(L) * x = e.
431*
432 rwork( 1 ) = one
433 DO 70 i = 2, n
434 rwork( i ) = one + rwork( i-1 )*abs( ef( i-1 ) )
435 70 CONTINUE
436*
437* Solve D * M(L)**H * x = b.
438*
439 rwork( n ) = rwork( n ) / df( n )
440 DO 80 i = n - 1, 1, -1
441 rwork( i ) = rwork( i ) / df( i ) +
442 $ rwork( i+1 )*abs( ef( i ) )
443 80 CONTINUE
444*
445* Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
446*
447 ix = idamax( n, rwork, 1 )
448 ferr( j ) = ferr( j )*abs( rwork( ix ) )
449*
450* Normalize error.
451*
452 lstres = zero
453 DO 90 i = 1, n
454 lstres = max( lstres, abs( x( i, j ) ) )
455 90 CONTINUE
456 IF( lstres.NE.zero )
457 $ ferr( j ) = ferr( j ) / lstres
458*
459 100 CONTINUE
460*
461 RETURN
462*
463* End of ZPTRFS
464*
zpttrs
subroutine zpttrs(uplo, n, nrhs, d, e, b, ldb, info)
ZPTTRS
Definition zpttrs.f:121
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

◆ zpttrf()

subroutine zpttrf ( integer n,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
integer info )

ZPTTRF

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

Purpose:
!>
!> ZPTTRF computes the L*D*L**H factorization of a complex Hermitian
!> positive definite tridiagonal matrix A.  The factorization may also
!> be regarded as having the form A = U**H *D*U.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, the n diagonal elements of the tridiagonal matrix
!>          A.  On exit, the n diagonal elements of the diagonal matrix
!>          D from the L*D*L**H factorization of A.
!>
[in,out]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          On entry, the (n-1) subdiagonal elements of the tridiagonal
!>          matrix A.  On exit, the (n-1) subdiagonal elements of the
!>          unit bidiagonal factor L from the L*D*L**H factorization of A.
!>          E can also be regarded as the superdiagonal of the unit
!>          bidiagonal factor U from the U**H *D*U factorization of A.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!>          > 0: if INFO = k, the leading minor of order k is not
!>               positive definite; if k < N, the factorization could not
!>               be completed, while if k = N, the factorization was
!>               completed, but D(N) <= 0.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 91 of file zpttrf.f.

92*
93* -- LAPACK computational routine --
94* -- LAPACK is a software package provided by Univ. of Tennessee, --
95* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96*
97* .. Scalar Arguments ..
98 INTEGER INFO, N
99* ..
100* .. Array Arguments ..
101 DOUBLE PRECISION D( * )
102 COMPLEX*16 E( * )
103* ..
104*
105* =====================================================================
106*
107* .. Parameters ..
108 DOUBLE PRECISION ZERO
109 parameter( zero = 0.0d+0 )
110* ..
111* .. Local Scalars ..
112 INTEGER I, I4
113 DOUBLE PRECISION EII, EIR, F, G
114* ..
115* .. External Subroutines ..
116 EXTERNAL xerbla
117* ..
118* .. Intrinsic Functions ..
119 INTRINSIC dble, dcmplx, dimag, mod
120* ..
121* .. Executable Statements ..
122*
123* Test the input parameters.
124*
125 info = 0
126 IF( n.LT.0 ) THEN
127 info = -1
128 CALL xerbla( 'ZPTTRF', -info )
129 RETURN
130 END IF
131*
132* Quick return if possible
133*
134 IF( n.EQ.0 )
135 $ RETURN
136*
137* Compute the L*D*L**H (or U**H *D*U) factorization of A.
138*
139 i4 = mod( n-1, 4 )
140 DO 10 i = 1, i4
141 IF( d( i ).LE.zero ) THEN
142 info = i
143 GO TO 30
144 END IF
145 eir = dble( e( i ) )
146 eii = dimag( e( i ) )
147 f = eir / d( i )
148 g = eii / d( i )
149 e( i ) = dcmplx( f, g )
150 d( i+1 ) = d( i+1 ) - f*eir - g*eii
151 10 CONTINUE
152*
153 DO 20 i = i4 + 1, n - 4, 4
154*
155* Drop out of the loop if d(i) <= 0: the matrix is not positive
156* definite.
157*
158 IF( d( i ).LE.zero ) THEN
159 info = i
160 GO TO 30
161 END IF
162*
163* Solve for e(i) and d(i+1).
164*
165 eir = dble( e( i ) )
166 eii = dimag( e( i ) )
167 f = eir / d( i )
168 g = eii / d( i )
169 e( i ) = dcmplx( f, g )
170 d( i+1 ) = d( i+1 ) - f*eir - g*eii
171*
172 IF( d( i+1 ).LE.zero ) THEN
173 info = i + 1
174 GO TO 30
175 END IF
176*
177* Solve for e(i+1) and d(i+2).
178*
179 eir = dble( e( i+1 ) )
180 eii = dimag( e( i+1 ) )
181 f = eir / d( i+1 )
182 g = eii / d( i+1 )
183 e( i+1 ) = dcmplx( f, g )
184 d( i+2 ) = d( i+2 ) - f*eir - g*eii
185*
186 IF( d( i+2 ).LE.zero ) THEN
187 info = i + 2
188 GO TO 30
189 END IF
190*
191* Solve for e(i+2) and d(i+3).
192*
193 eir = dble( e( i+2 ) )
194 eii = dimag( e( i+2 ) )
195 f = eir / d( i+2 )
196 g = eii / d( i+2 )
197 e( i+2 ) = dcmplx( f, g )
198 d( i+3 ) = d( i+3 ) - f*eir - g*eii
199*
200 IF( d( i+3 ).LE.zero ) THEN
201 info = i + 3
202 GO TO 30
203 END IF
204*
205* Solve for e(i+3) and d(i+4).
206*
207 eir = dble( e( i+3 ) )
208 eii = dimag( e( i+3 ) )
209 f = eir / d( i+3 )
210 g = eii / d( i+3 )
211 e( i+3 ) = dcmplx( f, g )
212 d( i+4 ) = d( i+4 ) - f*eir - g*eii
213 20 CONTINUE
214*
215* Check d(n) for positive definiteness.
216*
217 IF( d( n ).LE.zero )
218 $ info = n
219*
220 30 CONTINUE
221 RETURN
222*
223* End of ZPTTRF
224*

◆ zpttrs()

subroutine zpttrs ( character uplo,
integer n,
integer nrhs,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZPTTRS

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

Purpose:
!>
!> ZPTTRS solves a tridiagonal system of the form
!>    A * X = B
!> using the factorization A = U**H *D* U or A = L*D*L**H computed by ZPTTRF.
!> D is a diagonal matrix specified in the vector D, U (or L) is a unit
!> bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
!> the vector E, and X and B are N by NRHS matrices.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          Specifies the form of the factorization and whether the
!>          vector E is the superdiagonal of the upper bidiagonal factor
!>          U or the subdiagonal of the lower bidiagonal factor L.
!>          = 'U':  A = U**H *D*U, E is the superdiagonal of U
!>          = 'L':  A = L*D*L**H, E is the subdiagonal of L
!>
[in]N
!>          N is INTEGER
!>          The order of the tridiagonal 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]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from the
!>          factorization A = U**H *D*U or A = L*D*L**H.
!>
[in]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          If UPLO = 'U', the (n-1) superdiagonal elements of the unit
!>          bidiagonal factor U from the factorization A = U**H*D*U.
!>          If UPLO = 'L', the (n-1) subdiagonal elements of the unit
!>          bidiagonal factor L from the factorization A = L*D*L**H.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the right hand side vectors B for the system of
!>          linear equations.
!>          On exit, 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 120 of file zpttrs.f.

121*
122* -- LAPACK computational routine --
123* -- LAPACK is a software package provided by Univ. of Tennessee, --
124* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
125*
126* .. Scalar Arguments ..
127 CHARACTER UPLO
128 INTEGER INFO, LDB, N, NRHS
129* ..
130* .. Array Arguments ..
131 DOUBLE PRECISION D( * )
132 COMPLEX*16 B( LDB, * ), E( * )
133* ..
134*
135* =====================================================================
136*
137* .. Local Scalars ..
138 LOGICAL UPPER
139 INTEGER IUPLO, J, JB, NB
140* ..
141* .. External Functions ..
142 INTEGER ILAENV
143 EXTERNAL ilaenv
144* ..
145* .. External Subroutines ..
146 EXTERNAL xerbla, zptts2
147* ..
148* .. Intrinsic Functions ..
149 INTRINSIC max, min
150* ..
151* .. Executable Statements ..
152*
153* Test the input arguments.
154*
155 info = 0
156 upper = ( uplo.EQ.'U' .OR. uplo.EQ.'u' )
157 IF( .NOT.upper .AND. .NOT.( uplo.EQ.'L' .OR. uplo.EQ.'l' ) ) THEN
158 info = -1
159 ELSE IF( n.LT.0 ) THEN
160 info = -2
161 ELSE IF( nrhs.LT.0 ) THEN
162 info = -3
163 ELSE IF( ldb.LT.max( 1, n ) ) THEN
164 info = -7
165 END IF
166 IF( info.NE.0 ) THEN
167 CALL xerbla( 'ZPTTRS', -info )
168 RETURN
169 END IF
170*
171* Quick return if possible
172*
173 IF( n.EQ.0 .OR. nrhs.EQ.0 )
174 $ RETURN
175*
176* Determine the number of right-hand sides to solve at a time.
177*
178 IF( nrhs.EQ.1 ) THEN
179 nb = 1
180 ELSE
181 nb = max( 1, ilaenv( 1, 'ZPTTRS', uplo, n, nrhs, -1, -1 ) )
182 END IF
183*
184* Decode UPLO
185*
186 IF( upper ) THEN
187 iuplo = 1
188 ELSE
189 iuplo = 0
190 END IF
191*
192 IF( nb.GE.nrhs ) THEN
193 CALL zptts2( iuplo, n, nrhs, d, e, b, ldb )
194 ELSE
195 DO 10 j = 1, nrhs, nb
196 jb = min( nrhs-j+1, nb )
197 CALL zptts2( iuplo, n, jb, d, e, b( 1, j ), ldb )
198 10 CONTINUE
199 END IF
200*
201 RETURN
202*
203* End of ZPTTRS
204*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
zptts2
subroutine zptts2(iuplo, n, nrhs, d, e, b, ldb)
ZPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf...
Definition zptts2.f:113
min
#define min(a, b)
Definition macros.h:20

◆ zptts2()

subroutine zptts2 ( integer iuplo,
integer n,
integer nrhs,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
complex*16, dimension( ldb, * ) b,
integer ldb )

ZPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf.

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

Purpose:
!>
!> ZPTTS2 solves a tridiagonal system of the form
!>    A * X = B
!> using the factorization A = U**H *D*U or A = L*D*L**H computed by ZPTTRF.
!> D is a diagonal matrix specified in the vector D, U (or L) is a unit
!> bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
!> the vector E, and X and B are N by NRHS matrices.
!>
Parameters
[in]IUPLO
!>          IUPLO is INTEGER
!>          Specifies the form of the factorization and whether the
!>          vector E is the superdiagonal of the upper bidiagonal factor
!>          U or the subdiagonal of the lower bidiagonal factor L.
!>          = 1:  A = U**H *D*U, E is the superdiagonal of U
!>          = 0:  A = L*D*L**H, E is the subdiagonal of L
!>
[in]N
!>          N is INTEGER
!>          The order of the tridiagonal 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]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from the
!>          factorization A = U**H *D*U or A = L*D*L**H.
!>
[in]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          If IUPLO = 1, the (n-1) superdiagonal elements of the unit
!>          bidiagonal factor U from the factorization A = U**H*D*U.
!>          If IUPLO = 0, the (n-1) subdiagonal elements of the unit
!>          bidiagonal factor L from the factorization A = L*D*L**H.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the right hand side vectors B for the system of
!>          linear equations.
!>          On exit, 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 112 of file zptts2.f.

113*
114* -- LAPACK computational routine --
115* -- LAPACK is a software package provided by Univ. of Tennessee, --
116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117*
118* .. Scalar Arguments ..
119 INTEGER IUPLO, LDB, N, NRHS
120* ..
121* .. Array Arguments ..
122 DOUBLE PRECISION D( * )
123 COMPLEX*16 B( LDB, * ), E( * )
124* ..
125*
126* =====================================================================
127*
128* .. Local Scalars ..
129 INTEGER I, J
130* ..
131* .. External Subroutines ..
132 EXTERNAL zdscal
133* ..
134* .. Intrinsic Functions ..
135 INTRINSIC dconjg
136* ..
137* .. Executable Statements ..
138*
139* Quick return if possible
140*
141 IF( n.LE.1 ) THEN
142 IF( n.EQ.1 )
143 $ CALL zdscal( nrhs, 1.d0 / d( 1 ), b, ldb )
144 RETURN
145 END IF
146*
147 IF( iuplo.EQ.1 ) THEN
148*
149* Solve A * X = B using the factorization A = U**H *D*U,
150* overwriting each right hand side vector with its solution.
151*
152 IF( nrhs.LE.2 ) THEN
153 j = 1
154 10 CONTINUE
155*
156* Solve U**H * x = b.
157*
158 DO 20 i = 2, n
159 b( i, j ) = b( i, j ) - b( i-1, j )*dconjg( e( i-1 ) )
160 20 CONTINUE
161*
162* Solve D * U * x = b.
163*
164 DO 30 i = 1, n
165 b( i, j ) = b( i, j ) / d( i )
166 30 CONTINUE
167 DO 40 i = n - 1, 1, -1
168 b( i, j ) = b( i, j ) - b( i+1, j )*e( i )
169 40 CONTINUE
170 IF( j.LT.nrhs ) THEN
171 j = j + 1
172 GO TO 10
173 END IF
174 ELSE
175 DO 70 j = 1, nrhs
176*
177* Solve U**H * x = b.
178*
179 DO 50 i = 2, n
180 b( i, j ) = b( i, j ) - b( i-1, j )*dconjg( e( i-1 ) )
181 50 CONTINUE
182*
183* Solve D * U * x = b.
184*
185 b( n, j ) = b( n, j ) / d( n )
186 DO 60 i = n - 1, 1, -1
187 b( i, j ) = b( i, j ) / d( i ) - b( i+1, j )*e( i )
188 60 CONTINUE
189 70 CONTINUE
190 END IF
191 ELSE
192*
193* Solve A * X = B using the factorization A = L*D*L**H,
194* overwriting each right hand side vector with its solution.
195*
196 IF( nrhs.LE.2 ) THEN
197 j = 1
198 80 CONTINUE
199*
200* Solve L * x = b.
201*
202 DO 90 i = 2, n
203 b( i, j ) = b( i, j ) - b( i-1, j )*e( i-1 )
204 90 CONTINUE
205*
206* Solve D * L**H * x = b.
207*
208 DO 100 i = 1, n
209 b( i, j ) = b( i, j ) / d( i )
210 100 CONTINUE
211 DO 110 i = n - 1, 1, -1
212 b( i, j ) = b( i, j ) - b( i+1, j )*dconjg( e( i ) )
213 110 CONTINUE
214 IF( j.LT.nrhs ) THEN
215 j = j + 1
216 GO TO 80
217 END IF
218 ELSE
219 DO 140 j = 1, nrhs
220*
221* Solve L * x = b.
222*
223 DO 120 i = 2, n
224 b( i, j ) = b( i, j ) - b( i-1, j )*e( i-1 )
225 120 CONTINUE
226*
227* Solve D * L**H * x = b.
228*
229 b( n, j ) = b( n, j ) / d( n )
230 DO 130 i = n - 1, 1, -1
231 b( i, j ) = b( i, j ) / d( i ) -
232 $ b( i+1, j )*dconjg( e( i ) )
233 130 CONTINUE
234 140 CONTINUE
235 END IF
236 END IF
237*
238 RETURN
239*
240* End of ZPTTS2
241*
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78