LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for double:

Functions

subroutine dptcon (N, D, E, ANORM, RCOND, WORK, INFO)
 DPTCON More...
 
subroutine dpteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO)
 DPTEQR More...
 
subroutine dptrfs (N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, INFO)
 DPTRFS More...
 
subroutine dpttrf (N, D, E, INFO)
 DPTTRF More...
 
subroutine dpttrs (N, NRHS, D, E, B, LDB, INFO)
 DPTTRS More...
 
subroutine dptts2 (N, NRHS, D, E, B, LDB)
 DPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf. More...
 

Detailed Description

This is the group of double computational functions for PT matrices

Function Documentation

subroutine dptcon ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision  ANORM,
double precision  RCOND,
double precision, dimension( * )  WORK,
integer  INFO 
)

DPTCON

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

Purpose:
 DPTCON computes the reciprocal of the condition number (in the
 1-norm) of a real symmetric positive definite tridiagonal matrix
 using the factorization A = L*D*L**T or A = U**T*D*U computed by
 DPTTRF.

 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 DPTTRF.
[in]E
          E is DOUBLE PRECISION 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 DPTTRF.
[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]WORK
          WORK 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.
Date
September 2012
Further Details:
  The method used is described in Nicholas J. Higham, "Efficient
  Algorithms for Computing the Condition Number of a Tridiagonal
  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.

Definition at line 120 of file dptcon.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpteqr ( character  COMPZ,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer  INFO 
)

DPTEQR

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

Purpose:
 DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
 symmetric positive definite tridiagonal matrix by first factoring the
 matrix using DPTTRF, and then calling DBDSQR 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 symmetric positive definite matrix
 can also be found if DSYTRD, DSPTRD, or DSBTRD 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 symmetric
                  matrix also.  Array Z contains the orthogonal
                  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 DOUBLE PRECISION array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix used in the
          reduction to tridiagonal form.
          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
          original symmetric 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.
Date
September 2012

Definition at line 147 of file dpteqr.f.

147 *
148 * -- LAPACK computational routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  CHARACTER compz
155  INTEGER info, ldz, n
156 * ..
157 * .. Array Arguments ..
158  DOUBLE PRECISION d( * ), e( * ), work( * ), z( ldz, * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  DOUBLE PRECISION zero, one
165  parameter( zero = 0.0d0, one = 1.0d0 )
166 * ..
167 * .. External Functions ..
168  LOGICAL lsame
169  EXTERNAL lsame
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL dbdsqr, dlaset, dpttrf, xerbla
173 * ..
174 * .. Local Arrays ..
175  DOUBLE PRECISION c( 1, 1 ), vt( 1, 1 )
176 * ..
177 * .. Local Scalars ..
178  INTEGER i, icompz, nru
179 * ..
180 * .. Intrinsic Functions ..
181  INTRINSIC max, sqrt
182 * ..
183 * .. Executable Statements ..
184 *
185 * Test the input parameters.
186 *
187  info = 0
188 *
189  IF( lsame( compz, 'N' ) ) THEN
190  icompz = 0
191  ELSE IF( lsame( compz, 'V' ) ) THEN
192  icompz = 1
193  ELSE IF( lsame( compz, 'I' ) ) THEN
194  icompz = 2
195  ELSE
196  icompz = -1
197  END IF
198  IF( icompz.LT.0 ) THEN
199  info = -1
200  ELSE IF( n.LT.0 ) THEN
201  info = -2
202  ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
203  $ n ) ) ) THEN
204  info = -6
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'DPTEQR', -info )
208  RETURN
209  END IF
210 *
211 * Quick return if possible
212 *
213  IF( n.EQ.0 )
214  $ RETURN
215 *
216  IF( n.EQ.1 ) THEN
217  IF( icompz.GT.0 )
218  $ z( 1, 1 ) = one
219  RETURN
220  END IF
221  IF( icompz.EQ.2 )
222  $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
223 *
224 * Call DPTTRF to factor the matrix.
225 *
226  CALL dpttrf( n, d, e, info )
227  IF( info.NE.0 )
228  $ RETURN
229  DO 10 i = 1, n
230  d( i ) = sqrt( d( i ) )
231  10 CONTINUE
232  DO 20 i = 1, n - 1
233  e( i ) = e( i )*d( i )
234  20 CONTINUE
235 *
236 * Call DBDSQR to compute the singular values/vectors of the
237 * bidiagonal factor.
238 *
239  IF( icompz.GT.0 ) THEN
240  nru = n
241  ELSE
242  nru = 0
243  END IF
244  CALL dbdsqr( 'Lower', n, 0, nru, 0, d, e, vt, 1, z, ldz, c, 1,
245  $ work, info )
246 *
247 * Square the singular values.
248 *
249  IF( info.EQ.0 ) THEN
250  DO 30 i = 1, n
251  d( i ) = d( i )*d( i )
252  30 CONTINUE
253  ELSE
254  info = n + info
255  END IF
256 *
257  RETURN
258 *
259 * End of DPTEQR
260 *
subroutine dpttrf(N, D, E, INFO)
DPTTRF
Definition: dpttrf.f:93
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
Definition: dbdsqr.f:232

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dptrfs ( integer  N,
integer  NRHS,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  DF,
double precision, dimension( * )  EF,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
double precision, dimension( * )  WORK,
integer  INFO 
)

