LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dsposv()

subroutine dsposv ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( n, * )  WORK,
real, dimension( * )  SWORK,
integer  ITER,
integer  INFO 
)

DSPOSV computes the solution to system of linear equations A * X = B for PO matrices

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

Purpose:
 DSPOSV computes the solution to a real system of linear equations
    A * X = B,
 where A is an N-by-N symmetric positive definite matrix and X and B
 are N-by-NRHS matrices.

 DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
 and use this factorization within an iterative refinement procedure
 to produce a solution with DOUBLE PRECISION normwise backward error
 quality (see below). If the approach fails the method switches to a
 DOUBLE PRECISION factorization and solve.

 The iterative refinement is not going to be a winning strategy if
 the ratio SINGLE PRECISION performance over DOUBLE PRECISION
 performance is too small. A reasonable strategy should take the
 number of right-hand sides and the size of the matrix into account.
 This might be done with a call to ILAENV in the future. Up to now, we
 always try iterative refinement.

 The iterative refinement process is stopped if
     ITER > ITERMAX
 or for all the RHS we have:
     RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
 where
     o ITER is the number of the current iteration in the iterative
       refinement process
     o RNRM is the infinity-norm of the residual
     o XNRM is the infinity-norm of the solution
     o ANRM is the infinity-operator-norm of the matrix A
     o EPS is the machine epsilon returned by DLAMCH('Epsilon')
 The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
 respectively.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The number of linear equations, i.e., 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,out]A
          A is DOUBLE PRECISION array,
          dimension (LDA,N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
          N-by-N upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading N-by-N lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.
          On exit, if iterative refinement has been successfully used
          (INFO = 0 and ITER >= 0, see description below), then A is
          unchanged, if double precision factorization has been used
          (INFO = 0 and ITER < 0, see description below), then the
          array A contains the factor U or L from the Cholesky
          factorization A = U**T*U or A = L*L**T.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          The N-by-NRHS right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]X
          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
          If INFO = 0, the N-by-NRHS solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N,NRHS)
          This array is used to hold the residual vectors.
[out]SWORK
          SWORK is REAL array, dimension (N*(N+NRHS))
          This array is used to use the single precision matrix and the
          right-hand sides or solutions in single precision.
[out]ITER
          ITER is INTEGER
          < 0: iterative refinement has failed, double precision
               factorization has been performed
               -1 : the routine fell back to full precision for
                    implementation- or machine-specific reasons
               -2 : narrowing the precision induced an overflow,
                    the routine fell back to full precision
               -3 : failure of SPOTRF
               -31: stop the iterative refinement after the 30th
                    iterations
          > 0: iterative refinement has been successfully used.
               Returns the number of iterations
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the leading minor of order i of (DOUBLE
                PRECISION) A is not positive definite, so the
                factorization could not be completed, and the solution
                has not been computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 197 of file dsposv.f.

199 *
200 * -- LAPACK driver routine --
201 * -- LAPACK is a software package provided by Univ. of Tennessee, --
202 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203 *
204 * .. Scalar Arguments ..
205  CHARACTER UPLO
206  INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
207 * ..
208 * .. Array Arguments ..
209  REAL SWORK( * )
210  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
211  $ X( LDX, * )
212 * ..
213 *
214 * =====================================================================
215 *
216 * .. Parameters ..
217  LOGICAL DOITREF
218  parameter( doitref = .true. )
219 *
220  INTEGER ITERMAX
221  parameter( itermax = 30 )
222 *
223  DOUBLE PRECISION BWDMAX
224  parameter( bwdmax = 1.0e+00 )
225 *
226  DOUBLE PRECISION NEGONE, ONE
227  parameter( negone = -1.0d+0, one = 1.0d+0 )
228 *
229 * .. Local Scalars ..
230  INTEGER I, IITER, PTSA, PTSX
231  DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
232 *
233 * .. External Subroutines ..
234  EXTERNAL daxpy, dsymm, dlacpy, dlat2s, dlag2s, slag2d,
236 * ..
237 * .. External Functions ..
238  INTEGER IDAMAX
239  DOUBLE PRECISION DLAMCH, DLANSY
240  LOGICAL LSAME
241  EXTERNAL idamax, dlamch, dlansy, lsame
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC abs, dble, max, sqrt
245 * ..
246 * .. Executable Statements ..
247 *
248  info = 0
249  iter = 0
250 *
251 * Test the input parameters.
252 *
253  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
254  info = -1
255  ELSE IF( n.LT.0 ) THEN
256  info = -2
257  ELSE IF( nrhs.LT.0 ) THEN
258  info = -3
259  ELSE IF( lda.LT.max( 1, n ) ) THEN
260  info = -5
261  ELSE IF( ldb.LT.max( 1, n ) ) THEN
262  info = -7
263  ELSE IF( ldx.LT.max( 1, n ) ) THEN
264  info = -9
265  END IF
266  IF( info.NE.0 ) THEN
267  CALL xerbla( 'DSPOSV', -info )
268  RETURN
269  END IF
270 *
271 * Quick return if (N.EQ.0).
272 *
273  IF( n.EQ.0 )
274  $ RETURN
275 *
276 * Skip single precision iterative refinement if a priori slower
277 * than double precision factorization.
278 *
279  IF( .NOT.doitref ) THEN
280  iter = -1
281  GO TO 40
282  END IF
283 *
284 * Compute some constants.
285 *
286  anrm = dlansy( 'I', uplo, n, a, lda, work )
287  eps = dlamch( 'Epsilon' )
288  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
289 *
290 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
291 *
292  ptsa = 1
293  ptsx = ptsa + n*n
294 *
295 * Convert B from double precision to single precision and store the
296 * result in SX.
297 *
298  CALL dlag2s( n, nrhs, b, ldb, swork( ptsx ), n, info )
299 *
300  IF( info.NE.0 ) THEN
301  iter = -2
302  GO TO 40
303  END IF
304 *
305 * Convert A from double precision to single precision and store the
306 * result in SA.
307 *
308  CALL dlat2s( uplo, n, a, lda, swork( ptsa ), n, info )
309 *
310  IF( info.NE.0 ) THEN
311  iter = -2
312  GO TO 40
313  END IF
314 *
315 * Compute the Cholesky factorization of SA.
316 *
317  CALL spotrf( uplo, n, swork( ptsa ), n, info )
318 *
319  IF( info.NE.0 ) THEN
320  iter = -3
321  GO TO 40
322  END IF
323 *
324 * Solve the system SA*SX = SB.
325 *
326  CALL spotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
327  $ info )
328 *
329 * Convert SX back to double precision
330 *
331  CALL slag2d( n, nrhs, swork( ptsx ), n, x, ldx, info )
332 *
333 * Compute R = B - AX (R is WORK).
334 *
335  CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
336 *
337  CALL dsymm( 'Left', uplo, n, nrhs, negone, a, lda, x, ldx, one,
338  $ work, n )
339 *
340 * Check whether the NRHS normwise backward errors satisfy the
341 * stopping criterion. If yes, set ITER=0 and return.
342 *
343  DO i = 1, nrhs
344  xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
345  rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
346  IF( rnrm.GT.xnrm*cte )
347  $ GO TO 10
348  END DO
349 *
350 * If we are here, the NRHS normwise backward errors satisfy the
351 * stopping criterion. We are good to exit.
352 *
353  iter = 0
354  RETURN
355 *
356  10 CONTINUE
357 *
358  DO 30 iiter = 1, itermax
359 *
360 * Convert R (in WORK) from double precision to single precision
361 * and store the result in SX.
362 *
363  CALL dlag2s( n, nrhs, work, n, swork( ptsx ), n, info )
364 *
365  IF( info.NE.0 ) THEN
366  iter = -2
367  GO TO 40
368  END IF
369 *
370 * Solve the system SA*SX = SR.
371 *
372  CALL spotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
373  $ info )
374 *
375 * Convert SX back to double precision and update the current
376 * iterate.
377 *
378  CALL slag2d( n, nrhs, swork( ptsx ), n, work, n, info )
379 *
380  DO i = 1, nrhs
381  CALL daxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
382  END DO
383 *
384 * Compute R = B - AX (R is WORK).
385 *
386  CALL dlacpy( 'All', n, nrhs, b, ldb, work, n )
387 *
388  CALL dsymm( 'L', uplo, n, nrhs, negone, a, lda, x, ldx, one,
389  $ work, n )
390 *
391 * Check whether the NRHS normwise backward errors satisfy the
392 * stopping criterion. If yes, set ITER=IITER>0 and return.
393 *
394  DO i = 1, nrhs
395  xnrm = abs( x( idamax( n, x( 1, i ), 1 ), i ) )
396  rnrm = abs( work( idamax( n, work( 1, i ), 1 ), i ) )
397  IF( rnrm.GT.xnrm*cte )
398  $ GO TO 20
399  END DO
400 *
401 * If we are here, the NRHS normwise backward errors satisfy the
402 * stopping criterion, we are good to exit.
403 *
404  iter = iiter
405 *
406  RETURN
407 *
408  20 CONTINUE
409 *
410  30 CONTINUE
411 *
412 * If we are at this place of the code, this is because we have
413 * performed ITER=ITERMAX iterations and never satisfied the
414 * stopping criterion, set up the ITER flag accordingly and follow
415 * up on double precision routine.
416 *
417  iter = -itermax - 1
418 *
419  40 CONTINUE
420 *
421 * Single-precision iterative refinement failed to converge to a
422 * satisfactory solution, so we resort to double precision.
423 *
424  CALL dpotrf( uplo, n, a, lda, info )
425 *
426  IF( info.NE.0 )
427  $ RETURN
428 *
429  CALL dlacpy( 'All', n, nrhs, b, ldb, x, ldx )
430  CALL dpotrs( uplo, n, nrhs, a, lda, x, ldx, info )
431 *
432  RETURN
433 *
434 * End of DSPOSV
435 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine slag2d(M, N, SA, LDSA, A, LDA, INFO)
SLAG2D converts a single precision matrix to a double precision matrix.
Definition: slag2d.f:104
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dsymm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DSYMM
Definition: dsymm.f:189
subroutine dlat2s(UPLO, N, A, LDA, SA, LDSA, INFO)
DLAT2S converts a double-precision triangular matrix to a single-precision triangular matrix.
Definition: dlat2s.f:111
subroutine dlag2s(M, N, A, LDA, SA, LDSA, INFO)
DLAG2S converts a double precision matrix to a single precision matrix.
Definition: dlag2s.f:108
subroutine dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:110
subroutine dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:107
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlansy.f:122
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:107
subroutine spotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
SPOTRS
Definition: spotrs.f:110
Here is the call graph for this function:
Here is the caller graph for this function: