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, N2, NP1,i
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 = m
123 n2 = m+n
124 IF( n.GT.0 ) THEN
125 np1 = m+1
126 ELSE
127 np1 = 1
128 END IF
129 lwork = n2*n2*nb
130
131
132
133 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
134 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
135 $ d(m,n2),df(m,n2) )
136
137
138
139 ldt=nb
140 CALL claset(
'Full', m, n2, czero, czero, a, m )
141 CALL claset(
'Full', nb, m, czero, czero, t, nb )
142 DO j=1,m
143 CALL clarnv( 2, iseed, m-j+1, a( j, j ) )
144 END DO
145 IF( n.GT.0 ) THEN
146 DO j=1,n-l
147 CALL clarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
148 END DO
149 END IF
150 IF( l.GT.0 ) THEN
151 DO j=1,l
152 CALL clarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
153 $ + j - 1 ) )
154 END DO
155 END IF
156
157
158
159 CALL clacpy(
'Full', m, n2, a, m, af, m )
160
161
162
163 CALL ctplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
164
165
166
167 CALL claset(
'Full', n2, n2, czero, one, q, n2 )
168 CALL cgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
169 $ work, info )
170
171
172
173 CALL claset(
'Full', n2, n2, czero, czero, r, n2 )
174 CALL clacpy(
'Lower', m, n2, af, m, r, n2 )
175
176
177
178 CALL cgemm(
'N',
'C', m, n2, n2, -one, a, m, q, n2, one, r, n2)
179 anorm =
clange(
'1', m, n2, a, m, rwork )
180 resid =
clange(
'1', m, n2, r, n2, rwork )
181 IF( anorm.GT.zero ) THEN
182 result( 1 ) = resid / (eps*anorm*max(1,n2))
183 ELSE
184 result( 1 ) = zero
185 END IF
186
187
188
189 CALL claset(
'Full', n2, n2, czero, one, r, n2 )
190 CALL cherk(
'U',
'N', n2, n2, real(-one), q, n2, real(one),
191 $ r, n2 )
192 resid =
clansy(
'1',
'Upper', n2, r, n2, rwork )
193 result( 2 ) = resid / (eps*max(1,n2))
194
195
196
197 CALL claset(
'Full', n2, m, czero, one, c, n2 )
198 DO j=1,m
199 CALL clarnv( 2, iseed, n2, c( 1, j ) )
200 END DO
201 cnorm =
clange(
'1', n2, m, c, n2, rwork)
202 CALL clacpy(
'Full', n2, m, c, n2, cf, n2 )
203
204
205
206 CALL ctpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
207 $ cf(np1,1),n2,work,info)
208
209
210
211 CALL cgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
212 resid =
clange(
'1', n2, m, cf, n2, rwork )
213 IF( cnorm.GT.zero ) THEN
214 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
215 ELSE
216 result( 3 ) = zero
217 END IF
218
219
220
221
222 CALL clacpy(
'Full', n2, m, c, n2, cf, n2 )
223
224
225
226 CALL ctpmlqt(
'L',
'C',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
227 $ cf(np1,1),n2,work,info)
228
229
230
231 CALL cgemm(
'C',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
232 resid =
clange(
'1', n2, m, cf, n2, rwork )
233
234 IF( cnorm.GT.zero ) THEN
235 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
236 ELSE
237 result( 4 ) = zero
238 END IF
239
240
241
242 DO j=1,n2
243 CALL clarnv( 2, iseed, m, d( 1, j ) )
244 END DO
245 dnorm =
clange(
'1', m, n2, d, m, rwork)
246 CALL clacpy(
'Full', m, n2, d, m, df, m )
247
248
249
250 CALL ctpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
251 $ df(1,np1),m,work,info)
252
253
254
255 CALL cgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
256 resid =
clange(
'1',m, n2,df,m,rwork )
257 IF( cnorm.GT.zero ) THEN
258 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
259 ELSE
260 result( 5 ) = zero
261 END IF
262
263
264
265 CALL clacpy(
'Full',m,n2,d,m,df,m )
266
267
268
269 CALL ctpmlqt(
'R',
'C',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
270 $ df(1,np1),m,work,info)
271
272
273
274
275 CALL cgemm(
'N',
'C', m, n2, n2, -one, d, m, q, n2, one, df, m )
276 resid =
clange(
'1', m, n2, df, m, rwork )
277 IF( cnorm.GT.zero ) THEN
278 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
279 ELSE
280 result( 6 ) = zero
281 END IF
282
283
284
285 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
286 RETURN
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
subroutine ctplqt(m, n, l, mb, a, lda, b, ldb, t, ldt, work, info)
CTPLQT
subroutine ctpmlqt(side, trans, m, n, k, l, mb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
CTPMLQT