91 COMPLEX,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ l(:,:), rwork(:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
98 parameter( zero = 0.0)
99 parameter( one = (1.0,0.0), czero=(0.0,0.0) )
102 INTEGER info, j, k, ll, lwork, ldt
103 REAL anorm, eps, resid, cnorm, dnorm
118 DATA iseed / 1988, 1989, 1990, 1991 /
123 lwork = max(2,ll)*max(2,ll)*nb
127 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
128 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
135 CALL clarnv( 2, iseed, m, a( 1, j ) )
137 CALL clacpy(
'Full', m, n, a, m, af, m )
141 CALL cgelqt( m, n, nb, af, m, t, ldt, work, info )
145 CALL claset(
'Full', n, n, czero, one, q, n )
146 CALL cgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
151 CALL claset(
'Full', ll, n, czero, czero, l, ll )
152 CALL clacpy(
'Lower', m, n, af, m, l, ll )
156 CALL cgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, l, ll )
157 anorm =
clange(
'1', m, n, a, m, rwork )
158 resid =
clange(
'1', m, n, l, ll, rwork )
159 IF( anorm.GT.zero )
THEN
160 result( 1 ) = resid / (eps*max(1,m)*anorm)
167 CALL claset(
'Full', n, n, czero, one, l, ll )
168 CALL cherk(
'U',
'C', n, n,
REAL(-ONE), q, n,
REAL(ONE), l, ll)
169 resid =
clansy(
'1',
'Upper', n, l, ll, rwork )
170 result( 2 ) = resid / (eps*max(1,n))
175 CALL clarnv( 2, iseed, n, d( 1, j ) )
177 dnorm =
clange(
'1', n, m, d, n, rwork)
178 CALL clacpy(
'Full', n, m, d, n, df, n )
182 CALL cgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
187 CALL cgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
188 resid =
clange(
'1', n, m, df, n, rwork )
189 IF( dnorm.GT.zero )
THEN
190 result( 3 ) = resid / (eps*max(1,m)*dnorm)
197 CALL clacpy(
'Full', n, m, d, n, df, n )
201 CALL cgemlqt(
'L',
'C', n, m, k, nb, af, m, t, nb, df, n,
206 CALL cgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
207 resid =
clange(
'1', n, m, df, n, rwork )
208 IF( dnorm.GT.zero )
THEN
209 result( 4 ) = resid / (eps*max(1,m)*dnorm)
217 CALL clarnv( 2, iseed, m, c( 1, j ) )
219 cnorm =
clange(
'1', m, n, c, m, rwork)
220 CALL clacpy(
'Full', m, n, c, m, cf, m )
224 CALL cgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
229 CALL cgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
230 resid =
clange(
'1', n, m, df, n, rwork )
231 IF( cnorm.GT.zero )
THEN
232 result( 5 ) = resid / (eps*max(1,m)*dnorm)
239 CALL clacpy(
'Full', m, n, c, m, cf, m )
243 CALL cgemlqt(
'R',
'C', m, n, k, nb, af, m, t, nb, cf, m,
248 CALL cgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
249 resid =
clange(
'1', m, n, cf, m, rwork )
250 IF( cnorm.GT.zero )
THEN
251 result( 6 ) = resid / (eps*max(1,m)*dnorm)
258 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
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...
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
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 ...
subroutine cgelqt(M, N, MB, A, LDA, T, LDT, WORK, INFO)
subroutine cgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
logical function lsame(CA, CB)
LSAME
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
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, or the element of largest absolute value of a complex symmetric matrix.
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM