80 IMPLICIT NONE
81
82
83
84
85
86
87 INTEGER LWORK, M, N, L, NB, LDT
88
89 REAL RESULT(6)
90
91
92
93
94
95 COMPLEX, ALLOCATABLE :: AF(:,:), Q(:,:),
96 $ R(:,:), WORK( : ), T(:,:),
97 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
98 REAL, ALLOCATABLE :: RWORK(:)
99
100
101 REAL ZERO
102 COMPLEX ONE, CZERO
103 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
104
105
106 INTEGER INFO, J, K, M2, NP1
107 REAL ANORM, EPS, RESID, CNORM, DNORM
108
109
110 INTEGER ISEED( 4 )
111
112
113 REAL SLAMCH
114 REAL CLANGE, CLANSY
115 LOGICAL LSAME
117
118
119 DATA iseed / 1988, 1989, 1990, 1991 /
120
122 k = n
123 m2 = m+n
124 IF( m.GT.0 ) THEN
125 np1 = n+1
126 ELSE
127 np1 = 1
128 END IF
129 lwork = m2*m2*nb
130
131
132
133 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
134 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
135 $ d(n,m2),df(n,m2) )
136
137
138
139 ldt=nb
140 CALL claset(
'Full', m2, n, czero, czero, a, m2 )
141 CALL claset(
'Full', nb, n, czero, czero, t, nb )
142 DO j=1,n
143 CALL clarnv( 2, iseed, j, a( 1, j ) )
144 END DO
145 IF( m.GT.0 ) THEN
146 DO j=1,n
147 CALL clarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
148 END DO
149 END IF
150 IF( l.GT.0 ) THEN
151 DO j=1,n
152 CALL clarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
153 END DO
154 END IF
155
156
157
158 CALL clacpy(
'Full', m2, n, a, m2, af, m2 )
159
160
161
162 CALL ctpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
163
164
165
166 CALL claset(
'Full', m2, m2, czero, one, q, m2 )
167 CALL cgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
168 $ work, info )
169
170
171
172 CALL claset(
'Full', m2, n, czero, czero, r, m2 )
173 CALL clacpy(
'Upper', m2, n, af, m2, r, m2 )
174
175
176
177 CALL cgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
178 anorm =
clange(
'1', m2, n, a, m2, rwork )
179 resid =
clange(
'1', m2, n, r, m2, rwork )
180 IF( anorm.GT.zero ) THEN
181 result( 1 ) = resid / (eps*anorm*max(1,m2))
182 ELSE
183 result( 1 ) = zero
184 END IF
185
186
187
188 CALL claset(
'Full', m2, m2, czero, one, r, m2 )
189 CALL cherk(
'U',
'C', m2, m2, real(-one), q, m2, real(one),
190 $ r, m2 )
191 resid =
clansy(
'1',
'Upper', m2, r, m2, rwork )
192 result( 2 ) = resid / (eps*max(1,m2))
193
194
195
196 DO j=1,n
197 CALL clarnv( 2, iseed, m2, c( 1, j ) )
198 END DO
199 cnorm =
clange(
'1', m2, n, c, m2, rwork)
200 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
201
202
203
204 CALL ctpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
205 $ cf(np1,1),m2,work,info)
206
207
208
209 CALL cgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
210 resid =
clange(
'1', m2, n, cf, m2, rwork )
211 IF( cnorm.GT.zero ) THEN
212 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
213 ELSE
214 result( 3 ) = zero
215 END IF
216
217
218
219 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
220
221
222
223 CALL ctpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
224 $ cf(np1,1),m2,work,info)
225
226
227
228 CALL cgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
229 resid =
clange(
'1', m2, n, cf, m2, rwork )
230 IF( cnorm.GT.zero ) THEN
231 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
232 ELSE
233 result( 4 ) = zero
234 END IF
235
236
237
238 DO j=1,m2
239 CALL clarnv( 2, iseed, n, d( 1, j ) )
240 END DO
241 dnorm =
clange(
'1', n, m2, d, n, rwork)
242 CALL clacpy(
'Full', n, m2, d, n, df, n )
243
244
245
246 CALL ctpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
247 $ df(1,np1),n,work,info)
248
249
250
251 CALL cgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
252 resid =
clange(
'1',n, m2,df,n,rwork )
253 IF( cnorm.GT.zero ) THEN
254 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
255 ELSE
256 result( 5 ) = zero
257 END IF
258
259
260
261 CALL clacpy(
'Full',n,m2,d,n,df,n )
262
263
264
265 CALL ctpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
266 $ df(1,np1),n,work,info)
267
268
269
270
271 CALL cgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
272 resid =
clange(
'1', n, m2, df, n, rwork )
273 IF( cnorm.GT.zero ) THEN
274 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
275 ELSE
276 result( 6 ) = zero
277 END IF
278
279
280
281 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
282 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 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
subroutine ctpmqrt(side, trans, m, n, k, l, nb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
CTPMQRT
subroutine ctpqrt(m, n, l, nb, a, lda, b, ldb, t, ldt, work, info)
CTPQRT