PLASMA  2.4.5
PLASMA - Parallel Linear Algebra for Scalable Multi-core Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
dget02.f
Go to the documentation of this file.
1  SUBROUTINE dget02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB,
2  $ rwork, 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 trans
10  INTEGER lda, ldb, ldx, m, 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 * DGET02 computes the residual for a solution of a system of linear
22 * equations A*x = b or A'*x = b:
23 * RESID = norm( B - A*X ) / ( norm(A) * norm(X) + norm(RHS))* N * EPS )
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANS (input) CHARACTER*1
30 * Specifies the form of the system of equations:
31 * = 'N': A *x = b
32 * = 'T': A'*x = b, where A' is the transpose of A
33 * = 'C': A'*x = b, where A' is the transpose of A
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of columns of the matrix A. N >= 0.
40 *
41 * NRHS (input) INTEGER
42 * The number of columns of B, the matrix of right hand sides.
43 * NRHS >= 0.
44 *
45 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
46 * The original M x N matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
52 * The computed solution vectors for the system of linear
53 * equations.
54 *
55 * LDX (input) INTEGER
56 * The leading dimension of the array X. If TRANS = 'N',
57 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
58 *
59 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
60 * On entry, the right hand side vectors for the system of
61 * linear equations.
62 * On exit, B is overwritten with the difference B - A*X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. IF TRANS = 'N',
66 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * The maximum over the number of right hand sides of
72 * norm( B - A*X ) / ( norm(A) * norm(X) + norm(RHS))* N * EPS )
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77  DOUBLE PRECISION zero, one
78  parameter( zero = 0.0d+0, one = 1.0d+0 )
79 * ..
80 * .. Local Scalars ..
81  INTEGER j, n1, n2
82  DOUBLE PRECISION anorm, bnorm, rhsnorm, eps, xnorm
83 * ..
84 * .. External Functions ..
85  LOGICAL lsame
86  DOUBLE PRECISION dasum, dlamch, dlange
87  EXTERNAL lsame, dasum, dlamch, dlange
88 * ..
89 * .. External Subroutines ..
90  EXTERNAL dgemm
91 * ..
92 * .. Intrinsic Functions ..
93  INTRINSIC max
94 * ..
95 * .. Executable Statements ..
96 *
97 * Quick exit if M = 0 or N = 0 or NRHS = 0
98 *
99  IF( m.LE.0 .OR. n.LE.0 .OR. nrhs.EQ.0 ) THEN
100  resid = zero
101  return
102  END IF
103 *
104  IF( lsame( trans, 'T' ) .OR. lsame( trans, 'C' ) ) THEN
105  n1 = n
106  n2 = m
107  ELSE
108  n1 = m
109  n2 = n
110  END IF
111 *
112 * Exit with RESID = 1/EPS if ANORM = 0.
113 *
114  eps = dlamch( 'Epsilon' )
115  anorm = dlange( '1', n1, n2, a, lda, rwork )
116  rhsnorm = dlange( '1', n1, nrhs, b, ldb, rwork )
117  IF( anorm.LE.zero ) THEN
118  resid = one / eps
119  return
120  END IF
121 *
122 * Compute B - A*X (or B - A'*X ) and store in B.
123 *
124  CALL dgemm( trans, 'No transpose', n1, nrhs, n2, -one, a, lda, x,
125  $ ldx, one, b, ldb )
126 *
127 * Compute the maximum over the number of right hand sides of
128 * norm( B - A*X ) / ( norm(A) * norm(X) + norm(RHS))* N * EPS ) .
129 *
130  resid = zero
131  DO 10 j = 1, nrhs
132  bnorm = dasum( n1, b( 1, j ), 1 )
133  xnorm = dasum( n2, x( 1, j ), 1 )
134  IF( xnorm.LE.zero ) THEN
135  resid = one / eps
136  ELSE
137  resid = max( resid, ( bnorm) / ((anorm * xnorm + rhsnorm)*
138  $ n1 *eps ))
139  END IF
140  10 continue
141 *
142  return
143 *
144 * End of DGET02
145 *
146  END