DPTRFS

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

Purpose:
 DPTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 and tridiagonal, and provides error bounds and backward error
 estimates for the solution.
Parameters
[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 diagonal elements of the tridiagonal matrix A.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix A.
[in]DF
          DF is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the diagonal matrix D from the
          factorization computed by DPTTRF.
[in]EF
          EF is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the unit bidiagonal factor
          L from the factorization computed by DPTTRF.
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DPTTRS.
          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 DOUBLE PRECISION array, dimension (2*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.
Date
September 2012

Definition at line 165 of file dptrfs.f.

165 *
166 * -- LAPACK computational routine (version 3.4.2) --
167 * -- LAPACK is a software package provided by Univ. of Tennessee, --
168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 * September 2012
170 *
171 * .. Scalar Arguments ..
172  INTEGER info, ldb, ldx, n, nrhs
173 * ..
174 * .. Array Arguments ..
175  DOUBLE PRECISION b( ldb, * ), berr( * ), d( * ), df( * ),
176  $ e( * ), ef( * ), ferr( * ), work( * ),
177  $ x( ldx, * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  INTEGER itmax
184  parameter( itmax = 5 )
185  DOUBLE PRECISION zero
186  parameter( zero = 0.0d+0 )
187  DOUBLE PRECISION one
188  parameter( one = 1.0d+0 )
189  DOUBLE PRECISION two
190  parameter( two = 2.0d+0 )
191  DOUBLE PRECISION three
192  parameter( three = 3.0d+0 )
193 * ..
194 * .. Local Scalars ..
195  INTEGER count, i, ix, j, nz
196  DOUBLE PRECISION bi, cx, dx, eps, ex, lstres, s, safe1, safe2,
197  $ safmin
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL daxpy, dpttrs, xerbla
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC abs, max
204 * ..
205 * .. External Functions ..
206  INTEGER idamax
207  DOUBLE PRECISION dlamch
208  EXTERNAL idamax, dlamch
209 * ..
210 * .. Executable Statements ..
211 *
212 * Test the input parameters.
213 *
214  info = 0
215  IF( n.LT.0 ) THEN
216  info = -1
217  ELSE IF( nrhs.LT.0 ) THEN
218  info = -2
219  ELSE IF( ldb.LT.max( 1, n ) ) THEN
220  info = -8
221  ELSE IF( ldx.LT.max( 1, n ) ) THEN
222  info = -10
223  END IF
224  IF( info.NE.0 ) THEN
225  CALL xerbla( 'DPTRFS', -info )
226  RETURN
227  END IF
228 *
229 * Quick return if possible
230 *
231  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
232  DO 10 j = 1, nrhs
233  ferr( j ) = zero
234  berr( j ) = zero
235  10 CONTINUE
236  RETURN
237  END IF
238 *
239 * NZ = maximum number of nonzero elements in each row of A, plus 1
240 *
241  nz = 4
242  eps = dlamch( 'Epsilon' )
243  safmin = dlamch( 'Safe minimum' )
244  safe1 = nz*safmin
245  safe2 = safe1 / eps
246 *
247 * Do for each right hand side
248 *
249  DO 90 j = 1, nrhs
250 *
251  count = 1
252  lstres = three
253  20 CONTINUE
254 *
255 * Loop until stopping criterion is satisfied.
256 *
257 * Compute residual R = B - A * X. Also compute
258 * abs(A)*abs(x) + abs(b) for use in the backward error bound.
259 *
260  IF( n.EQ.1 ) THEN
261  bi = b( 1, j )
262  dx = d( 1 )*x( 1, j )
263  work( n+1 ) = bi - dx
264  work( 1 ) = abs( bi ) + abs( dx )
265  ELSE
266  bi = b( 1, j )
267  dx = d( 1 )*x( 1, j )
268  ex = e( 1 )*x( 2, j )
269  work( n+1 ) = bi - dx - ex
270  work( 1 ) = abs( bi ) + abs( dx ) + abs( ex )
271  DO 30 i = 2, n - 1
272  bi = b( i, j )
273  cx = e( i-1 )*x( i-1, j )
274  dx = d( i )*x( i, j )
275  ex = e( i )*x( i+1, j )
276  work( n+i ) = bi - cx - dx - ex
277  work( i ) = abs( bi ) + abs( cx ) + abs( dx ) + abs( ex )
278  30 CONTINUE
279  bi = b( n, j )
280  cx = e( n-1 )*x( n-1, j )
281  dx = d( n )*x( n, j )
282  work( n+n ) = bi - cx - dx
283  work( n ) = abs( bi ) + abs( cx ) + abs( dx )
284  END IF
285 *
286 * Compute componentwise relative backward error from formula
287 *
288 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
289 *
290 * where abs(Z) is the componentwise absolute value of the matrix
291 * or vector Z. If the i-th component of the denominator is less
292 * than SAFE2, then SAFE1 is added to the i-th components of the
293 * numerator and denominator before dividing.
294 *
295  s = zero
296  DO 40 i = 1, n
297  IF( work( i ).GT.safe2 ) THEN
298  s = max( s, abs( work( n+i ) ) / work( i ) )
299  ELSE
300  s = max( s, ( abs( work( n+i ) )+safe1 ) /
301  $ ( work( i )+safe1 ) )
302  END IF
303  40 CONTINUE
304  berr( j ) = s
305 *
306 * Test stopping criterion. Continue iterating if
307 * 1) The residual BERR(J) is larger than machine epsilon, and
308 * 2) BERR(J) decreased by at least a factor of 2 during the
309 * last iteration, and
310 * 3) At most ITMAX iterations tried.
311 *
312  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
313  $ count.LE.itmax ) THEN
314 *
315 * Update solution and try again.
316 *
317  CALL dpttrs( n, 1, df, ef, work( n+1 ), n, info )
318  CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
319  lstres = berr( j )
320  count = count + 1
321  GO TO 20
322  END IF
323 *
324 * Bound error from formula
325 *
326 * norm(X - XTRUE) / norm(X) .le. FERR =
327 * norm( abs(inv(A))*
328 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
329 *
330 * where
331 * norm(Z) is the magnitude of the largest component of Z
332 * inv(A) is the inverse of A
333 * abs(Z) is the componentwise absolute value of the matrix or
334 * vector Z
335 * NZ is the maximum number of nonzeros in any row of A, plus 1
336 * EPS is machine epsilon
337 *
338 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
339 * is incremented by SAFE1 if the i-th component of
340 * abs(A)*abs(X) + abs(B) is less than SAFE2.
341 *
342  DO 50 i = 1, n
343  IF( work( i ).GT.safe2 ) THEN
344  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
345  ELSE
346  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
347  END IF
348  50 CONTINUE
349  ix = idamax( n, work, 1 )
350  ferr( j ) = work( ix )
351 *
352 * Estimate the norm of inv(A).
353 *
354 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
355 *
356 * m(i,j) = abs(A(i,j)), i = j,
357 * m(i,j) = -abs(A(i,j)), i .ne. j,
358 *
359 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
360 *
361 * Solve M(L) * x = e.
362 *
363  work( 1 ) = one
364  DO 60 i = 2, n
365  work( i ) = one + work( i-1 )*abs( ef( i-1 ) )
366  60 CONTINUE
367 *
368 * Solve D * M(L)**T * x = b.
369 *
370  work( n ) = work( n ) / df( n )
371  DO 70 i = n - 1, 1, -1
372  work( i ) = work( i ) / df( i ) + work( i+1 )*abs( ef( i ) )
373  70 CONTINUE
374 *
375 * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
376 *
377  ix = idamax( n, work, 1 )
378  ferr( j ) = ferr( j )*abs( work( ix ) )
379 *
380 * Normalize error.
381 *
382  lstres = zero
383  DO 80 i = 1, n
384  lstres = max( lstres, abs( x( i, j ) ) )
385  80 CONTINUE
386  IF( lstres.NE.zero )
387  $ ferr( j ) = ferr( j ) / lstres
388 *
389  90 CONTINUE
390 *
391  RETURN
392 *
393 * End of DPTRFS
394 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dpttrs(N, NRHS, D, E, B, LDB, INFO)
DPTTRS
Definition: dpttrs.f:111
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpttrf ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
integer  INFO 
)

DPTTRF

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

Purpose:
 DPTTRF computes the L*D*L**T factorization of a real symmetric
 positive definite tridiagonal matrix A.  The factorization may also
 be regarded as having the form A = U**T*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**T factorization of A.
[in,out]E
          E is DOUBLE PRECISION 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**T factorization of A.
          E can also be regarded as the superdiagonal of the unit
          bidiagonal factor U from the U**T*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.
Date
September 2012

Definition at line 93 of file dpttrf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dpttrs ( integer  N,
integer  NRHS,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

DPTTRS

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

Purpose:
 DPTTRS solves a tridiagonal system of the form
    A * X = B
 using the L*D*L**T factorization of A computed by DPTTRF.  D is a
 diagonal matrix specified in the vector D, L is a unit bidiagonal
 matrix whose subdiagonal is specified in the vector E, and X and B
 are N by NRHS matrices.
Parameters
[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
          L*D*L**T factorization of A.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the unit bidiagonal factor
          L from the L*D*L**T factorization of A.  E can also be regarded
          as the superdiagonal of the unit bidiagonal factor U from the
          factorization A = U**T*D*U.
[in,out]B
          B is DOUBLE PRECISION 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.
Date
September 2012

Definition at line 111 of file dpttrs.f.

111 *
112 * -- LAPACK computational routine (version 3.4.2) --
113 * -- LAPACK is a software package provided by Univ. of Tennessee, --
114 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115 * September 2012
116 *
117 * .. Scalar Arguments ..
118  INTEGER info, ldb, n, nrhs
119 * ..
120 * .. Array Arguments ..
121  DOUBLE PRECISION b( ldb, * ), d( * ), e( * )
122 * ..
123 *
124 * =====================================================================
125 *
126 * .. Local Scalars ..
127  INTEGER j, jb, nb
128 * ..
129 * .. External Functions ..
130  INTEGER ilaenv
131  EXTERNAL ilaenv
132 * ..
133 * .. External Subroutines ..
134  EXTERNAL dptts2, xerbla
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC max, min
138 * ..
139 * .. Executable Statements ..
140 *
141 * Test the input arguments.
142 *
143  info = 0
144  IF( n.LT.0 ) THEN
145  info = -1
146  ELSE IF( nrhs.LT.0 ) THEN
147  info = -2
148  ELSE IF( ldb.LT.max( 1, n ) ) THEN
149  info = -6
150  END IF
151  IF( info.NE.0 ) THEN
152  CALL xerbla( 'DPTTRS', -info )
153  RETURN
154  END IF
155 *
156 * Quick return if possible
157 *
158  IF( n.EQ.0 .OR. nrhs.EQ.0 )
159  $ RETURN
160 *
161 * Determine the number of right-hand sides to solve at a time.
162 *
163  IF( nrhs.EQ.1 ) THEN
164  nb = 1
165  ELSE
166  nb = max( 1, ilaenv( 1, 'DPTTRS', ' ', n, nrhs, -1, -1 ) )
167  END IF
168 *
169  IF( nb.GE.nrhs ) THEN
170  CALL dptts2( n, nrhs, d, e, b, ldb )
171  ELSE
172  DO 10 j = 1, nrhs, nb
173  jb = min( nrhs-j+1, nb )
174  CALL dptts2( n, jb, d, e, b( 1, j ), ldb )
175  10 CONTINUE
176  END IF
177 *
178  RETURN
179 *
180 * End of DPTTRS
181 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dptts2(N, NRHS, D, E, B, LDB)
DPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf...
Definition: dptts2.f:104
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dptts2 ( integer  N,
integer  NRHS,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldb, * )  B,
integer  LDB 
)

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

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

Purpose:
 DPTTS2 solves a tridiagonal system of the form
    A * X = B
 using the L*D*L**T factorization of A computed by DPTTRF.  D is a
 diagonal matrix specified in the vector D, L is a unit bidiagonal
 matrix whose subdiagonal is specified in the vector E, and X and B
 are N by NRHS matrices.
Parameters
[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
          L*D*L**T factorization of A.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the unit bidiagonal factor
          L from the L*D*L**T factorization of A.  E can also be regarded
          as the superdiagonal of the unit bidiagonal factor U from the
          factorization A = U**T*D*U.
[in,out]B
          B is DOUBLE PRECISION 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.
Date
September 2012

Definition at line 104 of file dptts2.f.

104 *
105 * -- LAPACK computational routine (version 3.4.2) --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 * September 2012
109 *
110 * .. Scalar Arguments ..
111  INTEGER ldb, n, nrhs
112 * ..
113 * .. Array Arguments ..
114  DOUBLE PRECISION b( ldb, * ), d( * ), e( * )
115 * ..
116 *
117 * =====================================================================
118 *
119 * .. Local Scalars ..
120  INTEGER i, j
121 * ..
122 * .. External Subroutines ..
123  EXTERNAL dscal
124 * ..
125 * .. Executable Statements ..
126 *
127 * Quick return if possible
128 *
129  IF( n.LE.1 ) THEN
130  IF( n.EQ.1 )
131  $ CALL dscal( nrhs, 1.d0 / d( 1 ), b, ldb )
132  RETURN
133  END IF
134 *
135 * Solve A * X = B using the factorization A = L*D*L**T,
136 * overwriting each right hand side vector with its solution.
137 *
138  DO 30 j = 1, nrhs
139 *
140 * Solve L * x = b.
141 *
142  DO 10 i = 2, n
143  b( i, j ) = b( i, j ) - b( i-1, j )*e( i-1 )
144  10 CONTINUE
145 *
146 * Solve D * L**T * x = b.
147 *
148  b( n, j ) = b( n, j ) / d( n )
149  DO 20 i = n - 1, 1, -1
150  b( i, j ) = b( i, j ) / d( i ) - b( i+1, j )*e( i )
151  20 CONTINUE
152  30 CONTINUE
153 *
154  RETURN
155 *
156 * End of DPTTS2
157 *
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: