MAGMA  1.2.0 MatrixAlgebraonGPUandMulticoreArchitectures
dpot02.f
Go to the documentation of this file.
1  SUBROUTINE dpot02( UPLO, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK,
2  \$ resid )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9  CHARACTER uplo
10  INTEGER lda, ldb, ldx, n, nrhs
11  DOUBLE PRECISION resid
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION a( lda, * ), b( ldb, * ), rwork( * ),
15  \$ x( ldx, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DPOT02 computes the residual for the solution of a symmetric system
22 * of linear equations A*x = b:
23 *
24 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
25 *
26 * where EPS is the machine epsilon.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the upper or lower triangular part of the
33 * symmetric matrix A is stored:
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The number of rows and columns of the matrix A. N >= 0.
39 *
40 * NRHS (input) INTEGER
41 * The number of columns of B, the matrix of right hand sides.
42 * NRHS >= 0.
43 *
44 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
45 * The original symmetric matrix A.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,N)
49 *
50 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
51 * The computed solution vectors for the system of linear
52 * equations.
53 *
54 * LDX (input) INTEGER
55 * The leading dimension of the array X. LDX >= max(1,N).
56 *
57 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
58 * On entry, the right hand side vectors for the system of
59 * linear equations.
60 * On exit, B is overwritten with the difference B - A*X.
61 *
62 * LDB (input) INTEGER
63 * The leading dimension of the array B. LDB >= max(1,N).
64 *
65 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
66 *
67 * RESID (output) DOUBLE PRECISION
68 * The maximum over the number of right hand sides of
69 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
70 *
71 * =====================================================================
72 *
73 * .. Parameters ..
74  DOUBLE PRECISION zero, one
75  parameter( zero = 0.0d+0, one = 1.0d+0 )
76 * ..
77 * .. Local Scalars ..
78  INTEGER j
79  DOUBLE PRECISION anorm, bnorm, eps, xnorm
80 * ..
81 * .. External Functions ..
82  DOUBLE PRECISION dasum, dlamch, dlansy
83  EXTERNAL dasum, dlamch, dlansy
84 * ..
85 * .. External Subroutines ..
86  EXTERNAL dsymm
87 * ..
88 * .. Intrinsic Functions ..
89  INTRINSIC max
90 * ..
91 * .. Executable Statements ..
92 *
93 * Quick exit if N = 0 or NRHS = 0.
94 *
95  IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
96  resid = zero
97  return
98  END IF
99 *
100 * Exit with RESID = 1/EPS if ANORM = 0.
101 *
102  eps = dlamch( 'Epsilon' )
103  anorm = dlansy( '1', uplo, n, a, lda, rwork )
104  IF( anorm.LE.zero ) THEN
105  resid = one / eps
106  return
107  END IF
108 *
109 * Compute B - A*X
110 *
111  CALL dsymm( 'Left', uplo, n, nrhs, -one, a, lda, x, ldx, one, b,
112  \$ ldb )
113 *
114 * Compute the maximum over the number of right hand sides of
115 * norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
116 *
117  resid = zero
118  DO 10 j = 1, nrhs
119  bnorm = dasum( n, b( 1, j ), 1 )
120  xnorm = dasum( n, x( 1, j ), 1 )
121  IF( xnorm.LE.zero ) THEN
122  resid = one / eps
123  ELSE
124  resid = max( resid, ( ( bnorm / anorm ) / xnorm ) / eps )
125  END IF
126  10 continue
127 *
128  return
129 *
130 * End of DPOT02
131 *
132  END