80 IMPLICIT NONE
81
82
83
84
85
86
87 INTEGER LWORK, M, N, L, NB, LDT
88
89 DOUBLE PRECISION RESULT(6)
90
91
92
93
94
95 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
96 $ R(:,:), WORK( : ), T(:,:),
97 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
98 DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
99
100
101 DOUBLE PRECISION ZERO
102 COMPLEX*16 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 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
108
109
110 INTEGER ISEED( 4 )
111
112
113 DOUBLE PRECISION DLAMCH
114 DOUBLE PRECISION ZLANGE, ZLANSY
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 zlaset(
'Full', m, n2, czero, czero, a, m )
141 CALL zlaset(
'Full', nb, m, czero, czero, t, nb )
142 DO j=1,m
143 CALL zlarnv( 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 zlarnv( 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 zlarnv( 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 zlacpy(
'Full', m, n2, a, m, af, m )
160
161
162
163 CALL ztplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
164
165
166
167 CALL zlaset(
'Full', n2, n2, czero, one, q, n2 )
168 CALL zgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
169 $ work, info )
170
171
172
173 CALL zlaset(
'Full', n2, n2, czero, czero, r, n2 )
174 CALL zlacpy(
'Lower', m, n2, af, m, r, n2 )
175
176
177
178 CALL zgemm(
'N',
'C', m, n2, n2, -one, a, m, q, n2, one, r, n2)
179 anorm =
zlange(
'1', m, n2, a, m, rwork )
180 resid =
zlange(
'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 zlaset(
'Full', n2, n2, czero, one, r, n2 )
190 CALL zherk(
'U',
'N', n2, n2, dreal(-one), q, n2, dreal(one),
191 $ r, n2 )
192 resid =
zlansy(
'1',
'Upper', n2, r, n2, rwork )
193 result( 2 ) = resid / (eps*max(1,n2))
194
195
196
197 CALL zlaset(
'Full', n2, m, czero, one, c, n2 )
198 DO j=1,m
199 CALL zlarnv( 2, iseed, n2, c( 1, j ) )
200 END DO
201 cnorm =
zlange(
'1', n2, m, c, n2, rwork)
202 CALL zlacpy(
'Full', n2, m, c, n2, cf, n2 )
203
204
205
206 CALL ztpmlqt(
'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 zgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
212 resid =
zlange(
'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 zlacpy(
'Full', n2, m, c, n2, cf, n2 )
223
224
225
226 CALL ztpmlqt(
'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 zgemm(
'C',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
232 resid =
zlange(
'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 zlarnv( 2, iseed, m, d( 1, j ) )
244 END DO
245 dnorm =
zlange(
'1', m, n2, d, m, rwork)
246 CALL zlacpy(
'Full', m, n2, d, m, df, m )
247
248
249
250 CALL ztpmlqt(
'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 zgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
256 resid =
zlange(
'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 zlacpy(
'Full',m,n2,d,m,df,m )
266
267
268
269 CALL ztpmlqt(
'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 zgemm(
'N',
'C', m, n2, n2, -one, d, m, q, n2, one, df, m )
276 resid =
zlange(
'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
double precision function dlamch(CMACH)
DLAMCH
logical function lsame(CA, CB)
LSAME
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
double precision function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine zgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMLQT
subroutine ztplqt(M, N, L, MB, A, LDA, B, LDB, T, LDT, WORK, INFO)
ZTPLQT
subroutine ztpmlqt(SIDE, TRANS, M, N, K, L, MB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
ZTPMLQT