73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB
81
82 REAL RESULT(6)
83
84
85
86
87
88 COMPLEX, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ L(:,:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
91 REAL, ALLOCATABLE :: RWORK(:)
92
93
94 REAL ZERO
95 COMPLEX ONE, CZERO
96 parameter( zero = 0.0)
97 parameter( one = (1.0,0.0), czero=(0.0,0.0) )
98
99
100 INTEGER INFO, J, K, LL, LWORK, LDT
101 REAL ANORM, EPS, RESID, CNORM, DNORM
102
103
104 INTEGER ISEED( 4 )
105
106
107 REAL SLAMCH
108 REAL CLANGE, CLANSY
109 LOGICAL LSAME
111
112
113 INTRINSIC max, min
114
115
116 DATA iseed / 1988, 1989, 1990, 1991 /
117
119 k = min(m,n)
120 ll = max(m,n)
121 lwork = max(2,ll)*max(2,ll)*nb
122
123
124
125 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
126 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
127 $ d(n,m), df(n,m) )
128
129
130
131 ldt=nb
132 DO j=1,n
133 CALL clarnv( 2, iseed, m, a( 1, j ) )
134 END DO
135 CALL clacpy(
'Full', m, n, a, m, af, m )
136
137
138
139 CALL cgelqt( m, n, nb, af, m, t, ldt, work, info )
140
141
142
143 CALL claset(
'Full', n, n, czero, one, q, n )
144 CALL cgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
145 $ work, info )
146
147
148
149 CALL claset(
'Full', ll, n, czero, czero, l, ll )
150 CALL clacpy(
'Lower', m, n, af, m, l, ll )
151
152
153
154 CALL cgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, l, ll )
155 anorm =
clange(
'1', m, n, a, m, rwork )
156 resid =
clange(
'1', m, n, l, ll, rwork )
157 IF( anorm.GT.zero ) THEN
158 result( 1 ) = resid / (eps*max(1,m)*anorm)
159 ELSE
160 result( 1 ) = zero
161 END IF
162
163
164
165 CALL claset(
'Full', n, n, czero, one, l, ll )
166 CALL cherk(
'U',
'C', n, n, real(-one), q, n, real(one), l, ll)
167 resid =
clansy(
'1',
'Upper', n, l, ll, rwork )
168 result( 2 ) = resid / (eps*max(1,n))
169
170
171
172 DO j=1,m
173 CALL clarnv( 2, iseed, n, d( 1, j ) )
174 END DO
175 dnorm =
clange(
'1', n, m, d, n, rwork)
176 CALL clacpy(
'Full', n, m, d, n, df, n )
177
178
179
180 CALL cgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
181 $ work, info)
182
183
184
185 CALL cgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
186 resid =
clange(
'1', n, m, df, n, rwork )
187 IF( dnorm.GT.zero ) THEN
188 result( 3 ) = resid / (eps*max(1,m)*dnorm)
189 ELSE
190 result( 3 ) = zero
191 END IF
192
193
194
195 CALL clacpy(
'Full', n, m, d, n, df, n )
196
197
198
199 CALL cgemlqt(
'L',
'C', n, m, k, nb, af, m, t, nb, df, n,
200 $ work, info)
201
202
203
204 CALL cgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
205 resid =
clange(
'1', n, m, df, n, rwork )
206 IF( dnorm.GT.zero ) THEN
207 result( 4 ) = resid / (eps*max(1,m)*dnorm)
208 ELSE
209 result( 4 ) = zero
210 END IF
211
212
213
214 DO j=1,n
215 CALL clarnv( 2, iseed, m, c( 1, j ) )
216 END DO
217 cnorm =
clange(
'1', m, n, c, m, rwork)
218 CALL clacpy(
'Full', m, n, c, m, cf, m )
219
220
221
222 CALL cgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
223 $ work, info)
224
225
226
227 CALL cgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
228 resid =
clange(
'1', n, m, df, n, rwork )
229 IF( cnorm.GT.zero ) THEN
230 result( 5 ) = resid / (eps*max(1,m)*dnorm)
231 ELSE
232 result( 5 ) = zero
233 END IF
234
235
236
237 CALL clacpy(
'Full', m, n, c, m, cf, m )
238
239
240
241 CALL cgemlqt(
'R',
'C', m, n, k, nb, af, m, t, nb, cf, m,
242 $ work, info)
243
244
245
246 CALL cgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
247 resid =
clange(
'1', m, n, cf, m, rwork )
248 IF( cnorm.GT.zero ) THEN
249 result( 6 ) = resid / (eps*max(1,m)*dnorm)
250 ELSE
251 result( 6 ) = zero
252 END IF
253
254
255
256 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
257
258 RETURN
subroutine cgelqt(m, n, mb, a, lda, t, ldt, work, info)
CGELQT
subroutine cgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
CGEMLQT
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function clansy(norm, uplo, n, a, lda, work)
CLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME