1 SUBROUTINE DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
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 KL, KU, LDA, LDAFAC, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGBT01 reconstructs a band 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 * The expression L*U - A is computed one column at a time, so A and
26 * AFAC are not modified.
27 *
28 * Arguments
29 * =========
30 *
31 * M (input) INTEGER
32 * The number of rows of the matrix A. M >= 0.
33 *
34 * N (input) INTEGER
35 * The number of columns of the matrix A. N >= 0.
36 *
37 * KL (input) INTEGER
38 * The number of subdiagonals within the band of A. KL >= 0.
39 *
40 * KU (input) INTEGER
41 * The number of superdiagonals within the band of A. KU >= 0.
42 *
43 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
44 * The original matrix A in band storage, stored in rows 1 to
45 * KL+KU+1.
46 *
47 * LDA (input) INTEGER.
48 * The leading dimension of the array A. LDA >= max(1,KL+KU+1).
49 *
50 * AFAC (input) DOUBLE PRECISION array, dimension (LDAFAC,N)
51 * The factored form of the matrix A. AFAC contains the banded
52 * factors L and U from the L*U factorization, as computed by
53 * DGBTRF. U is stored as an upper triangular band matrix with
54 * KL+KU superdiagonals in rows 1 to KL+KU+1, and the
55 * multipliers used during the factorization are stored in rows
56 * KL+KU+2 to 2*KL+KU+1. See DGBTRF for further details.
57 *
58 * LDAFAC (input) INTEGER
59 * The leading dimension of the array AFAC.
60 * LDAFAC >= max(1,2*KL*KU+1).
61 *
62 * IPIV (input) INTEGER array, dimension (min(M,N))
63 * The pivot indices from DGBTRF.
64 *
65 * WORK (workspace) DOUBLE PRECISION array, dimension (2*KL+KU+1)
66 *
67 * RESID (output) DOUBLE PRECISION
68 * norm(L*U - A) / ( N * norm(A) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 * ..
76 * .. Local Scalars ..
77 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
78 DOUBLE PRECISION ANORM, EPS, T
79 * ..
80 * .. External Functions ..
81 DOUBLE PRECISION DASUM, DLAMCH
82 EXTERNAL DASUM, DLAMCH
83 * ..
84 * .. External Subroutines ..
85 EXTERNAL DAXPY, DCOPY
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC DBLE, MAX, MIN
89 * ..
90 * .. Executable Statements ..
91 *
92 * Quick exit if M = 0 or N = 0.
93 *
94 RESID = ZERO
95 IF( M.LE.0 .OR. N.LE.0 )
96 $ RETURN
97 *
98 * Determine EPS and the norm of A.
99 *
100 EPS = DLAMCH( 'Epsilon' )
101 KD = KU + 1
102 ANORM = ZERO
103 DO 10 J = 1, N
104 I1 = MAX( KD+1-J, 1 )
105 I2 = MIN( KD+M-J, KL+KD )
106 IF( I2.GE.I1 )
107 $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) )
108 10 CONTINUE
109 *
110 * Compute one column at a time of L*U - A.
111 *
112 KD = KL + KU + 1
113 DO 40 J = 1, N
114 *
115 * Copy the J-th column of U to WORK.
116 *
117 JU = MIN( KL+KU, J-1 )
118 JL = MIN( KL, M-J )
119 LENJ = MIN( M, J ) - J + JU + 1
120 IF( LENJ.GT.0 ) THEN
121 CALL DCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
122 DO 20 I = LENJ + 1, JU + JL + 1
123 WORK( I ) = ZERO
124 20 CONTINUE
125 *
126 * Multiply by the unit lower triangular matrix L. Note that L
127 * is stored as a product of transformations and permutations.
128 *
129 DO 30 I = MIN( M-1, J ), J - JU, -1
130 IL = MIN( KL, M-I )
131 IF( IL.GT.0 ) THEN
132 IW = I - J + JU + 1
133 T = WORK( IW )
134 CALL DAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
135 $ 1 )
136 IP = IPIV( I )
137 IF( I.NE.IP ) THEN
138 IP = IP - J + JU + 1
139 WORK( IW ) = WORK( IP )
140 WORK( IP ) = T
141 END IF
142 END IF
143 30 CONTINUE
144 *
145 * Subtract the corresponding column of A.
146 *
147 JUA = MIN( JU, KU )
148 IF( JUA+JL+1.GT.0 )
149 $ CALL DAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1,
150 $ WORK( JU+1-JUA ), 1 )
151 *
152 * Compute the 1-norm of the column.
153 *
154 RESID = MAX( RESID, DASUM( JU+JL+1, WORK, 1 ) )
155 END IF
156 40 CONTINUE
157 *
158 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
159 *
160 IF( ANORM.LE.ZERO ) THEN
161 IF( RESID.NE.ZERO )
162 $ RESID = ONE / EPS
163 ELSE
164 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
165 END IF
166 *
167 RETURN
168 *
169 * End of DGBT01
170 *
171 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 KL, KU, LDA, LDAFAC, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGBT01 reconstructs a band 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 * The expression L*U - A is computed one column at a time, so A and
26 * AFAC are not modified.
27 *
28 * Arguments
29 * =========
30 *
31 * M (input) INTEGER
32 * The number of rows of the matrix A. M >= 0.
33 *
34 * N (input) INTEGER
35 * The number of columns of the matrix A. N >= 0.
36 *
37 * KL (input) INTEGER
38 * The number of subdiagonals within the band of A. KL >= 0.
39 *
40 * KU (input) INTEGER
41 * The number of superdiagonals within the band of A. KU >= 0.
42 *
43 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
44 * The original matrix A in band storage, stored in rows 1 to
45 * KL+KU+1.
46 *
47 * LDA (input) INTEGER.
48 * The leading dimension of the array A. LDA >= max(1,KL+KU+1).
49 *
50 * AFAC (input) DOUBLE PRECISION array, dimension (LDAFAC,N)
51 * The factored form of the matrix A. AFAC contains the banded
52 * factors L and U from the L*U factorization, as computed by
53 * DGBTRF. U is stored as an upper triangular band matrix with
54 * KL+KU superdiagonals in rows 1 to KL+KU+1, and the
55 * multipliers used during the factorization are stored in rows
56 * KL+KU+2 to 2*KL+KU+1. See DGBTRF for further details.
57 *
58 * LDAFAC (input) INTEGER
59 * The leading dimension of the array AFAC.
60 * LDAFAC >= max(1,2*KL*KU+1).
61 *
62 * IPIV (input) INTEGER array, dimension (min(M,N))
63 * The pivot indices from DGBTRF.
64 *
65 * WORK (workspace) DOUBLE PRECISION array, dimension (2*KL+KU+1)
66 *
67 * RESID (output) DOUBLE PRECISION
68 * norm(L*U - A) / ( N * norm(A) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 * ..
76 * .. Local Scalars ..
77 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
78 DOUBLE PRECISION ANORM, EPS, T
79 * ..
80 * .. External Functions ..
81 DOUBLE PRECISION DASUM, DLAMCH
82 EXTERNAL DASUM, DLAMCH
83 * ..
84 * .. External Subroutines ..
85 EXTERNAL DAXPY, DCOPY
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC DBLE, MAX, MIN
89 * ..
90 * .. Executable Statements ..
91 *
92 * Quick exit if M = 0 or N = 0.
93 *
94 RESID = ZERO
95 IF( M.LE.0 .OR. N.LE.0 )
96 $ RETURN
97 *
98 * Determine EPS and the norm of A.
99 *
100 EPS = DLAMCH( 'Epsilon' )
101 KD = KU + 1
102 ANORM = ZERO
103 DO 10 J = 1, N
104 I1 = MAX( KD+1-J, 1 )
105 I2 = MIN( KD+M-J, KL+KD )
106 IF( I2.GE.I1 )
107 $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) )
108 10 CONTINUE
109 *
110 * Compute one column at a time of L*U - A.
111 *
112 KD = KL + KU + 1
113 DO 40 J = 1, N
114 *
115 * Copy the J-th column of U to WORK.
116 *
117 JU = MIN( KL+KU, J-1 )
118 JL = MIN( KL, M-J )
119 LENJ = MIN( M, J ) - J + JU + 1
120 IF( LENJ.GT.0 ) THEN
121 CALL DCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
122 DO 20 I = LENJ + 1, JU + JL + 1
123 WORK( I ) = ZERO
124 20 CONTINUE
125 *
126 * Multiply by the unit lower triangular matrix L. Note that L
127 * is stored as a product of transformations and permutations.
128 *
129 DO 30 I = MIN( M-1, J ), J - JU, -1
130 IL = MIN( KL, M-I )
131 IF( IL.GT.0 ) THEN
132 IW = I - J + JU + 1
133 T = WORK( IW )
134 CALL DAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
135 $ 1 )
136 IP = IPIV( I )
137 IF( I.NE.IP ) THEN
138 IP = IP - J + JU + 1
139 WORK( IW ) = WORK( IP )
140 WORK( IP ) = T
141 END IF
142 END IF
143 30 CONTINUE
144 *
145 * Subtract the corresponding column of A.
146 *
147 JUA = MIN( JU, KU )
148 IF( JUA+JL+1.GT.0 )
149 $ CALL DAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1,
150 $ WORK( JU+1-JUA ), 1 )
151 *
152 * Compute the 1-norm of the column.
153 *
154 RESID = MAX( RESID, DASUM( JU+JL+1, WORK, 1 ) )
155 END IF
156 40 CONTINUE
157 *
158 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
159 *
160 IF( ANORM.LE.ZERO ) THEN
161 IF( RESID.NE.ZERO )
162 $ RESID = ONE / EPS
163 ELSE
164 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
165 END IF
166 *
167 RETURN
168 *
169 * End of DGBT01
170 *
171 END