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