1 SUBROUTINE DGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, 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 INTEGER LDA, LDAFAC, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET01 reconstructs a matrix A from its L*U factorization and
21 * computes the residual
22 * norm(L*U - A) / ( N * norm(A) * EPS ),
23 * where EPS is the machine epsilon.
24 *
25 * Arguments
26 * ==========
27 *
28 * M (input) INTEGER
29 * The number of rows of the matrix A. M >= 0.
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A. N >= 0.
33 *
34 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
35 * The original M x N matrix A.
36 *
37 * LDA (input) INTEGER
38 * The leading dimension of the array A. LDA >= max(1,M).
39 *
40 * AFAC (input/output) DOUBLE PRECISION array, dimension (LDAFAC,N)
41 * The factored form of the matrix A. AFAC contains the factors
42 * L and U from the L*U factorization as computed by DGETRF.
43 * Overwritten with the reconstructed matrix, and then with the
44 * difference L*U - A.
45 *
46 * LDAFAC (input) INTEGER
47 * The leading dimension of the array AFAC. LDAFAC >= max(1,M).
48 *
49 * IPIV (input) INTEGER array, dimension (N)
50 * The pivot indices from DGETRF.
51 *
52 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
53 *
54 * RESID (output) DOUBLE PRECISION
55 * norm(L*U - A) / ( N * norm(A) * EPS )
56 *
57 * =====================================================================
58 *
59 *
60 * .. Parameters ..
61 DOUBLE PRECISION ZERO, ONE
62 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER I, J, K
66 DOUBLE PRECISION ANORM, EPS, T
67 * ..
68 * .. External Functions ..
69 DOUBLE PRECISION DDOT, DLAMCH, DLANGE
70 EXTERNAL DDOT, DLAMCH, DLANGE
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL DGEMV, DLASWP, DSCAL, DTRMV
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE, MIN
77 * ..
78 * .. Executable Statements ..
79 *
80 * Quick exit if M = 0 or N = 0.
81 *
82 IF( M.LE.0 .OR. N.LE.0 ) THEN
83 RESID = ZERO
84 RETURN
85 END IF
86 *
87 * Determine EPS and the norm of A.
88 *
89 EPS = DLAMCH( 'Epsilon' )
90 ANORM = DLANGE( '1', M, N, A, LDA, RWORK )
91 *
92 * Compute the product L*U and overwrite AFAC with the result.
93 * A column at a time of the product is obtained, starting with
94 * column N.
95 *
96 DO 10 K = N, 1, -1
97 IF( K.GT.M ) THEN
98 CALL DTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC,
99 $ LDAFAC, AFAC( 1, K ), 1 )
100 ELSE
101 *
102 * Compute elements (K+1:M,K)
103 *
104 T = AFAC( K, K )
105 IF( K+1.LE.M ) THEN
106 CALL DSCAL( M-K, T, AFAC( K+1, K ), 1 )
107 CALL DGEMV( 'No transpose', M-K, K-1, ONE,
108 $ AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1, ONE,
109 $ AFAC( K+1, K ), 1 )
110 END IF
111 *
112 * Compute the (K,K) element
113 *
114 AFAC( K, K ) = T + DDOT( K-1, AFAC( K, 1 ), LDAFAC,
115 $ AFAC( 1, K ), 1 )
116 *
117 * Compute elements (1:K-1,K)
118 *
119 CALL DTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC,
120 $ LDAFAC, AFAC( 1, K ), 1 )
121 END IF
122 10 CONTINUE
123 CALL DLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 )
124 *
125 * Compute the difference L*U - A and store in AFAC.
126 *
127 DO 30 J = 1, N
128 DO 20 I = 1, M
129 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
130 20 CONTINUE
131 30 CONTINUE
132 *
133 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
134 *
135 RESID = DLANGE( '1', M, N, AFAC, LDAFAC, RWORK )
136 *
137 IF( ANORM.LE.ZERO ) THEN
138 IF( RESID.NE.ZERO )
139 $ RESID = ONE / EPS
140 ELSE
141 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
142 END IF
143 *
144 RETURN
145 *
146 * End of DGET01
147 *
148 END
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 INTEGER LDA, LDAFAC, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET01 reconstructs a matrix A from its L*U factorization and
21 * computes the residual
22 * norm(L*U - A) / ( N * norm(A) * EPS ),
23 * where EPS is the machine epsilon.
24 *
25 * Arguments
26 * ==========
27 *
28 * M (input) INTEGER
29 * The number of rows of the matrix A. M >= 0.
30 *
31 * N (input) INTEGER
32 * The number of columns of the matrix A. N >= 0.
33 *
34 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
35 * The original M x N matrix A.
36 *
37 * LDA (input) INTEGER
38 * The leading dimension of the array A. LDA >= max(1,M).
39 *
40 * AFAC (input/output) DOUBLE PRECISION array, dimension (LDAFAC,N)
41 * The factored form of the matrix A. AFAC contains the factors
42 * L and U from the L*U factorization as computed by DGETRF.
43 * Overwritten with the reconstructed matrix, and then with the
44 * difference L*U - A.
45 *
46 * LDAFAC (input) INTEGER
47 * The leading dimension of the array AFAC. LDAFAC >= max(1,M).
48 *
49 * IPIV (input) INTEGER array, dimension (N)
50 * The pivot indices from DGETRF.
51 *
52 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
53 *
54 * RESID (output) DOUBLE PRECISION
55 * norm(L*U - A) / ( N * norm(A) * EPS )
56 *
57 * =====================================================================
58 *
59 *
60 * .. Parameters ..
61 DOUBLE PRECISION ZERO, ONE
62 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER I, J, K
66 DOUBLE PRECISION ANORM, EPS, T
67 * ..
68 * .. External Functions ..
69 DOUBLE PRECISION DDOT, DLAMCH, DLANGE
70 EXTERNAL DDOT, DLAMCH, DLANGE
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL DGEMV, DLASWP, DSCAL, DTRMV
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE, MIN
77 * ..
78 * .. Executable Statements ..
79 *
80 * Quick exit if M = 0 or N = 0.
81 *
82 IF( M.LE.0 .OR. N.LE.0 ) THEN
83 RESID = ZERO
84 RETURN
85 END IF
86 *
87 * Determine EPS and the norm of A.
88 *
89 EPS = DLAMCH( 'Epsilon' )
90 ANORM = DLANGE( '1', M, N, A, LDA, RWORK )
91 *
92 * Compute the product L*U and overwrite AFAC with the result.
93 * A column at a time of the product is obtained, starting with
94 * column N.
95 *
96 DO 10 K = N, 1, -1
97 IF( K.GT.M ) THEN
98 CALL DTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC,
99 $ LDAFAC, AFAC( 1, K ), 1 )
100 ELSE
101 *
102 * Compute elements (K+1:M,K)
103 *
104 T = AFAC( K, K )
105 IF( K+1.LE.M ) THEN
106 CALL DSCAL( M-K, T, AFAC( K+1, K ), 1 )
107 CALL DGEMV( 'No transpose', M-K, K-1, ONE,
108 $ AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1, ONE,
109 $ AFAC( K+1, K ), 1 )
110 END IF
111 *
112 * Compute the (K,K) element
113 *
114 AFAC( K, K ) = T + DDOT( K-1, AFAC( K, 1 ), LDAFAC,
115 $ AFAC( 1, K ), 1 )
116 *
117 * Compute elements (1:K-1,K)
118 *
119 CALL DTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC,
120 $ LDAFAC, AFAC( 1, K ), 1 )
121 END IF
122 10 CONTINUE
123 CALL DLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 )
124 *
125 * Compute the difference L*U - A and store in AFAC.
126 *
127 DO 30 J = 1, N
128 DO 20 I = 1, M
129 AFAC( I, J ) = AFAC( I, J ) - A( I, J )
130 20 CONTINUE
131 30 CONTINUE
132 *
133 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
134 *
135 RESID = DLANGE( '1', M, N, AFAC, LDAFAC, RWORK )
136 *
137 IF( ANORM.LE.ZERO ) THEN
138 IF( RESID.NE.ZERO )
139 $ RESID = ONE / EPS
140 ELSE
141 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
142 END IF
143 *
144 RETURN
145 *
146 * End of DGET01
147 *
148 END