1 SUBROUTINE ZGBT01( 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 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGBT01 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) COMPLEX*16 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) COMPLEX*16 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 * ZGBTRF. 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 ZGBTRF 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 ZGBTRF.
64 *
65 * WORK (workspace) COMPLEX*16 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
79 COMPLEX*16 T
80 * ..
81 * .. External Functions ..
82 DOUBLE PRECISION DLAMCH, DZASUM
83 EXTERNAL DLAMCH, DZASUM
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL ZAXPY, ZCOPY
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC DBLE, DCMPLX, MAX, MIN
90 * ..
91 * .. Executable Statements ..
92 *
93 * Quick exit if M = 0 or N = 0.
94 *
95 RESID = ZERO
96 IF( M.LE.0 .OR. N.LE.0 )
97 $ RETURN
98 *
99 * Determine EPS and the norm of A.
100 *
101 EPS = DLAMCH( 'Epsilon' )
102 KD = KU + 1
103 ANORM = ZERO
104 DO 10 J = 1, N
105 I1 = MAX( KD+1-J, 1 )
106 I2 = MIN( KD+M-J, KL+KD )
107 IF( I2.GE.I1 )
108 $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
109 10 CONTINUE
110 *
111 * Compute one column at a time of L*U - A.
112 *
113 KD = KL + KU + 1
114 DO 40 J = 1, N
115 *
116 * Copy the J-th column of U to WORK.
117 *
118 JU = MIN( KL+KU, J-1 )
119 JL = MIN( KL, M-J )
120 LENJ = MIN( M, J ) - J + JU + 1
121 IF( LENJ.GT.0 ) THEN
122 CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
123 DO 20 I = LENJ + 1, JU + JL + 1
124 WORK( I ) = ZERO
125 20 CONTINUE
126 *
127 * Multiply by the unit lower triangular matrix L. Note that L
128 * is stored as a product of transformations and permutations.
129 *
130 DO 30 I = MIN( M-1, J ), J - JU, -1
131 IL = MIN( KL, M-I )
132 IF( IL.GT.0 ) THEN
133 IW = I - J + JU + 1
134 T = WORK( IW )
135 CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
136 $ 1 )
137 IP = IPIV( I )
138 IF( I.NE.IP ) THEN
139 IP = IP - J + JU + 1
140 WORK( IW ) = WORK( IP )
141 WORK( IP ) = T
142 END IF
143 END IF
144 30 CONTINUE
145 *
146 * Subtract the corresponding column of A.
147 *
148 JUA = MIN( JU, KU )
149 IF( JUA+JL+1.GT.0 )
150 $ CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ),
151 $ 1, WORK( JU+1-JUA ), 1 )
152 *
153 * Compute the 1-norm of the column.
154 *
155 RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) )
156 END IF
157 40 CONTINUE
158 *
159 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
160 *
161 IF( ANORM.LE.ZERO ) THEN
162 IF( RESID.NE.ZERO )
163 $ RESID = ONE / EPS
164 ELSE
165 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
166 END IF
167 *
168 RETURN
169 *
170 * End of ZGBT01
171 *
172 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 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGBT01 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) COMPLEX*16 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) COMPLEX*16 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 * ZGBTRF. 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 ZGBTRF 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 ZGBTRF.
64 *
65 * WORK (workspace) COMPLEX*16 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
79 COMPLEX*16 T
80 * ..
81 * .. External Functions ..
82 DOUBLE PRECISION DLAMCH, DZASUM
83 EXTERNAL DLAMCH, DZASUM
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL ZAXPY, ZCOPY
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC DBLE, DCMPLX, MAX, MIN
90 * ..
91 * .. Executable Statements ..
92 *
93 * Quick exit if M = 0 or N = 0.
94 *
95 RESID = ZERO
96 IF( M.LE.0 .OR. N.LE.0 )
97 $ RETURN
98 *
99 * Determine EPS and the norm of A.
100 *
101 EPS = DLAMCH( 'Epsilon' )
102 KD = KU + 1
103 ANORM = ZERO
104 DO 10 J = 1, N
105 I1 = MAX( KD+1-J, 1 )
106 I2 = MIN( KD+M-J, KL+KD )
107 IF( I2.GE.I1 )
108 $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
109 10 CONTINUE
110 *
111 * Compute one column at a time of L*U - A.
112 *
113 KD = KL + KU + 1
114 DO 40 J = 1, N
115 *
116 * Copy the J-th column of U to WORK.
117 *
118 JU = MIN( KL+KU, J-1 )
119 JL = MIN( KL, M-J )
120 LENJ = MIN( M, J ) - J + JU + 1
121 IF( LENJ.GT.0 ) THEN
122 CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
123 DO 20 I = LENJ + 1, JU + JL + 1
124 WORK( I ) = ZERO
125 20 CONTINUE
126 *
127 * Multiply by the unit lower triangular matrix L. Note that L
128 * is stored as a product of transformations and permutations.
129 *
130 DO 30 I = MIN( M-1, J ), J - JU, -1
131 IL = MIN( KL, M-I )
132 IF( IL.GT.0 ) THEN
133 IW = I - J + JU + 1
134 T = WORK( IW )
135 CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
136 $ 1 )
137 IP = IPIV( I )
138 IF( I.NE.IP ) THEN
139 IP = IP - J + JU + 1
140 WORK( IW ) = WORK( IP )
141 WORK( IP ) = T
142 END IF
143 END IF
144 30 CONTINUE
145 *
146 * Subtract the corresponding column of A.
147 *
148 JUA = MIN( JU, KU )
149 IF( JUA+JL+1.GT.0 )
150 $ CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ),
151 $ 1, WORK( JU+1-JUA ), 1 )
152 *
153 * Compute the 1-norm of the column.
154 *
155 RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) )
156 END IF
157 40 CONTINUE
158 *
159 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
160 *
161 IF( ANORM.LE.ZERO ) THEN
162 IF( RESID.NE.ZERO )
163 $ RESID = ONE / EPS
164 ELSE
165 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
166 END IF
167 *
168 RETURN
169 *
170 * End of ZGBT01
171 *
172 END