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

Functions

subroutine sptcon (N, D, E, ANORM, RCOND, WORK, INFO)
 SPTCON More...
 
subroutine spteqr (COMPZ, N, D, E, Z, LDZ, WORK, INFO)
 SPTEQR More...
 
subroutine sptrfs (N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, INFO)
 SPTRFS More...
 
subroutine spttrs (N, NRHS, D, E, B, LDB, INFO)
 SPTTRS More...
 
subroutine sptts2 (N, NRHS, D, E, B, LDB)
 SPTTS2 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 real computational functions for PT matrices

Function Documentation

subroutine sptcon ( integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real  ANORM,
real  RCOND,
real, dimension( * )  WORK,
integer  INFO 
)

SPTCON

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

Purpose:
 SPTCON 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
 SPTTRF.

 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 SPTTRF.
[in]E
          E is REAL 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 SPTTRF.
[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]WORK
          WORK 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.
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 sptcon.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  REAL anorm, rcond
129 * ..
130 * .. Array Arguments ..
131  REAL d( * ), e( * ), work( * )
132 * ..
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  REAL one, zero
138  parameter( one = 1.0e+0, zero = 0.0e+0 )
139 * ..
140 * .. Local Scalars ..
141  INTEGER i, ix
142  REAL ainvnm
143 * ..
144 * .. External Functions ..
145  INTEGER isamax
146  EXTERNAL isamax
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( 'SPTCON', -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 = isamax( 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 SPTCON
220 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
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 spteqr ( character  COMPZ,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  INFO 
)

SPTEQR

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

Purpose:
 SPTEQR computes all eigenvalues and, optionally, eigenvectors of a
 symmetric positive definite tridiagonal matrix by first factoring the
 matrix using SPTTRF, and then calling SBDSQR 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 SSYTRD, SSPTRD, or SSBTRD 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 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 REAL 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 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.
Date
September 2012

Definition at line 147 of file spteqr.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  REAL d( * ), e( * ), work( * ), z( ldz, * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  REAL zero, one
165  parameter( zero = 0.0e0, one = 1.0e0 )
166 * ..
167 * .. External Functions ..
168  LOGICAL lsame
169  EXTERNAL lsame
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL sbdsqr, slaset, spttrf, xerbla
173 * ..
174 * .. Local Arrays ..
175  REAL 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( 'SPTEQR', -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 slaset( 'Full', n, n, zero, one, z, ldz )
223 *
224 * Call SPTTRF to factor the matrix.
225 *
226  CALL spttrf( 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 SBDSQR 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 sbdsqr( '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 SPTEQR
260 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:232
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine spttrf(N, D, E, INFO)
SPTTRF
Definition: spttrf.f:93
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112

Here is the call graph for this function:

Here is the caller graph for this function:

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

SPTRFS

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

Purpose:
 SPTRFS 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 REAL array, dimension (N)
          The n diagonal elements of the tridiagonal matrix A.
[in]E
          E is REAL array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix A.
[in]DF
          DF is REAL array, dimension (N)
          The n diagonal elements of the diagonal matrix D from the
          factorization computed by SPTTRF.
[in]EF
          EF is REAL array, dimension (N-1)
          The (n-1) subdiagonal elements of the unit bidiagonal factor
          L from the factorization computed by SPTTRF.
[in]B
          B is REAL 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 REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SPTTRS.
          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 REAL 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 sptrfs.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  REAL 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  REAL zero
186  parameter( zero = 0.0e+0 )
187  REAL one
188  parameter( one = 1.0e+0 )
189  REAL two
190  parameter( two = 2.0e+0 )
191  REAL three
192  parameter( three = 3.0e+0 )
193 * ..
194 * .. Local Scalars ..
195  INTEGER count, i, ix, j, nz
196  REAL bi, cx, dx, eps, ex, lstres, s, safe1, safe2,
197  $ safmin
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL saxpy, spttrs, xerbla
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC abs, max
204 * ..
205 * .. External Functions ..
206  INTEGER isamax
207  REAL slamch
208  EXTERNAL isamax, slamch
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( 'SPTRFS', -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 = slamch( 'Epsilon' )
243  safmin = slamch( '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 spttrs( n, 1, df, ef, work( n+1 ), n, info )
318  CALL saxpy( 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 = isamax( 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 = isamax( 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 SPTRFS
394 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
subroutine spttrs(N, NRHS, D, E, B, LDB, INFO)
SPTTRS
Definition: spttrs.f:111

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine spttrs ( integer  N,
integer  NRHS,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

SPTTRS

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

Purpose:
 SPTTRS solves a tridiagonal system of the form
    A * X = B
 using the L*D*L**T factorization of A computed by SPTTRF.  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 REAL 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 REAL 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 REAL 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 spttrs.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  REAL 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 sptts2, 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( 'SPTTRS', -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, 'SPTTRS', ' ', n, nrhs, -1, -1 ) )
167  END IF
168 *
169  IF( nb.GE.nrhs ) THEN
170  CALL sptts2( n, nrhs, d, e, b, ldb )
171  ELSE
172  DO 10 j = 1, nrhs, nb
173  jb = min( nrhs-j+1, nb )
174  CALL sptts2( n, jb, d, e, b( 1, j ), ldb )
175  10 CONTINUE
176  END IF
177 *
178  RETURN
179 *
180 * End of SPTTRS
181 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine sptts2(N, NRHS, D, E, B, LDB)
SPTTS2 solves a tridiagonal system of the form AX=B using the L D LH factorization computed by spttrf...
Definition: sptts2.f:104

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sptts2 ( integer  N,
integer  NRHS,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldb, * )  B,
integer  LDB 
)

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

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

Purpose:
 SPTTS2 solves a tridiagonal system of the form
    A * X = B
 using the L*D*L**T factorization of A computed by SPTTRF.  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 REAL 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 REAL 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 REAL 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 sptts2.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  REAL b( ldb, * ), d( * ), e( * )
115 * ..
116 *
117 * =====================================================================
118 *
119 * .. Local Scalars ..
120  INTEGER i, j
121 * ..
122 * .. External Subroutines ..
123  EXTERNAL sscal
124 * ..
125 * .. Executable Statements ..
126 *
127 * Quick return if possible
128 *
129  IF( n.LE.1 ) THEN
130  IF( n.EQ.1 )
131  $ CALL sscal( nrhs, 1. / 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 SPTTS2
157 *
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: