LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
slqt05.f
Go to the documentation of this file.
1 *> \brief \b SLQT05
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE SLQT05(M,N,L,NB,RESULT)
12 *
13 * .. Scalar Arguments ..
14 * INTEGER LWORK, M, N, L, NB, LDT
15 * .. Return values ..
16 * REAL RESULT(6)
17 *
18 *
19 *> \par Purpose:
20 * =============
21 *>
22 *> \verbatim
23 *> SQRT05 tests STPLQT and STPMLQT.
24 *> \endverbatim
25 *
26 * Arguments:
27 * ==========
28 *
29 *> \param[in] M
30 *> \verbatim
31 *> M is INTEGER
32 *> Number of rows in lower part of the test matrix.
33 *> \endverbatim
34 *>
35 *> \param[in] N
36 *> \verbatim
37 *> N is INTEGER
38 *> Number of columns in test matrix.
39 *> \endverbatim
40 *>
41 *> \param[in] L
42 *> \verbatim
43 *> L is INTEGER
44 *> The number of rows of the upper trapezoidal part the
45 *> lower test matrix. 0 <= L <= M.
46 *> \endverbatim
47 *>
48 *> \param[in] NB
49 *> \verbatim
50 *> NB is INTEGER
51 *> Block size of test matrix. NB <= N.
52 *> \endverbatim
53 *>
54 *> \param[out] RESULT
55 *> \verbatim
56 *> RESULT is REAL array, dimension (6)
57 *> Results of each of the six tests below.
58 *>
59 *> RESULT(1) = | A - Q R |
60 *> RESULT(2) = | I - Q^H Q |
61 *> RESULT(3) = | Q C - Q C |
62 *> RESULT(4) = | Q^H C - Q^H C |
63 *> RESULT(5) = | C Q - C Q |
64 *> RESULT(6) = | C Q^H - C Q^H |
65 *> \endverbatim
66 *
67 * Authors:
68 * ========
69 *
70 *> \author Univ. of Tennessee
71 *> \author Univ. of California Berkeley
72 *> \author Univ. of Colorado Denver
73 *> \author NAG Ltd.
74 *
75 *> \ingroup double_lin
76 *
77 * =====================================================================
78  SUBROUTINE slqt05(M,N,L,NB,RESULT)
79  IMPLICIT NONE
80 *
81 * -- LAPACK test routine --
82 * -- LAPACK is a software package provided by Univ. of Tennessee, --
83 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
84 *
85 * .. Scalar Arguments ..
86  INTEGER LWORK, M, N, L, NB, LDT
87 * .. Return values ..
88  REAL RESULT(6)
89 *
90 * =====================================================================
91 *
92 * ..
93 * .. Local allocatable arrays
94  REAL, ALLOCATABLE :: AF(:,:), Q(:,:),
95  $ R(:,:), RWORK(:), WORK( : ), T(:,:),
96  $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
97 *
98 * .. Parameters ..
99  REAL ONE, ZERO
100  parameter( zero = 0.0, one = 1.0 )
101 * ..
102 * .. Local Scalars ..
103  INTEGER INFO, J, K, N2, NP1,i
104  REAL ANORM, EPS, RESID, CNORM, DNORM
105 * ..
106 * .. Local Arrays ..
107  INTEGER ISEED( 4 )
108 * ..
109 * .. External Functions ..
110  REAL SLAMCH, SLANGE, SLANSY
111  LOGICAL LSAME
112  EXTERNAL slamch, slange, slansy, lsame
113 * ..
114 * .. Data statements ..
115  DATA iseed / 1988, 1989, 1990, 1991 /
116 *
117  eps = slamch( 'Epsilon' )
118  k = m
119  n2 = m+n
120  IF( n.GT.0 ) THEN
121  np1 = m+1
122  ELSE
123  np1 = 1
124  END IF
125  lwork = n2*n2*nb
126 *
127 * Dynamically allocate all arrays
128 *
129  ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
130  $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
131  $ d(m,n2),df(m,n2) )
132 *
133 * Put random stuff into A
134 *
135  ldt=nb
136  CALL slaset( 'Full', m, n2, zero, zero, a, m )
137  CALL slaset( 'Full', nb, m, zero, zero, t, nb )
138  DO j=1,m
139  CALL slarnv( 2, iseed, m-j+1, a( j, j ) )
140  END DO
141  IF( n.GT.0 ) THEN
142  DO j=1,n-l
143  CALL slarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
144  END DO
145  END IF
146  IF( l.GT.0 ) THEN
147  DO j=1,l
148  CALL slarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
149  $ + j - 1 ) )
150  END DO
151  END IF
152 *
153 * Copy the matrix A to the array AF.
154 *
155  CALL slacpy( 'Full', m, n2, a, m, af, m )
156 *
157 * Factor the matrix A in the array AF.
158 *
159  CALL stplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
160 *
161 * Generate the (M+N)-by-(M+N) matrix Q by applying H to I
162 *
163  CALL slaset( 'Full', n2, n2, zero, one, q, n2 )
164  CALL sgemlqt( 'L', 'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
165  $ work, info )
166 *
167 * Copy L
168 *
169  CALL slaset( 'Full', n2, n2, zero, zero, r, n2 )
170  CALL slacpy( 'Lower', m, n2, af, m, r, n2 )
171 *
172 * Compute |L - A*Q*T| / |A| and store in RESULT(1)
173 *
174  CALL sgemm( 'N', 'T', m, n2, n2, -one, a, m, q, n2, one, r, n2)
175  anorm = slange( '1', m, n2, a, m, rwork )
176  resid = slange( '1', m, n2, r, n2, rwork )
177  IF( anorm.GT.zero ) THEN
178  result( 1 ) = resid / (eps*anorm*max(1,n2))
179  ELSE
180  result( 1 ) = zero
181  END IF
182 *
183 * Compute |I - Q*Q'| and store in RESULT(2)
184 *
185  CALL slaset( 'Full', n2, n2, zero, one, r, n2 )
186  CALL ssyrk( 'U', 'N', n2, n2, -one, q, n2, one, r, n2 )
187  resid = slansy( '1', 'Upper', n2, r, n2, rwork )
188  result( 2 ) = resid / (eps*max(1,n2))
189 *
190 * Generate random m-by-n matrix C and a copy CF
191 *
192  CALL slaset( 'Full', n2, m, zero, one, c, n2 )
193  DO j=1,m
194  CALL slarnv( 2, iseed, n2, c( 1, j ) )
195  END DO
196  cnorm = slange( '1', n2, m, c, n2, rwork)
197  CALL slacpy( 'Full', n2, m, c, n2, cf, n2 )
198 *
199 * Apply Q to C as Q*C
200 *
201  CALL stpmlqt( 'L','N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
202  $ cf(np1,1),n2,work,info)
203 *
204 * Compute |Q*C - Q*C| / |C|
205 *
206  CALL sgemm( 'N', 'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
207  resid = slange( '1', n2, m, cf, n2, rwork )
208  IF( cnorm.GT.zero ) THEN
209  result( 3 ) = resid / (eps*max(1,n2)*cnorm)
210  ELSE
211  result( 3 ) = zero
212  END IF
213 
214 *
215 * Copy C into CF again
216 *
217  CALL slacpy( 'Full', n2, m, c, n2, cf, n2 )
218 *
219 * Apply Q to C as QT*C
220 *
221  CALL stpmlqt( 'L','T',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
222  $ cf(np1,1),n2,work,info)
223 *
224 * Compute |QT*C - QT*C| / |C|
225 *
226  CALL sgemm('T','N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
227  resid = slange( '1', n2, m, cf, n2, rwork )
228 
229  IF( cnorm.GT.zero ) THEN
230  result( 4 ) = resid / (eps*max(1,n2)*cnorm)
231  ELSE
232  result( 4 ) = zero
233  END IF
234 *
235 * Generate random m-by-n matrix D and a copy DF
236 *
237  DO j=1,n2
238  CALL slarnv( 2, iseed, m, d( 1, j ) )
239  END DO
240  dnorm = slange( '1', m, n2, d, m, rwork)
241  CALL slacpy( 'Full', m, n2, d, m, df, m )
242 *
243 * Apply Q to D as D*Q
244 *
245  CALL stpmlqt('R','N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
246  $ df(1,np1),m,work,info)
247 *
248 * Compute |D*Q - D*Q| / |D|
249 *
250  CALL sgemm('N','N',m,n2,n2,-one,d,m,q,n2,one,df,m)
251  resid = slange('1',m, n2,df,m,rwork )
252  IF( cnorm.GT.zero ) THEN
253  result( 5 ) = resid / (eps*max(1,n2)*dnorm)
254  ELSE
255  result( 5 ) = zero
256  END IF
257 *
258 * Copy D into DF again
259 *
260  CALL slacpy('Full',m,n2,d,m,df,m )
261 *
262 * Apply Q to D as D*QT
263 *
264  CALL stpmlqt('R','T',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
265  $ df(1,np1),m,work,info)
266 
267 *
268 * Compute |D*QT - D*QT| / |D|
269 *
270  CALL sgemm( 'N', 'T', m, n2, n2, -one, d, m, q, n2, one, df, m )
271  resid = slange( '1', m, n2, df, m, rwork )
272  IF( cnorm.GT.zero ) THEN
273  result( 6 ) = resid / (eps*max(1,n2)*dnorm)
274  ELSE
275  result( 6 ) = zero
276  END IF
277 *
278 * Deallocate all arrays
279 *
280  DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
281  RETURN
282  END
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:97
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slqt05(M, N, L, NB, RESULT)
SLQT05
Definition: slqt05.f:79
subroutine sgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMLQT
Definition: sgemlqt.f:153
subroutine stplqt(M, N, L, MB, A, LDA, B, LDB, T, LDT, WORK, INFO)
STPLQT
Definition: stplqt.f:189
subroutine stpmlqt(SIDE, TRANS, M, N, K, L, MB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
STPMLQT
Definition: stpmlqt.f:214
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187