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

Functions

subroutine cptcon (n, d, e, anorm, rcond, rwork, info)
 CPTCON
subroutine cpteqr (compz, n, d, e, z, ldz, work, info)
 CPTEQR
subroutine cptrfs (uplo, n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 CPTRFS
subroutine cpttrf (n, d, e, info)
 CPTTRF
subroutine cpttrs (uplo, n, nrhs, d, e, b, ldb, info)
 CPTTRS
subroutine cptts2 (iuplo, n, nrhs, d, e, b, ldb)
 CPTTS2 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 complex computational functions for PT matrices

Function Documentation

◆ cptcon()

subroutine cptcon ( integer n,
real, dimension( * ) d,
complex, dimension( * ) e,
real anorm,
real rcond,
real, dimension( * ) rwork,
integer info )

CPTCON

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

Purpose:
!>
!> CPTCON 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
!> CPTTRF.
!>
!> 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 REAL array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from the
!>          factorization of A, as computed by CPTTRF.
!>
[in]E
!>          E is COMPLEX 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 CPTTRF.
!>
[in]ANORM
!>          ANORM is REAL
!>          The 1-norm of the original matrix A.
!>
[out]RCOND
!>          RCOND is REAL
!>          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 REAL 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 cptcon.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 REAL ANORM, RCOND
127* ..
128* .. Array Arguments ..
129 REAL D( * ), RWORK( * )
130 COMPLEX E( * )
131* ..
132*
133* =====================================================================
134*
135* .. Parameters ..
136 REAL ONE, ZERO
137 parameter( one = 1.0e+0, zero = 0.0e+0 )
138* ..
139* .. Local Scalars ..
140 INTEGER I, IX
141 REAL AINVNM
142* ..
143* .. External Functions ..
144 INTEGER ISAMAX
145 EXTERNAL isamax
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( 'CPTCON', -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 = isamax( 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 CPTCON
219*
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71

◆ cpteqr()

subroutine cpteqr ( character compz,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
complex, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer info )

CPTEQR

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

Purpose:
!>
!> CPTEQR computes all eigenvalues and, optionally, eigenvectors of a
!> symmetric positive definite tridiagonal matrix by first factoring the
!> matrix using SPTTRF and then calling CBDSQR 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 CHETRD, CHPTRD, or CHBTRD 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 REAL 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 REAL 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 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 REAL 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 cpteqr.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 REAL D( * ), E( * ), WORK( * )
156 COMPLEX Z( LDZ, * )
157* ..
158*
159* ====================================================================
160*
161* .. Parameters ..
162 COMPLEX CZERO, CONE
163 parameter( czero = ( 0.0e+0, 0.0e+0 ),
164 $ cone = ( 1.0e+0, 0.0e+0 ) )
165* ..
166* .. External Functions ..
167 LOGICAL LSAME
168 EXTERNAL lsame
169* ..
170* .. External Subroutines ..
171 EXTERNAL cbdsqr, claset, spttrf, xerbla
172* ..
173* .. Local Arrays ..
174 COMPLEX 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( 'CPTEQR', -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 claset( 'Full', n, n, czero, cone, z, ldz )
222*
223* Call SPTTRF to factor the matrix.
224*
225 CALL spttrf( 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 CBDSQR 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 cbdsqr( '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 CPTEQR
259*
spttrf
subroutine spttrf(n, d, e, info)
SPTTRF
Definition spttrf.f:91
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
cbdsqr
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
Definition cbdsqr.f:222
max
#define max(a, b)
Definition macros.h:21

◆ cptrfs()

subroutine cptrfs ( character uplo,
integer n,
integer nrhs,
real, dimension( * ) d,
complex, dimension( * ) e,
real, dimension( * ) df,
complex, dimension( * ) ef,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( ldx, * ) x,
integer ldx,
real, dimension( * ) ferr,
real, dimension( * ) berr,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CPTRFS

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

Purpose:
!>
!> CPTRFS 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 REAL array, dimension (N)
!>          The n real diagonal elements of the tridiagonal matrix A.
!>
[in]E
!>          E is COMPLEX array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the tridiagonal matrix A
!>          (see UPLO).
!>
[in]DF
!>          DF is REAL array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from
!>          the factorization computed by CPTTRF.
!>
[in]EF
!>          EF is COMPLEX array, dimension (N-1)
!>          The (n-1) off-diagonal elements of the unit bidiagonal
!>          factor U or L from the factorization computed by CPTTRF
!>          (see UPLO).
!>
[in]B
!>          B is COMPLEX 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 array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by CPTTRS.
!>          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 REAL 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 REAL 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 array, dimension (N)
!>
[out]RWORK
!>          RWORK is REAL 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 cptrfs.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 REAL BERR( * ), D( * ), DF( * ), FERR( * ),
194 $ RWORK( * )
195 COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
196 $ X( LDX, * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 INTEGER ITMAX
203 parameter( itmax = 5 )
204 REAL ZERO
205 parameter( zero = 0.0e+0 )
206 REAL ONE
207 parameter( one = 1.0e+0 )
208 REAL TWO
209 parameter( two = 2.0e+0 )
210 REAL THREE
211 parameter( three = 3.0e+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL UPPER
215 INTEGER COUNT, I, IX, J, NZ
216 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
217 COMPLEX BI, CX, DX, EX, ZDUM
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 INTEGER ISAMAX
222 REAL SLAMCH
223 EXTERNAL lsame, isamax, slamch
224* ..
225* .. External Subroutines ..
226 EXTERNAL caxpy, cpttrs, xerbla
227* ..
228* .. Intrinsic Functions ..
229 INTRINSIC abs, aimag, cmplx, conjg, max, real
230* ..
231* .. Statement Functions ..
232 REAL CABS1
233* ..
234* .. Statement Function definitions ..
235 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CPTRFS', -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 = slamch( 'Epsilon' )
273 safmin = slamch( '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 = conjg( 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 = conjg( 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 = conjg( 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 = conjg( 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 cpttrs( uplo, n, 1, df, ef, work, n, info )
386 CALL caxpy( n, cmplx( 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 = isamax( 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 = isamax( 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 CPTRFS
464*
cmplx
float cmplx[2]
Definition pblas.h:136
cpttrs
subroutine cpttrs(uplo, n, nrhs, d, e, b, ldb, info)
CPTTRS
Definition cpttrs.f:121
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ cpttrf()

subroutine cpttrf ( integer n,
real, dimension( * ) d,
complex, dimension( * ) e,
integer info )

CPTTRF

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

Purpose:
!>
!> CPTTRF 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 REAL 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 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 cpttrf.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 REAL D( * )
102 COMPLEX E( * )
103* ..
104*
105* =====================================================================
106*
107* .. Parameters ..
108 REAL ZERO
109 parameter( zero = 0.0e+0 )
110* ..
111* .. Local Scalars ..
112 INTEGER I, I4
113 REAL EII, EIR, F, G
114* ..
115* .. External Subroutines ..
116 EXTERNAL xerbla
117* ..
118* .. Intrinsic Functions ..
119 INTRINSIC aimag, cmplx, mod, real
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( 'CPTTRF', -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 20
144 END IF
145 eir = real( e( i ) )
146 eii = aimag( e( i ) )
147 f = eir / d( i )
148 g = eii / d( i )
149 e( i ) = cmplx( f, g )
150 d( i+1 ) = d( i+1 ) - f*eir - g*eii
151 10 CONTINUE
152*
153 DO 110 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 20
161 END IF
162*
163* Solve for e(i) and d(i+1).
164*
165 eir = real( e( i ) )
166 eii = aimag( e( i ) )
167 f = eir / d( i )
168 g = eii / d( i )
169 e( i ) = cmplx( 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 20
175 END IF
176*
177* Solve for e(i+1) and d(i+2).
178*
179 eir = real( e( i+1 ) )
180 eii = aimag( e( i+1 ) )
181 f = eir / d( i+1 )
182 g = eii / d( i+1 )
183 e( i+1 ) = cmplx( 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 20
189 END IF
190*
191* Solve for e(i+2) and d(i+3).
192*
193 eir = real( e( i+2 ) )
194 eii = aimag( e( i+2 ) )
195 f = eir / d( i+2 )
196 g = eii / d( i+2 )
197 e( i+2 ) = cmplx( 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 20
203 END IF
204*
205* Solve for e(i+3) and d(i+4).
206*
207 eir = real( e( i+3 ) )
208 eii = aimag( e( i+3 ) )
209 f = eir / d( i+3 )
210 g = eii / d( i+3 )
211 e( i+3 ) = cmplx( f, g )
212 d( i+4 ) = d( i+4 ) - f*eir - g*eii
213 110 CONTINUE
214*
215* Check d(n) for positive definiteness.
216*
217 IF( d( n ).LE.zero )
218 $ info = n
219*
220 20 CONTINUE
221 RETURN
222*
223* End of CPTTRF
224*

◆ cpttrs()

subroutine cpttrs ( character uplo,
integer n,
integer nrhs,
real, dimension( * ) d,
complex, dimension( * ) e,
complex, dimension( ldb, * ) b,
integer ldb,
integer info )

CPTTRS

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

Purpose:
!>
!> CPTTRS 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 CPTTRF.
!> 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 REAL 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 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 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 cpttrs.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 REAL D( * )
132 COMPLEX 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 cptts2, xerbla
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( 'CPTTRS', -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, 'CPTTRS', 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 cptts2( 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 cptts2( iuplo, n, jb, d, e, b( 1, j ), ldb )
198 10 CONTINUE
199 END IF
200*
201 RETURN
202*
203* End of CPTTRS
204*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
cptts2
subroutine cptts2(iuplo, n, nrhs, d, e, b, ldb)
CPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf...
Definition cptts2.f:113
min
#define min(a, b)
Definition macros.h:20

◆ cptts2()

subroutine cptts2 ( integer iuplo,
integer n,
integer nrhs,
real, dimension( * ) d,
complex, dimension( * ) e,
complex, dimension( ldb, * ) b,
integer ldb )

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

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

Purpose:
!>
!> CPTTS2 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 CPTTRF.
!> 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 REAL 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 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 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 cptts2.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 REAL D( * )
123 COMPLEX B( LDB, * ), E( * )
124* ..
125*
126* =====================================================================
127*
128* .. Local Scalars ..
129 INTEGER I, J
130* ..
131* .. External Subroutines ..
132 EXTERNAL csscal
133* ..
134* .. Intrinsic Functions ..
135 INTRINSIC conjg
136* ..
137* .. Executable Statements ..
138*
139* Quick return if possible
140*
141 IF( n.LE.1 ) THEN
142 IF( n.EQ.1 )
143 $ CALL csscal( nrhs, 1. / 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 5 CONTINUE
155*
156* Solve U**H * x = b.
157*
158 DO 10 i = 2, n
159 b( i, j ) = b( i, j ) - b( i-1, j )*conjg( e( i-1 ) )
160 10 CONTINUE
161*
162* Solve D * U * x = b.
163*
164 DO 20 i = 1, n
165 b( i, j ) = b( i, j ) / d( i )
166 20 CONTINUE
167 DO 30 i = n - 1, 1, -1
168 b( i, j ) = b( i, j ) - b( i+1, j )*e( i )
169 30 CONTINUE
170 IF( j.LT.nrhs ) THEN
171 j = j + 1
172 GO TO 5
173 END IF
174 ELSE
175 DO 60 j = 1, nrhs
176*
177* Solve U**H * x = b.
178*
179 DO 40 i = 2, n
180 b( i, j ) = b( i, j ) - b( i-1, j )*conjg( e( i-1 ) )
181 40 CONTINUE
182*
183* Solve D * U * x = b.
184*
185 b( n, j ) = b( n, j ) / d( n )
186 DO 50 i = n - 1, 1, -1
187 b( i, j ) = b( i, j ) / d( i ) - b( i+1, j )*e( i )
188 50 CONTINUE
189 60 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 65 CONTINUE
199*
200* Solve L * x = b.
201*
202 DO 70 i = 2, n
203 b( i, j ) = b( i, j ) - b( i-1, j )*e( i-1 )
204 70 CONTINUE
205*
206* Solve D * L**H * x = b.
207*
208 DO 80 i = 1, n
209 b( i, j ) = b( i, j ) / d( i )
210 80 CONTINUE
211 DO 90 i = n - 1, 1, -1
212 b( i, j ) = b( i, j ) - b( i+1, j )*conjg( e( i ) )
213 90 CONTINUE
214 IF( j.LT.nrhs ) THEN
215 j = j + 1
216 GO TO 65
217 END IF
218 ELSE
219 DO 120 j = 1, nrhs
220*
221* Solve L * x = b.
222*
223 DO 100 i = 2, n
224 b( i, j ) = b( i, j ) - b( i-1, j )*e( i-1 )
225 100 CONTINUE
226*
227* Solve D * L**H * x = b.
228*
229 b( n, j ) = b( n, j ) / d( n )
230 DO 110 i = n - 1, 1, -1
231 b( i, j ) = b( i, j ) / d( i ) -
232 $ b( i+1, j )*conjg( e( i ) )
233 110 CONTINUE
234 120 CONTINUE
235 END IF
236 END IF
237*
238 RETURN
239*
240* End of CPTTS2
241*
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78