73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB, LDT
81
82 REAL RESULT(6)
83
84
85
86
87
88 COMPLEX, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ R(:,:), 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, one = (1.0,0.0), czero=(0.0,0.0) )
97
98
99 INTEGER INFO, J, K, L, LWORK
100 REAL ANORM, EPS, RESID, CNORM, DNORM
101
102
103 INTEGER ISEED( 4 )
104
105
106 REAL SLAMCH
107 REAL CLANGE, CLANSY
108 LOGICAL LSAME
110
111
112 INTRINSIC max, min
113
114
115 DATA iseed / 1988, 1989, 1990, 1991 /
116
118 k = min(m,n)
119 l = max(m,n)
120 lwork = max(2,l)*max(2,l)*nb
121
122
123
124 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
125 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
126 $ d(n,m), df(n,m) )
127
128
129
130 ldt=nb
131 DO j=1,n
132 CALL clarnv( 2, iseed, m, a( 1, j ) )
133 END DO
134 CALL clacpy(
'Full', m, n, a, m, af, m )
135
136
137
138 CALL cgeqrt( m, n, nb, af, m, t, ldt, work, info )
139
140
141
142 CALL claset(
'Full', m, m, czero, one, q, m )
143 CALL cgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
144 $ work, info )
145
146
147
148 CALL claset(
'Full', m, n, czero, czero, r, m )
149 CALL clacpy(
'Upper', m, n, af, m, r, m )
150
151
152
153 CALL cgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
154 anorm =
clange(
'1', m, n, a, m, rwork )
155 resid =
clange(
'1', m, n, r, m, rwork )
156 IF( anorm.GT.zero ) THEN
157 result( 1 ) = resid / (eps*max(1,m)*anorm)
158 ELSE
159 result( 1 ) = zero
160 END IF
161
162
163
164 CALL claset(
'Full', m, m, czero, one, r, m )
165 CALL cherk(
'U',
'C', m, m, real(-one), q, m, real(one), r, m )
166 resid =
clansy(
'1',
'Upper', m, r, m, rwork )
167 result( 2 ) = resid / (eps*max(1,m))
168
169
170
171 DO j=1,n
172 CALL clarnv( 2, iseed, m, c( 1, j ) )
173 END DO
174 cnorm =
clange(
'1', m, n, c, m, rwork)
175 CALL clacpy(
'Full', m, n, c, m, cf, m )
176
177
178
179 CALL cgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
180 $ work, info)
181
182
183
184 CALL cgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
185 resid =
clange(
'1', m, n, cf, m, rwork )
186 IF( cnorm.GT.zero ) THEN
187 result( 3 ) = resid / (eps*max(1,m)*cnorm)
188 ELSE
189 result( 3 ) = zero
190 END IF
191
192
193
194 CALL clacpy(
'Full', m, n, c, m, cf, m )
195
196
197
198 CALL cgemqrt(
'L',
'C', m, n, k, nb, af, m, t, nb, cf, m,
199 $ work, info)
200
201
202
203 CALL cgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
204 resid =
clange(
'1', m, n, cf, m, rwork )
205 IF( cnorm.GT.zero ) THEN
206 result( 4 ) = resid / (eps*max(1,m)*cnorm)
207 ELSE
208 result( 4 ) = zero
209 END IF
210
211
212
213 DO j=1,m
214 CALL clarnv( 2, iseed, n, d( 1, j ) )
215 END DO
216 dnorm =
clange(
'1', n, m, d, n, rwork)
217 CALL clacpy(
'Full', n, m, d, n, df, n )
218
219
220
221 CALL cgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
222 $ work, info)
223
224
225
226 CALL cgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
227 resid =
clange(
'1', n, m, df, n, rwork )
228 IF( cnorm.GT.zero ) THEN
229 result( 5 ) = resid / (eps*max(1,m)*dnorm)
230 ELSE
231 result( 5 ) = zero
232 END IF
233
234
235
236 CALL clacpy(
'Full', n, m, d, n, df, n )
237
238
239
240 CALL cgemqrt(
'R',
'C', n, m, k, nb, af, m, t, nb, df, n,
241 $ work, info)
242
243
244
245 CALL cgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
246 resid =
clange(
'1', n, m, df, n, rwork )
247 IF( cnorm.GT.zero ) THEN
248 result( 6 ) = resid / (eps*max(1,m)*dnorm)
249 ELSE
250 result( 6 ) = zero
251 END IF
252
253
254
255 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
256
257 RETURN
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
CGEMQRT
subroutine cgeqrt(m, n, nb, a, lda, t, ldt, work, info)
CGEQRT
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