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 DOUBLE PRECISION, ALLOCATABLE :: AF(:,:), Q(:,:),
96 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
97 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
98
99
100 DOUBLE PRECISION ONE, ZERO
101 parameter( zero = 0.0, one = 1.0 )
102
103
104 INTEGER INFO, J, K, N2, NP1,i
105 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
106
107
108 INTEGER ISEED( 4 )
109
110
111 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
112 LOGICAL LSAME
114
115
116 DATA iseed / 1988, 1989, 1990, 1991 /
117
119 k = m
120 n2 = m+n
121 IF( n.GT.0 ) THEN
122 np1 = m+1
123 ELSE
124 np1 = 1
125 END IF
126 lwork = n2*n2*nb
127
128
129
130 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
131 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
132 $ d(m,n2),df(m,n2) )
133
134
135
136 ldt=nb
137 CALL dlaset(
'Full', m, n2, zero, zero, a, m )
138 CALL dlaset(
'Full', nb, m, zero, zero, t, nb )
139 DO j=1,m
140 CALL dlarnv( 2, iseed, m-j+1, a( j, j ) )
141 END DO
142 IF( n.GT.0 ) THEN
143 DO j=1,n-l
144 CALL dlarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
145 END DO
146 END IF
147 IF( l.GT.0 ) THEN
148 DO j=1,l
149 CALL dlarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
150 $ + j - 1 ) )
151 END DO
152 END IF
153
154
155
156 CALL dlacpy(
'Full', m, n2, a, m, af, m )
157
158
159
160 CALL dtplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
161
162
163
164 CALL dlaset(
'Full', n2, n2, zero, one, q, n2 )
165 CALL dgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
166 $ work, info )
167
168
169
170 CALL dlaset(
'Full', n2, n2, zero, zero, r, n2 )
171 CALL dlacpy(
'Lower', m, n2, af, m, r, n2 )
172
173
174
175 CALL dgemm(
'N',
'T', m, n2, n2, -one, a, m, q, n2, one, r, n2)
176 anorm =
dlange(
'1', m, n2, a, m, rwork )
177 resid =
dlange(
'1', m, n2, r, n2, rwork )
178 IF( anorm.GT.zero ) THEN
179 result( 1 ) = resid / (eps*anorm*max(1,n2))
180 ELSE
181 result( 1 ) = zero
182 END IF
183
184
185
186 CALL dlaset(
'Full', n2, n2, zero, one, r, n2 )
187 CALL dsyrk(
'U',
'N', n2, n2, -one, q, n2, one, r, n2 )
188 resid =
dlansy(
'1',
'Upper', n2, r, n2, rwork )
189 result( 2 ) = resid / (eps*max(1,n2))
190
191
192
193 CALL dlaset(
'Full', n2, m, zero, one, c, n2 )
194 DO j=1,m
195 CALL dlarnv( 2, iseed, n2, c( 1, j ) )
196 END DO
197 cnorm =
dlange(
'1', n2, m, c, n2, rwork)
198 CALL dlacpy(
'Full', n2, m, c, n2, cf, n2 )
199
200
201
202 CALL dtpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
203 $ cf(np1,1),n2,work,info)
204
205
206
207 CALL dgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
208 resid =
dlange(
'1', n2, m, cf, n2, rwork )
209 IF( cnorm.GT.zero ) THEN
210 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
211 ELSE
212 result( 3 ) = zero
213 END IF
214
215
216
217
218 CALL dlacpy(
'Full', n2, m, c, n2, cf, n2 )
219
220
221
222 CALL dtpmlqt(
'L',
'T',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
223 $ cf(np1,1),n2,work,info)
224
225
226
227 CALL dgemm(
'T',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
228 resid =
dlange(
'1', n2, m, cf, n2, rwork )
229
230 IF( cnorm.GT.zero ) THEN
231 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
232 ELSE
233 result( 4 ) = zero
234 END IF
235
236
237
238 DO j=1,n2
239 CALL dlarnv( 2, iseed, m, d( 1, j ) )
240 END DO
241 dnorm =
dlange(
'1', m, n2, d, m, rwork)
242 CALL dlacpy(
'Full', m, n2, d, m, df, m )
243
244
245
246 CALL dtpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
247 $ df(1,np1),m,work,info)
248
249
250
251 CALL dgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
252 resid =
dlange(
'1',m, n2,df,m,rwork )
253 IF( cnorm.GT.zero ) THEN
254 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
255 ELSE
256 result( 5 ) = zero
257 END IF
258
259
260
261 CALL dlacpy(
'Full',m,n2,d,m,df,m )
262
263
264
265 CALL dtpmlqt(
'R',
'T',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
266 $ df(1,np1),m,work,info)
267
268
269
270
271 CALL dgemm(
'N',
'T', m, n2, n2, -one, d, m, q, n2, one, df, m )
272 resid =
dlange(
'1', m, n2, df, m, rwork )
273 IF( cnorm.GT.zero ) THEN
274 result( 6 ) = resid / (eps*max(1,n2)*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 dgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
DGEMLQT
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
subroutine dtplqt(m, n, l, mb, a, lda, b, ldb, t, ldt, work, info)
DTPLQT
subroutine dtpmlqt(side, trans, m, n, k, l, mb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
DTPMLQT