OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
double
LAPACK Testing » Matrix Generation

Functions

subroutine dlagge (m, n, kl, ku, d, a, lda, iseed, work, info)
 DLAGGE
subroutine dlagsy (n, k, d, a, lda, iseed, work, info)
 DLAGSY
subroutine dlahilb (n, nrhs, a, lda, x, ldx, b, ldb, work, info)
 DLAHILB
subroutine dlakf2 (m, n, a, lda, b, d, e, z, ldz)
 DLAKF2
subroutine dlarge (n, a, lda, iseed, work, info)
 DLARGE
double precision function dlarnd (idist, iseed)
 DLARND
subroutine dlaror (side, init, m, n, a, lda, iseed, x, info)
 DLAROR
subroutine dlarot (lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
 DLAROT
subroutine dlatm1 (mode, cond, irsign, idist, iseed, d, n, info)
 DLATM1
double precision function dlatm2 (m, n, i, j, kl, ku, idist, iseed, d, igrade, dl, dr, ipvtng, iwork, sparse)
 DLATM2
double precision function dlatm3 (m, n, i, j, isub, jsub, kl, ku, idist, iseed, d, igrade, dl, dr, ipvtng, iwork, sparse)
 DLATM3
subroutine dlatm5 (prtype, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, r, ldr, l, ldl, alpha, qblcka, qblckb)
 DLATM5
subroutine dlatm6 (type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
 DLATM6
subroutine dlatm7 (mode, cond, irsign, idist, iseed, d, n, rank, info)
 DLATM7
subroutine dlatme (n, dist, iseed, d, mode, cond, dmax, ei, rsign, upper, sim, ds, modes, conds, kl, ku, anorm, a, lda, work, info)
 DLATME
subroutine dlatmr (m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
 DLATMR
subroutine dlatms (m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
 DLATMS
subroutine dlatmt (m, n, dist, iseed, sym, d, mode, cond, dmax, rank, kl, ku, pack, a, lda, work, info)
 DLATMT

Detailed Description

This is the group of double LAPACK TESTING MATGEN routines.

Function Documentation

◆ dlagge()

subroutine dlagge ( integer m,
integer n,
integer kl,
integer ku,
double precision, dimension( * ) d,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) work,
integer info )

DLAGGE

Purpose:
!>
!> DLAGGE generates a real general m by n matrix A, by pre- and post-
!> multiplying a real diagonal matrix D with random orthogonal matrices:
!> A = U*D*V. The lower and upper bandwidths may then be reduced to
!> kl and ku by additional orthogonal transformations.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.
!>
[in]KL
!>          KL is INTEGER
!>          The number of nonzero subdiagonals within the band of A.
!>          0 <= KL <= M-1.
!>
[in]KU
!>          KU is INTEGER
!>          The number of nonzero superdiagonals within the band of A.
!>          0 <= KU <= N-1.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (min(M,N))
!>          The diagonal elements of the diagonal matrix D.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          The generated m by n matrix A.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= M.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry, the seed of the random number generator; the array
!>          elements must be between 0 and 4095, and ISEED(4) must be
!>          odd.
!>          On exit, the seed is updated.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (M+N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 112 of file dlagge.f.

113*
114* -- LAPACK auxiliary routine --
115* -- LAPACK is a software package provided by Univ. of Tennessee, --
116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117*
118* .. Scalar Arguments ..
119 INTEGER INFO, KL, KU, LDA, M, N
120* ..
121* .. Array Arguments ..
122 INTEGER ISEED( 4 )
123 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
124* ..
125*
126* =====================================================================
127*
128* .. Parameters ..
129 DOUBLE PRECISION ZERO, ONE
130 parameter( zero = 0.0d+0, one = 1.0d+0 )
131* ..
132* .. Local Scalars ..
133 INTEGER I, J
134 DOUBLE PRECISION TAU, WA, WB, WN
135* ..
136* .. External Subroutines ..
137 EXTERNAL dgemv, dger, dlarnv, dscal, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, min, sign
141* ..
142* .. External Functions ..
143 DOUBLE PRECISION DNRM2
144 EXTERNAL dnrm2
145* ..
146* .. Executable Statements ..
147*
148* Test the input arguments
149*
150 info = 0
151 IF( m.LT.0 ) THEN
152 info = -1
153 ELSE IF( n.LT.0 ) THEN
154 info = -2
155 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
156 info = -3
157 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
158 info = -4
159 ELSE IF( lda.LT.max( 1, m ) ) THEN
160 info = -7
161 END IF
162 IF( info.LT.0 ) THEN
163 CALL xerbla( 'DLAGGE', -info )
164 RETURN
165 END IF
166*
167* initialize A to diagonal matrix
168*
169 DO 20 j = 1, n
170 DO 10 i = 1, m
171 a( i, j ) = zero
172 10 CONTINUE
173 20 CONTINUE
174 DO 30 i = 1, min( m, n )
175 a( i, i ) = d( i )
176 30 CONTINUE
177*
178* Quick exit if the user wants a diagonal matrix
179*
180 IF(( kl .EQ. 0 ).AND.( ku .EQ. 0)) RETURN
181*
182* pre- and post-multiply A by random orthogonal matrices
183*
184 DO 40 i = min( m, n ), 1, -1
185 IF( i.LT.m ) THEN
186*
187* generate random reflection
188*
189 CALL dlarnv( 3, iseed, m-i+1, work )
190 wn = dnrm2( m-i+1, work, 1 )
191 wa = sign( wn, work( 1 ) )
192 IF( wn.EQ.zero ) THEN
193 tau = zero
194 ELSE
195 wb = work( 1 ) + wa
196 CALL dscal( m-i, one / wb, work( 2 ), 1 )
197 work( 1 ) = one
198 tau = wb / wa
199 END IF
200*
201* multiply A(i:m,i:n) by random reflection from the left
202*
203 CALL dgemv( 'Transpose', m-i+1, n-i+1, one, a( i, i ), lda,
204 $ work, 1, zero, work( m+1 ), 1 )
205 CALL dger( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
206 $ a( i, i ), lda )
207 END IF
208 IF( i.LT.n ) THEN
209*
210* generate random reflection
211*
212 CALL dlarnv( 3, iseed, n-i+1, work )
213 wn = dnrm2( n-i+1, work, 1 )
214 wa = sign( wn, work( 1 ) )
215 IF( wn.EQ.zero ) THEN
216 tau = zero
217 ELSE
218 wb = work( 1 ) + wa
219 CALL dscal( n-i, one / wb, work( 2 ), 1 )
220 work( 1 ) = one
221 tau = wb / wa
222 END IF
223*
224* multiply A(i:m,i:n) by random reflection from the right
225*
226 CALL dgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
227 $ lda, work, 1, zero, work( n+1 ), 1 )
228 CALL dger( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
229 $ a( i, i ), lda )
230 END IF
231 40 CONTINUE
232*
233* Reduce number of subdiagonals to KL and number of superdiagonals
234* to KU
235*
236 DO 70 i = 1, max( m-1-kl, n-1-ku )
237 IF( kl.LE.ku ) THEN
238*
239* annihilate subdiagonal elements first (necessary if KL = 0)
240*
241 IF( i.LE.min( m-1-kl, n ) ) THEN
242*
243* generate reflection to annihilate A(kl+i+1:m,i)
244*
245 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
246 wa = sign( wn, a( kl+i, i ) )
247 IF( wn.EQ.zero ) THEN
248 tau = zero
249 ELSE
250 wb = a( kl+i, i ) + wa
251 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
252 a( kl+i, i ) = one
253 tau = wb / wa
254 END IF
255*
256* apply reflection to A(kl+i:m,i+1:n) from the left
257*
258 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
259 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
260 $ work, 1 )
261 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
262 $ a( kl+i, i+1 ), lda )
263 a( kl+i, i ) = -wa
264 END IF
265*
266 IF( i.LE.min( n-1-ku, m ) ) THEN
267*
268* generate reflection to annihilate A(i,ku+i+1:n)
269*
270 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
271 wa = sign( wn, a( i, ku+i ) )
272 IF( wn.EQ.zero ) THEN
273 tau = zero
274 ELSE
275 wb = a( i, ku+i ) + wa
276 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
277 a( i, ku+i ) = one
278 tau = wb / wa
279 END IF
280*
281* apply reflection to A(i+1:m,ku+i:n) from the right
282*
283 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
284 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
285 $ work, 1 )
286 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
287 $ lda, a( i+1, ku+i ), lda )
288 a( i, ku+i ) = -wa
289 END IF
290 ELSE
291*
292* annihilate superdiagonal elements first (necessary if
293* KU = 0)
294*
295 IF( i.LE.min( n-1-ku, m ) ) THEN
296*
297* generate reflection to annihilate A(i,ku+i+1:n)
298*
299 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
300 wa = sign( wn, a( i, ku+i ) )
301 IF( wn.EQ.zero ) THEN
302 tau = zero
303 ELSE
304 wb = a( i, ku+i ) + wa
305 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
306 a( i, ku+i ) = one
307 tau = wb / wa
308 END IF
309*
310* apply reflection to A(i+1:m,ku+i:n) from the right
311*
312 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
313 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
314 $ work, 1 )
315 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
316 $ lda, a( i+1, ku+i ), lda )
317 a( i, ku+i ) = -wa
318 END IF
319*
320 IF( i.LE.min( m-1-kl, n ) ) THEN
321*
322* generate reflection to annihilate A(kl+i+1:m,i)
323*
324 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
325 wa = sign( wn, a( kl+i, i ) )
326 IF( wn.EQ.zero ) THEN
327 tau = zero
328 ELSE
329 wb = a( kl+i, i ) + wa
330 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
331 a( kl+i, i ) = one
332 tau = wb / wa
333 END IF
334*
335* apply reflection to A(kl+i:m,i+1:n) from the left
336*
337 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
338 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
339 $ work, 1 )
340 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
341 $ a( kl+i, i+1 ), lda )
342 a( kl+i, i ) = -wa
343 END IF
344 END IF
345*
346 IF (i .LE. n) THEN
347 DO 50 j = kl + i + 1, m
348 a( j, i ) = zero
349 50 CONTINUE
350 END IF
351*
352 IF (i .LE. m) THEN
353 DO 60 j = ku + i + 1, n
354 a( i, j ) = zero
355 60 CONTINUE
356 END IF
357 70 CONTINUE
358 RETURN
359*
360* End of DLAGGE
361*
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ dlagsy()

subroutine dlagsy ( integer n,
integer k,
double precision, dimension( * ) d,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) work,
integer info )

DLAGSY

Purpose:
!>
!> DLAGSY generates a real symmetric matrix A, by pre- and post-
!> multiplying a real diagonal matrix D with a random orthogonal matrix:
!> A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
!> orthogonal transformations.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]K
!>          K is INTEGER
!>          The number of nonzero subdiagonals within the band of A.
!>          0 <= K <= N-1.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The diagonal elements of the diagonal matrix D.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          The generated n by n symmetric matrix A (the full matrix is
!>          stored).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= N.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry, the seed of the random number generator; the array
!>          elements must be between 0 and 4095, and ISEED(4) must be
!>          odd.
!>          On exit, the seed is updated.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 100 of file dlagsy.f.

101*
102* -- LAPACK auxiliary routine --
103* -- LAPACK is a software package provided by Univ. of Tennessee, --
104* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
105*
106* .. Scalar Arguments ..
107 INTEGER INFO, K, LDA, N
108* ..
109* .. Array Arguments ..
110 INTEGER ISEED( 4 )
111 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
112* ..
113*
114* =====================================================================
115*
116* .. Parameters ..
117 DOUBLE PRECISION ZERO, ONE, HALF
118 parameter( zero = 0.0d+0, one = 1.0d+0, half = 0.5d+0 )
119* ..
120* .. Local Scalars ..
121 INTEGER I, J
122 DOUBLE PRECISION ALPHA, TAU, WA, WB, WN
123* ..
124* .. External Subroutines ..
125 EXTERNAL daxpy, dgemv, dger, dlarnv, dscal, dsymv,
126 $ dsyr2, xerbla
127* ..
128* .. External Functions ..
129 DOUBLE PRECISION DDOT, DNRM2
130 EXTERNAL ddot, dnrm2
131* ..
132* .. Intrinsic Functions ..
133 INTRINSIC max, sign
134* ..
135* .. Executable Statements ..
136*
137* Test the input arguments
138*
139 info = 0
140 IF( n.LT.0 ) THEN
141 info = -1
142 ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
143 info = -2
144 ELSE IF( lda.LT.max( 1, n ) ) THEN
145 info = -5
146 END IF
147 IF( info.LT.0 ) THEN
148 CALL xerbla( 'DLAGSY', -info )
149 RETURN
150 END IF
151*
152* initialize lower triangle of A to diagonal matrix
153*
154 DO 20 j = 1, n
155 DO 10 i = j + 1, n
156 a( i, j ) = zero
157 10 CONTINUE
158 20 CONTINUE
159 DO 30 i = 1, n
160 a( i, i ) = d( i )
161 30 CONTINUE
162*
163* Generate lower triangle of symmetric matrix
164*
165 DO 40 i = n - 1, 1, -1
166*
167* generate random reflection
168*
169 CALL dlarnv( 3, iseed, n-i+1, work )
170 wn = dnrm2( n-i+1, work, 1 )
171 wa = sign( wn, work( 1 ) )
172 IF( wn.EQ.zero ) THEN
173 tau = zero
174 ELSE
175 wb = work( 1 ) + wa
176 CALL dscal( n-i, one / wb, work( 2 ), 1 )
177 work( 1 ) = one
178 tau = wb / wa
179 END IF
180*
181* apply random reflection to A(i:n,i:n) from the left
182* and the right
183*
184* compute y := tau * A * u
185*
186 CALL dsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
187 $ work( n+1 ), 1 )
188*
189* compute v := y - 1/2 * tau * ( y, u ) * u
190*
191 alpha = -half*tau*ddot( n-i+1, work( n+1 ), 1, work, 1 )
192 CALL daxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
193*
194* apply the transformation as a rank-2 update to A(i:n,i:n)
195*
196 CALL dsyr2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
197 $ a( i, i ), lda )
198 40 CONTINUE
199*
200* Reduce number of subdiagonals to K
201*
202 DO 60 i = 1, n - 1 - k
203*
204* generate reflection to annihilate A(k+i+1:n,i)
205*
206 wn = dnrm2( n-k-i+1, a( k+i, i ), 1 )
207 wa = sign( wn, a( k+i, i ) )
208 IF( wn.EQ.zero ) THEN
209 tau = zero
210 ELSE
211 wb = a( k+i, i ) + wa
212 CALL dscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
213 a( k+i, i ) = one
214 tau = wb / wa
215 END IF
216*
217* apply reflection to A(k+i:n,i+1:k+i-1) from the left
218*
219 CALL dgemv( 'Transpose', n-k-i+1, k-1, one, a( k+i, i+1 ), lda,
220 $ a( k+i, i ), 1, zero, work, 1 )
221 CALL dger( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
222 $ a( k+i, i+1 ), lda )
223*
224* apply reflection to A(k+i:n,k+i:n) from the left and the right
225*
226* compute y := tau * A * u
227*
228 CALL dsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
229 $ a( k+i, i ), 1, zero, work, 1 )
230*
231* compute v := y - 1/2 * tau * ( y, u ) * u
232*
233 alpha = -half*tau*ddot( n-k-i+1, work, 1, a( k+i, i ), 1 )
234 CALL daxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
235*
236* apply symmetric rank-2 update to A(k+i:n,k+i:n)
237*
238 CALL dsyr2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
239 $ a( k+i, k+i ), lda )
240*
241 a( k+i, i ) = -wa
242 DO 50 j = k + i + 1, n
243 a( j, i ) = zero
244 50 CONTINUE
245 60 CONTINUE
246*
247* Store full symmetric matrix
248*
249 DO 80 j = 1, n
250 DO 70 i = j + 1, n
251 a( j, i ) = a( i, j )
252 70 CONTINUE
253 80 CONTINUE
254 RETURN
255*
256* End of DLAGSY
257*
alpha
#define alpha
Definition eval.h:35
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dsymv
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
dsyr2
subroutine dsyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
DSYR2
Definition dsyr2.f:147

◆ dlahilb()

subroutine dlahilb ( integer n,
integer nrhs,
double precision, dimension(lda, n) a,
integer lda,
double precision, dimension(ldx, nrhs) x,
integer ldx,
double precision, dimension(ldb, nrhs) b,
integer ldb,
double precision, dimension(n) work,
integer info )

DLAHILB

Purpose:
!>
!> DLAHILB generates an N by N scaled Hilbert matrix in A along with
!> NRHS right-hand sides in B and solutions in X such that A*X=B.
!>
!> The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
!> entries are integers.  The right-hand sides are the first NRHS
!> columns of M * the identity matrix, and the solutions are the
!> first NRHS columns of the inverse Hilbert matrix.
!>
!> The condition number of the Hilbert matrix grows exponentially with
!> its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
!> Hilbert matrices beyond a relatively small dimension cannot be
!> generated exactly without extra precision.  Precision is exhausted
!> when the largest entry in the inverse Hilbert matrix is greater than
!> 2 to the power of the number of bits in the fraction of the data type
!> used plus one, which is 24 for single precision.
!>
!> In single, the generated solution is exact for N <= 6 and has
!> small componentwise error for 7 <= N <= 11.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The dimension of the matrix A.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The requested number of right-hand sides.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA, N)
!>          The generated scaled Hilbert matrix.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= N.
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (LDX, NRHS)
!>          The generated exact solutions.  Currently, the first NRHS
!>          columns of the inverse Hilbert matrix.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= N.
!>
[out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, NRHS)
!>          The generated right-hand sides.  Currently, the first NRHS
!>          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= N.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          = 1: N is too large; the data is still generated but may not
!>               be not exact.
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 123 of file dlahilb.f.

124*
125* -- LAPACK test routine --
126* -- LAPACK is a software package provided by Univ. of Tennessee, --
127* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
128*
129* .. Scalar Arguments ..
130 INTEGER N, NRHS, LDA, LDX, LDB, INFO
131* .. Array Arguments ..
132 DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
133* ..
134*
135* =====================================================================
136* .. Local Scalars ..
137 INTEGER TM, TI, R
138 INTEGER M
139 INTEGER I, J
140* ..
141* .. Parameters ..
142* NMAX_EXACT the largest dimension where the generated data is
143* exact.
144* NMAX_APPROX the largest dimension where the generated data has
145* a small componentwise relative error.
146 INTEGER NMAX_EXACT, NMAX_APPROX
147 parameter(nmax_exact = 6, nmax_approx = 11)
148* ..
149* .. External Subroutines ..
150 EXTERNAL xerbla
151* ..
152* .. External Functions
153 EXTERNAL dlaset
154 INTRINSIC dble
155* ..
156* .. Executable Statements ..
157*
158* Test the input arguments
159*
160 info = 0
161 IF (n .LT. 0 .OR. n .GT. nmax_approx) THEN
162 info = -1
163 ELSE IF (nrhs .LT. 0) THEN
164 info = -2
165 ELSE IF (lda .LT. n) THEN
166 info = -4
167 ELSE IF (ldx .LT. n) THEN
168 info = -6
169 ELSE IF (ldb .LT. n) THEN
170 info = -8
171 END IF
172 IF (info .LT. 0) THEN
173 CALL xerbla('DLAHILB', -info)
174 RETURN
175 END IF
176 IF (n .GT. nmax_exact) THEN
177 info = 1
178 END IF
179*
180* Compute M = the LCM of the integers [1, 2*N-1]. The largest
181* reasonable N is small enough that integers suffice (up to N = 11).
182 m = 1
183 DO i = 2, (2*n-1)
184 tm = m
185 ti = i
186 r = mod(tm, ti)
187 DO WHILE (r .NE. 0)
188 tm = ti
189 ti = r
190 r = mod(tm, ti)
191 END DO
192 m = (m / ti) * i
193 END DO
194*
195* Generate the scaled Hilbert matrix in A
196 DO j = 1, n
197 DO i = 1, n
198 a(i, j) = dble(m) / (i + j - 1)
199 END DO
200 END DO
201*
202* Generate matrix B as simply the first NRHS columns of M * the
203* identity.
204 CALL dlaset('Full', n, nrhs, 0.0d+0, dble(m), b, ldb)
205
206* Generate the true solutions in X. Because B = the first NRHS
207* columns of M*I, the true solutions are just the first NRHS columns
208* of the inverse Hilbert matrix.
209 work(1) = n
210 DO j = 2, n
211 work(j) = ( ( (work(j-1)/(j-1)) * (j-1 - n) ) /(j-1) )
212 $ * (n +j -1)
213 END DO
214*
215 DO j = 1, nrhs
216 DO i = 1, n
217 x(i, j) = (work(i)*work(j)) / (i + j - 1)
218 END DO
219 END DO
220*
dlaset
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.
Definition dlaset.f:110

◆ dlakf2()

subroutine dlakf2 ( integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( lda, * ) b,
double precision, dimension( lda, * ) d,
double precision, dimension( lda, * ) e,
double precision, dimension( ldz, * ) z,
integer ldz )

DLAKF2

Purpose:
!>
!> Form the 2*M*N by 2*M*N matrix
!>
!>        Z = [ kron(In, A)  -kron(B', Im) ]
!>            [ kron(In, D)  -kron(E', Im) ],
!>
!> where In is the identity matrix of size n and X' is the transpose
!> of X. kron(X, Y) is the Kronecker product between the matrices X
!> and Y.
!>
Parameters
[in]M
!>          M is INTEGER
!>          Size of matrix, must be >= 1.
!>
[in]N
!>          N is INTEGER
!>          Size of matrix, must be >= 1.
!>
[in]A
!>          A is DOUBLE PRECISION, dimension ( LDA, M )
!>          The matrix A in the output matrix Z.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A, B, D, and E. ( LDA >= M+N )
!>
[in]B
!>          B is DOUBLE PRECISION, dimension ( LDA, N )
!>
[in]D
!>          D is DOUBLE PRECISION, dimension ( LDA, M )
!>
[in]E
!>          E is DOUBLE PRECISION, dimension ( LDA, N )
!>
!>          The matrices used in forming the output matrix Z.
!>
[out]Z
!>          Z is DOUBLE PRECISION, dimension ( LDZ, 2*M*N )
!>          The resultant Kronecker M*N*2 by M*N*2 matrix (see above.)
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of Z. ( LDZ >= 2*M*N )
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 104 of file dlakf2.f.

105*
106* -- LAPACK computational routine --
107* -- LAPACK is a software package provided by Univ. of Tennessee, --
108* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
109*
110* .. Scalar Arguments ..
111 INTEGER LDA, LDZ, M, N
112* ..
113* .. Array Arguments ..
114 DOUBLE PRECISION A( LDA, * ), B( LDA, * ), D( LDA, * ),
115 $ E( LDA, * ), Z( LDZ, * )
116* ..
117*
118* ====================================================================
119*
120* .. Parameters ..
121 DOUBLE PRECISION ZERO
122 parameter( zero = 0.0d+0 )
123* ..
124* .. Local Scalars ..
125 INTEGER I, IK, J, JK, L, MN, MN2
126* ..
127* .. External Subroutines ..
128 EXTERNAL dlaset
129* ..
130* .. Executable Statements ..
131*
132* Initialize Z
133*
134 mn = m*n
135 mn2 = 2*mn
136 CALL dlaset( 'Full', mn2, mn2, zero, zero, z, ldz )
137*
138 ik = 1
139 DO 50 l = 1, n
140*
141* form kron(In, A)
142*
143 DO 20 i = 1, m
144 DO 10 j = 1, m
145 z( ik+i-1, ik+j-1 ) = a( i, j )
146 10 CONTINUE
147 20 CONTINUE
148*
149* form kron(In, D)
150*
151 DO 40 i = 1, m
152 DO 30 j = 1, m
153 z( ik+mn+i-1, ik+j-1 ) = d( i, j )
154 30 CONTINUE
155 40 CONTINUE
156*
157 ik = ik + m
158 50 CONTINUE
159*
160 ik = 1
161 DO 90 l = 1, n
162 jk = mn + 1
163*
164 DO 80 j = 1, n
165*
166* form -kron(B', Im)
167*
168 DO 60 i = 1, m
169 z( ik+i-1, jk+i-1 ) = -b( j, l )
170 60 CONTINUE
171*
172* form -kron(E', Im)
173*
174 DO 70 i = 1, m
175 z( ik+mn+i-1, jk+i-1 ) = -e( j, l )
176 70 CONTINUE
177*
178 jk = jk + m
179 80 CONTINUE
180*
181 ik = ik + m
182 90 CONTINUE
183*
184 RETURN
185*
186* End of DLAKF2
187*

◆ dlarge()

subroutine dlarge ( integer n,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) work,
integer info )

DLARGE

Purpose:
!>
!> DLARGE pre- and post-multiplies a real general n by n matrix A
!> with a random orthogonal matrix: A = U*D*U'.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the original n by n matrix A.
!>          On exit, A is overwritten by U*A*U' for some random
!>          orthogonal matrix U.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= N.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry, the seed of the random number generator; the array
!>          elements must be between 0 and 4095, and ISEED(4) must be
!>          odd.
!>          On exit, the seed is updated.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 86 of file dlarge.f.

87*
88* -- LAPACK auxiliary routine --
89* -- LAPACK is a software package provided by Univ. of Tennessee, --
90* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
91*
92* .. Scalar Arguments ..
93 INTEGER INFO, LDA, N
94* ..
95* .. Array Arguments ..
96 INTEGER ISEED( 4 )
97 DOUBLE PRECISION A( LDA, * ), WORK( * )
98* ..
99*
100* =====================================================================
101*
102* .. Parameters ..
103 DOUBLE PRECISION ZERO, ONE
104 parameter( zero = 0.0d+0, one = 1.0d+0 )
105* ..
106* .. Local Scalars ..
107 INTEGER I
108 DOUBLE PRECISION TAU, WA, WB, WN
109* ..
110* .. External Subroutines ..
111 EXTERNAL dgemv, dger, dlarnv, dscal, xerbla
112* ..
113* .. Intrinsic Functions ..
114 INTRINSIC max, sign
115* ..
116* .. External Functions ..
117 DOUBLE PRECISION DNRM2
118 EXTERNAL dnrm2
119* ..
120* .. Executable Statements ..
121*
122* Test the input arguments
123*
124 info = 0
125 IF( n.LT.0 ) THEN
126 info = -1
127 ELSE IF( lda.LT.max( 1, n ) ) THEN
128 info = -3
129 END IF
130 IF( info.LT.0 ) THEN
131 CALL xerbla( 'DLARGE', -info )
132 RETURN
133 END IF
134*
135* pre- and post-multiply A by random orthogonal matrix
136*
137 DO 10 i = n, 1, -1
138*
139* generate random reflection
140*
141 CALL dlarnv( 3, iseed, n-i+1, work )
142 wn = dnrm2( n-i+1, work, 1 )
143 wa = sign( wn, work( 1 ) )
144 IF( wn.EQ.zero ) THEN
145 tau = zero
146 ELSE
147 wb = work( 1 ) + wa
148 CALL dscal( n-i, one / wb, work( 2 ), 1 )
149 work( 1 ) = one
150 tau = wb / wa
151 END IF
152*
153* multiply A(i:n,1:n) by random reflection from the left
154*
155 CALL dgemv( 'Transpose', n-i+1, n, one, a( i, 1 ), lda, work,
156 $ 1, zero, work( n+1 ), 1 )
157 CALL dger( n-i+1, n, -tau, work, 1, work( n+1 ), 1, a( i, 1 ),
158 $ lda )
159*
160* multiply A(1:n,i:n) by random reflection from the right
161*
162 CALL dgemv( 'No transpose', n, n-i+1, one, a( 1, i ), lda,
163 $ work, 1, zero, work( n+1 ), 1 )
164 CALL dger( n, n-i+1, -tau, work( n+1 ), 1, work, 1, a( 1, i ),
165 $ lda )
166 10 CONTINUE
167 RETURN
168*
169* End of DLARGE
170*

◆ dlarnd()

double precision function dlarnd ( integer idist,
integer, dimension( 4 ) iseed )

DLARND

Purpose:
!>
!> DLARND returns a random real number from a uniform or normal
!> distribution.
!>
Parameters
[in]IDIST
!>          IDIST is INTEGER
!>          Specifies the distribution of the random numbers:
!>          = 1:  uniform (0,1)
!>          = 2:  uniform (-1,1)
!>          = 3:  normal (0,1)
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry, the seed of the random number generator; the array
!>          elements must be between 0 and 4095, and ISEED(4) must be
!>          odd.
!>          On exit, the seed is updated.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  This routine calls the auxiliary routine DLARAN to generate a random
!>  real number from a uniform (0,1) distribution. The Box-Muller method
!>  is used to transform numbers from a uniform to a normal distribution.
!>

Definition at line 72 of file dlarnd.f.

73*
74* -- LAPACK auxiliary routine --
75* -- LAPACK is a software package provided by Univ. of Tennessee, --
76* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
77*
78* .. Scalar Arguments ..
79 INTEGER IDIST
80* ..
81* .. Array Arguments ..
82 INTEGER ISEED( 4 )
83* ..
84*
85* =====================================================================
86*
87* .. Parameters ..
88 DOUBLE PRECISION ONE, TWO
89 parameter( one = 1.0d+0, two = 2.0d+0 )
90 DOUBLE PRECISION TWOPI
91 parameter( twopi = 6.28318530717958647692528676655900576839d+0 )
92* ..
93* .. Local Scalars ..
94 DOUBLE PRECISION T1, T2
95* ..
96* .. External Functions ..
97 DOUBLE PRECISION DLARAN
98 EXTERNAL dlaran
99* ..
100* .. Intrinsic Functions ..
101 INTRINSIC cos, log, sqrt
102* ..
103* .. Executable Statements ..
104*
105* Generate a real random number from a uniform (0,1) distribution
106*
107 t1 = dlaran( iseed )
108*
109 IF( idist.EQ.1 ) THEN
110*
111* uniform (0,1)
112*
113 dlarnd = t1
114 ELSE IF( idist.EQ.2 ) THEN
115*
116* uniform (-1,1)
117*
118 dlarnd = two*t1 - one
119 ELSE IF( idist.EQ.3 ) THEN
120*
121* normal (0,1)
122*
123 t2 = dlaran( iseed )
124 dlarnd = sqrt( -two*log( t1 ) )*cos( twopi*t2 )
125 END IF
126 RETURN
127*
128* End of DLARND
129*
dlarnd
double precision function dlarnd(idist, iseed)
DLARND
Definition dlarnd.f:73
dlaran
double precision function dlaran(iseed)
DLARAN
Definition dlaran.f:67

◆ dlaror()

subroutine dlaror ( character side,
character init,
integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) x,
integer info )

DLAROR

Purpose:
!>
!> DLAROR pre- or post-multiplies an M by N matrix A by a random
!> orthogonal matrix U, overwriting A.  A may optionally be initialized
!> to the identity matrix before multiplying by U.  U is generated using
!> the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
!>
Parameters
[in]SIDE
!>          SIDE is CHARACTER*1
!>          Specifies whether A is multiplied on the left or right by U.
!>          = 'L':         Multiply A on the left (premultiply) by U
!>          = 'R':         Multiply A on the right (postmultiply) by U'
!>          = 'C' or 'T':  Multiply A on the left by U and the right
!>                          by U' (Here, U' means U-transpose.)
!>
[in]INIT
!>          INIT is CHARACTER*1
!>          Specifies whether or not A should be initialized to the
!>          identity matrix.
!>          = 'I':  Initialize A to (a section of) the identity matrix
!>                   before applying U.
!>          = 'N':  No initialization.  Apply U to the input matrix A.
!>
!>          INIT = 'I' may be used to generate square or rectangular
!>          orthogonal matrices:
!>
!>          For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
!>          to each other, as will the columns.
!>
!>          If M < N, SIDE = 'R' produces a dense matrix whose rows are
!>          orthogonal and whose columns are not, while SIDE = 'L'
!>          produces a matrix whose rows are orthogonal, and whose first
!>          M columns are orthogonal, and whose remaining columns are
!>          zero.
!>
!>          If M > N, SIDE = 'L' produces a dense matrix whose columns
!>          are orthogonal and whose rows are not, while SIDE = 'R'
!>          produces a matrix whose columns are orthogonal, and whose
!>          first M rows are orthogonal, and whose remaining rows are
!>          zero.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of A.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of A.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA, N)
!>          On entry, the array A.
!>          On exit, overwritten by U A ( if SIDE = 'L' ),
!>           or by A U ( if SIDE = 'R' ),
!>           or by U A U' ( if SIDE = 'C' or 'T').
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry ISEED specifies the seed of the random number
!>          generator. The array elements should be between 0 and 4095;
!>          if not they will be reduced mod 4096.  Also, ISEED(4) must
!>          be odd.  The random number generator uses a linear
!>          congruential sequence limited to small integers, and so
!>          should produce machine independent random numbers. The
!>          values of ISEED are changed on exit, and can be used in the
!>          next call to DLAROR to continue the same random number
!>          sequence.
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (3*MAX( M, N ))
!>          Workspace of length
!>              2*M + N if SIDE = 'L',
!>              2*N + M if SIDE = 'R',
!>              3*N     if SIDE = 'C' or 'T'.
!>
[out]INFO
!>          INFO is INTEGER
!>          An error flag.  It is set to:
!>          = 0:  normal return
!>          < 0:  if INFO = -k, the k-th argument had an illegal value
!>          = 1:  if the random numbers generated by DLARND are bad.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 145 of file dlaror.f.

146*
147* -- LAPACK auxiliary routine --
148* -- LAPACK is a software package provided by Univ. of Tennessee, --
149* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150*
151* .. Scalar Arguments ..
152 CHARACTER INIT, SIDE
153 INTEGER INFO, LDA, M, N
154* ..
155* .. Array Arguments ..
156 INTEGER ISEED( 4 )
157 DOUBLE PRECISION A( LDA, * ), X( * )
158* ..
159*
160* =====================================================================
161*
162* .. Parameters ..
163 DOUBLE PRECISION ZERO, ONE, TOOSML
164 parameter( zero = 0.0d+0, one = 1.0d+0,
165 $ toosml = 1.0d-20 )
166* ..
167* .. Local Scalars ..
168 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
169 DOUBLE PRECISION FACTOR, XNORM, XNORMS
170* ..
171* .. External Functions ..
172 LOGICAL LSAME
173 DOUBLE PRECISION DLARND, DNRM2
174 EXTERNAL lsame, dlarnd, dnrm2
175* ..
176* .. External Subroutines ..
177 EXTERNAL dgemv, dger, dlaset, dscal, xerbla
178* ..
179* .. Intrinsic Functions ..
180 INTRINSIC abs, sign
181* ..
182* .. Executable Statements ..
183*
184 info = 0
185 IF( n.EQ.0 .OR. m.EQ.0 )
186 $ RETURN
187*
188 itype = 0
189 IF( lsame( side, 'L' ) ) THEN
190 itype = 1
191 ELSE IF( lsame( side, 'R' ) ) THEN
192 itype = 2
193 ELSE IF( lsame( side, 'C' ) .OR. lsame( side, 'T' ) ) THEN
194 itype = 3
195 END IF
196*
197* Check for argument errors.
198*
199 IF( itype.EQ.0 ) THEN
200 info = -1
201 ELSE IF( m.LT.0 ) THEN
202 info = -3
203 ELSE IF( n.LT.0 .OR. ( itype.EQ.3 .AND. n.NE.m ) ) THEN
204 info = -4
205 ELSE IF( lda.LT.m ) THEN
206 info = -6
207 END IF
208 IF( info.NE.0 ) THEN
209 CALL xerbla( 'DLAROR', -info )
210 RETURN
211 END IF
212*
213 IF( itype.EQ.1 ) THEN
214 nxfrm = m
215 ELSE
216 nxfrm = n
217 END IF
218*
219* Initialize A to the identity matrix if desired
220*
221 IF( lsame( init, 'I' ) )
222 $ CALL dlaset( 'Full', m, n, zero, one, a, lda )
223*
224* If no rotation possible, multiply by random +/-1
225*
226* Compute rotation by computing Householder transformations
227* H(2), H(3), ..., H(nhouse)
228*
229 DO 10 j = 1, nxfrm
230 x( j ) = zero
231 10 CONTINUE
232*
233 DO 30 ixfrm = 2, nxfrm
234 kbeg = nxfrm - ixfrm + 1
235*
236* Generate independent normal( 0, 1 ) random numbers
237*
238 DO 20 j = kbeg, nxfrm
239 x( j ) = dlarnd( 3, iseed )
240 20 CONTINUE
241*
242* Generate a Householder transformation from the random vector X
243*
244 xnorm = dnrm2( ixfrm, x( kbeg ), 1 )
245 xnorms = sign( xnorm, x( kbeg ) )
246 x( kbeg+nxfrm ) = sign( one, -x( kbeg ) )
247 factor = xnorms*( xnorms+x( kbeg ) )
248 IF( abs( factor ).LT.toosml ) THEN
249 info = 1
250 CALL xerbla( 'DLAROR', info )
251 RETURN
252 ELSE
253 factor = one / factor
254 END IF
255 x( kbeg ) = x( kbeg ) + xnorms
256*
257* Apply Householder transformation to A
258*
259 IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
260*
261* Apply H(k) from the left.
262*
263 CALL dgemv( 'T', ixfrm, n, one, a( kbeg, 1 ), lda,
264 $ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
265 CALL dger( ixfrm, n, -factor, x( kbeg ), 1, x( 2*nxfrm+1 ),
266 $ 1, a( kbeg, 1 ), lda )
267*
268 END IF
269*
270 IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
271*
272* Apply H(k) from the right.
273*
274 CALL dgemv( 'N', m, ixfrm, one, a( 1, kbeg ), lda,
275 $ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
276 CALL dger( m, ixfrm, -factor, x( 2*nxfrm+1 ), 1, x( kbeg ),
277 $ 1, a( 1, kbeg ), lda )
278*
279 END IF
280 30 CONTINUE
281*
282 x( 2*nxfrm ) = sign( one, dlarnd( 3, iseed ) )
283*
284* Scale the matrix A by D.
285*
286 IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
287 DO 40 irow = 1, m
288 CALL dscal( n, x( nxfrm+irow ), a( irow, 1 ), lda )
289 40 CONTINUE
290 END IF
291*
292 IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
293 DO 50 jcol = 1, n
294 CALL dscal( m, x( nxfrm+jcol ), a( 1, jcol ), 1 )
295 50 CONTINUE
296 END IF
297 RETURN
298*
299* End of DLAROR
300*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53

◆ dlarot()

subroutine dlarot ( logical lrows,
logical lleft,
logical lright,
integer nl,
double precision c,
double precision s,
double precision, dimension( * ) a,
integer lda,
double precision xleft,
double precision xright )

DLAROT

Purpose:
!>
!>    DLAROT applies a (Givens) rotation to two adjacent rows or
!>    columns, where one element of the first and/or last column/row
!>    for use on matrices stored in some format other than GE, so
!>    that elements of the matrix may be used or modified for which
!>    no array element is provided.
!>
!>    One example is a symmetric matrix in SB format (bandwidth=4), for
!>    which UPLO='L':  Two adjacent rows will have the format:
!>
!>    row j:     C> C> C> C> C> .  .  .  .
!>    row j+1:      C> C> C> C> C> .  .  .  .
!>
!>    '*' indicates elements for which storage is provided,
!>    '.' indicates elements for which no storage is provided, but
!>    are not necessarily zero; their values are determined by
!>    symmetry.  ' ' indicates elements which are necessarily zero,
!>     and have no storage provided.
!>
!>    Those columns which have two '*'s can be handled by DROT.
!>    Those columns which have no '*'s can be ignored, since as long
!>    as the Givens rotations are carefully applied to preserve
!>    symmetry, their values are determined.
!>    Those columns which have one '*' have to be handled separately,
!>    by using separate variables  and :
!>
!>    row j:     C> C> C> C> C> p  .  .  .
!>    row j+1:   q  C> C> C> C> C> .  .  .  .
!>
!>    The element p would have to be set correctly, then that column
!>    is rotated, setting p to its new value.  The next call to
!>    DLAROT would rotate columns j and j+1, using p, and restore
!>    symmetry.  The element q would start out being zero, and be
!>    made non-zero by the rotation.  Later, rotations would presumably
!>    be chosen to zero q out.
!>
!>    Typical Calling Sequences: rotating the i-th and (i+1)-st rows.
!>    ------- ------- ---------
!>
!>      General dense matrix:
!>
!>              CALL DLAROT(.TRUE.,.FALSE.,.FALSE., N, C,S,
!>                      A(i,1),LDA, DUMMY, DUMMY)
!>
!>      General banded matrix in GB format:
!>
!>              j = MAX(1, i-KL )
!>              NL = MIN( N, i+KU+1 ) + 1-j
!>              CALL DLAROT( .TRUE., i-KL.GE.1, i+KU.LT.N, NL, C,S,
!>                      A(KU+i+1-j,j),LDA-1, XLEFT, XRIGHT )
!>
!>              [ note that i+1-j is just MIN(i,KL+1) ]
!>
!>      Symmetric banded matrix in SY format, bandwidth K,
!>      lower triangle only:
!>
!>              j = MAX(1, i-K )
!>              NL = MIN( K+1, i ) + 1
!>              CALL DLAROT( .TRUE., i-K.GE.1, .TRUE., NL, C,S,
!>                      A(i,j), LDA, XLEFT, XRIGHT )
!>
!>      Same, but upper triangle only:
!>
!>              NL = MIN( K+1, N-i ) + 1
!>              CALL DLAROT( .TRUE., .TRUE., i+K.LT.N, NL, C,S,
!>                      A(i,i), LDA, XLEFT, XRIGHT )
!>
!>      Symmetric banded matrix in SB format, bandwidth K,
!>      lower triangle only:
!>
!>              [ same as for SY, except:]
!>                  . . . .
!>                      A(i+1-j,j), LDA-1, XLEFT, XRIGHT )
!>
!>              [ note that i+1-j is just MIN(i,K+1) ]
!>
!>      Same, but upper triangle only:
!>                   . . .
!>                      A(K+1,i), LDA-1, XLEFT, XRIGHT )
!>
!>      Rotating columns is just the transpose of rotating rows, except
!>      for GB and SB: (rotating columns i and i+1)
!>
!>      GB:
!>              j = MAX(1, i-KU )
!>              NL = MIN( N, i+KL+1 ) + 1-j
!>              CALL DLAROT( .TRUE., i-KU.GE.1, i+KL.LT.N, NL, C,S,
!>                      A(KU+j+1-i,i),LDA-1, XTOP, XBOTTM )
!>
!>              [note that KU+j+1-i is just MAX(1,KU+2-i)]
!>
!>      SB: (upper triangle)
!>
!>                   . . . . . .
!>                      A(K+j+1-i,i),LDA-1, XTOP, XBOTTM )
!>
!>      SB: (lower triangle)
!>
!>                   . . . . . .
!>                      A(1,i),LDA-1, XTOP, XBOTTM )
!> 
!>  LROWS  - LOGICAL
!>           If .TRUE., then DLAROT will rotate two rows.  If .FALSE.,
!>           then it will rotate two columns.
!>           Not modified.
!>
!>  LLEFT  - LOGICAL
!>           If .TRUE., then XLEFT will be used instead of the
!>           corresponding element of A for the first element in the
!>           second row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.)
!>           If .FALSE., then the corresponding element of A will be
!>           used.
!>           Not modified.
!>
!>  LRIGHT - LOGICAL
!>           If .TRUE., then XRIGHT will be used instead of the
!>           corresponding element of A for the last element in the
!>           first row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) If
!>           .FALSE., then the corresponding element of A will be used.
!>           Not modified.
!>
!>  NL     - INTEGER
!>           The length of the rows (if LROWS=.TRUE.) or columns (if
!>           LROWS=.FALSE.) to be rotated.  If XLEFT and/or XRIGHT are
!>           used, the columns/rows they are in should be included in
!>           NL, e.g., if LLEFT = LRIGHT = .TRUE., then NL must be at
!>           least 2.  The number of rows/columns to be rotated
!>           exclusive of those involving XLEFT and/or XRIGHT may
!>           not be negative, i.e., NL minus how many of LLEFT and
!>           LRIGHT are .TRUE. must be at least zero; if not, XERBLA
!>           will be called.
!>           Not modified.
!>
!>  C, S   - DOUBLE PRECISION
!>           Specify the Givens rotation to be applied.  If LROWS is
!>           true, then the matrix ( c  s )
!>                                 (-s  c )  is applied from the left;
!>           if false, then the transpose thereof is applied from the
!>           right.  For a Givens rotation, C**2 + S**2 should be 1,
!>           but this is not checked.
!>           Not modified.
!>
!>  A      - DOUBLE PRECISION array.
!>           The array containing the rows/columns to be rotated.  The
!>           first element of A should be the upper left element to
!>           be rotated.
!>           Read and modified.
!>
!>  LDA    - INTEGER
!>           The  leading dimension of A.  If A contains
!>           a matrix stored in GE or SY format, then this is just
!>           the leading dimension of A as dimensioned in the calling
!>           routine.  If A contains a matrix stored in band (GB or SB)
!>           format, then this should be *one less* than the leading
!>           dimension used in the calling routine.  Thus, if
!>           A were dimensioned A(LDA,*) in DLAROT, then A(1,j) would
!>           be the j-th element in the first of the two rows
!>           to be rotated, and A(2,j) would be the j-th in the second,
!>           regardless of how the array may be stored in the calling
!>           routine.  [A cannot, however, actually be dimensioned thus,
!>           since for band format, the row number may exceed LDA, which
!>           is not legal FORTRAN.]
!>           If LROWS=.TRUE., then LDA must be at least 1, otherwise
!>           it must be at least NL minus the number of .TRUE. values
!>           in XLEFT and XRIGHT.
!>           Not modified.
!>
!>  XLEFT  - DOUBLE PRECISION
!>           If LLEFT is .TRUE., then XLEFT will be used and modified
!>           instead of A(2,1) (if LROWS=.TRUE.) or A(1,2)
!>           (if LROWS=.FALSE.).
!>           Read and modified.
!>
!>  XRIGHT - DOUBLE PRECISION
!>           If LRIGHT is .TRUE., then XRIGHT will be used and modified
!>           instead of A(1,NL) (if LROWS=.TRUE.) or A(NL,1)
!>           (if LROWS=.FALSE.).
!>           Read and modified.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file dlarot.f.

226*
227* -- LAPACK auxiliary routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 LOGICAL LLEFT, LRIGHT, LROWS
233 INTEGER LDA, NL
234 DOUBLE PRECISION C, S, XLEFT, XRIGHT
235* ..
236* .. Array Arguments ..
237 DOUBLE PRECISION A( * )
238* ..
239*
240* =====================================================================
241*
242* .. Local Scalars ..
243 INTEGER IINC, INEXT, IX, IY, IYT, NT
244* ..
245* .. Local Arrays ..
246 DOUBLE PRECISION XT( 2 ), YT( 2 )
247* ..
248* .. External Subroutines ..
249 EXTERNAL drot, xerbla
250* ..
251* .. Executable Statements ..
252*
253* Set up indices, arrays for ends
254*
255 IF( lrows ) THEN
256 iinc = lda
257 inext = 1
258 ELSE
259 iinc = 1
260 inext = lda
261 END IF
262*
263 IF( lleft ) THEN
264 nt = 1
265 ix = 1 + iinc
266 iy = 2 + lda
267 xt( 1 ) = a( 1 )
268 yt( 1 ) = xleft
269 ELSE
270 nt = 0
271 ix = 1
272 iy = 1 + inext
273 END IF
274*
275 IF( lright ) THEN
276 iyt = 1 + inext + ( nl-1 )*iinc
277 nt = nt + 1
278 xt( nt ) = xright
279 yt( nt ) = a( iyt )
280 END IF
281*
282* Check for errors
283*
284 IF( nl.LT.nt ) THEN
285 CALL xerbla( 'DLAROT', 4 )
286 RETURN
287 END IF
288 IF( lda.LE.0 .OR. ( .NOT.lrows .AND. lda.LT.nl-nt ) ) THEN
289 CALL xerbla( 'DLAROT', 8 )
290 RETURN
291 END IF
292*
293* Rotate
294*
295 CALL drot( nl-nt, a( ix ), iinc, a( iy ), iinc, c, s )
296 CALL drot( nt, xt, 1, yt, 1, c, s )
297*
298* Stuff values back into XLEFT, XRIGHT, etc.
299*
300 IF( lleft ) THEN
301 a( 1 ) = xt( 1 )
302 xleft = yt( 1 )
303 END IF
304*
305 IF( lright ) THEN
306 xright = xt( nt )
307 a( iyt ) = yt( nt )
308 END IF
309*
310 RETURN
311*
312* End of DLAROT
313*
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
nl
character *2 function nl()
Definition message.F:2354

◆ dlatm1()

subroutine dlatm1 ( integer mode,
double precision cond,
integer irsign,
integer idist,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) d,
integer n,
integer info )

DLATM1

Purpose:
!>
!>    DLATM1 computes the entries of D(1..N) as specified by
!>    MODE, COND and IRSIGN. IDIST and ISEED determine the generation
!>    of random numbers. DLATM1 is called by DLATMR to generate
!>    random test matrices for LAPACK programs.
!>
Parameters
[in]MODE
!>          MODE is INTEGER
!>           On entry describes how D is to be computed:
!>           MODE = 0 means do not change D.
!>           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
!>           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           Not modified.
!>
[in]COND
!>          COND is DOUBLE PRECISION
!>           On entry, used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
[in]IRSIGN
!>          IRSIGN is INTEGER
!>           On entry, if MODE neither -6, 0 nor 6, determines sign of
!>           entries of D
!>           0 => leave entries of D unchanged
!>           1 => multiply each entry of D by 1 or -1 with probability .5
!>
[in]IDIST
!>          IDIST is INTEGER
!>           On entry, IDIST specifies the type of distribution to be
!>           used to generate a random matrix .
!>           1 => UNIFORM( 0, 1 )
!>           2 => UNIFORM( -1, 1 )
!>           3 => NORMAL( 0, 1 )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. The random number generator uses a
!>           linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATM1
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension ( N )
!>           Array to be computed according to MODE, COND and IRSIGN.
!>           May be changed on exit if MODE is nonzero.
!>
[in]N
!>          N is INTEGER
!>           Number of entries of D. Not modified.
!>
[out]INFO
!>          INFO is INTEGER
!>            0  => normal termination
!>           -1  => if MODE not in range -6 to 6
!>           -2  => if MODE neither -6, 0 nor 6, and
!>                  IRSIGN neither 0 nor 1
!>           -3  => if MODE neither -6, 0 nor 6 and COND less than 1
!>           -4  => if MODE equals 6 or -6 and IDIST not in range 1 to 3
!>           -7  => if N negative
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 134 of file dlatm1.f.

135*
136* -- LAPACK auxiliary routine --
137* -- LAPACK is a software package provided by Univ. of Tennessee, --
138* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
139*
140* .. Scalar Arguments ..
141 INTEGER IDIST, INFO, IRSIGN, MODE, N
142 DOUBLE PRECISION COND
143* ..
144* .. Array Arguments ..
145 INTEGER ISEED( 4 )
146 DOUBLE PRECISION D( * )
147* ..
148*
149* =====================================================================
150*
151* .. Parameters ..
152 DOUBLE PRECISION ONE
153 parameter( one = 1.0d0 )
154 DOUBLE PRECISION HALF
155 parameter( half = 0.5d0 )
156* ..
157* .. Local Scalars ..
158 INTEGER I
159 DOUBLE PRECISION ALPHA, TEMP
160* ..
161* .. External Functions ..
162 DOUBLE PRECISION DLARAN
163 EXTERNAL dlaran
164* ..
165* .. External Subroutines ..
166 EXTERNAL dlarnv, xerbla
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC abs, dble, exp, log
170* ..
171* .. Executable Statements ..
172*
173* Decode and Test the input parameters. Initialize flags & seed.
174*
175 info = 0
176*
177* Quick return if possible
178*
179 IF( n.EQ.0 )
180 $ RETURN
181*
182* Set INFO if an error
183*
184 IF( mode.LT.-6 .OR. mode.GT.6 ) THEN
185 info = -1
186 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
187 $ ( irsign.NE.0 .AND. irsign.NE.1 ) ) THEN
188 info = -2
189 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
190 $ cond.LT.one ) THEN
191 info = -3
192 ELSE IF( ( mode.EQ.6 .OR. mode.EQ.-6 ) .AND.
193 $ ( idist.LT.1 .OR. idist.GT.3 ) ) THEN
194 info = -4
195 ELSE IF( n.LT.0 ) THEN
196 info = -7
197 END IF
198*
199 IF( info.NE.0 ) THEN
200 CALL xerbla( 'DLATM1', -info )
201 RETURN
202 END IF
203*
204* Compute D according to COND and MODE
205*
206 IF( mode.NE.0 ) THEN
207 GO TO ( 10, 30, 50, 70, 90, 110 )abs( mode )
208*
209* One large D value:
210*
211 10 CONTINUE
212 DO 20 i = 1, n
213 d( i ) = one / cond
214 20 CONTINUE
215 d( 1 ) = one
216 GO TO 120
217*
218* One small D value:
219*
220 30 CONTINUE
221 DO 40 i = 1, n
222 d( i ) = one
223 40 CONTINUE
224 d( n ) = one / cond
225 GO TO 120
226*
227* Exponentially distributed D values:
228*
229 50 CONTINUE
230 d( 1 ) = one
231 IF( n.GT.1 ) THEN
232 alpha = cond**( -one / dble( n-1 ) )
233 DO 60 i = 2, n
234 d( i ) = alpha**( i-1 )
235 60 CONTINUE
236 END IF
237 GO TO 120
238*
239* Arithmetically distributed D values:
240*
241 70 CONTINUE
242 d( 1 ) = one
243 IF( n.GT.1 ) THEN
244 temp = one / cond
245 alpha = ( one-temp ) / dble( n-1 )
246 DO 80 i = 2, n
247 d( i ) = dble( n-i )*alpha + temp
248 80 CONTINUE
249 END IF
250 GO TO 120
251*
252* Randomly distributed D values on ( 1/COND , 1):
253*
254 90 CONTINUE
255 alpha = log( one / cond )
256 DO 100 i = 1, n
257 d( i ) = exp( alpha*dlaran( iseed ) )
258 100 CONTINUE
259 GO TO 120
260*
261* Randomly distributed D values from IDIST
262*
263 110 CONTINUE
264 CALL dlarnv( idist, iseed, n, d )
265*
266 120 CONTINUE
267*
268* If MODE neither -6 nor 0 nor 6, and IRSIGN = 1, assign
269* random signs to D
270*
271 IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
272 $ irsign.EQ.1 ) THEN
273 DO 130 i = 1, n
274 temp = dlaran( iseed )
275 IF( temp.GT.half )
276 $ d( i ) = -d( i )
277 130 CONTINUE
278 END IF
279*
280* Reverse if MODE < 0
281*
282 IF( mode.LT.0 ) THEN
283 DO 140 i = 1, n / 2
284 temp = d( i )
285 d( i ) = d( n+1-i )
286 d( n+1-i ) = temp
287 140 CONTINUE
288 END IF
289*
290 END IF
291*
292 RETURN
293*
294* End of DLATM1
295*

◆ dlatm2()

double precision function dlatm2 ( integer m,
integer n,
integer i,
integer j,
integer kl,
integer ku,
integer idist,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) d,
integer igrade,
double precision, dimension( * ) dl,
double precision, dimension( * ) dr,
integer ipvtng,
integer, dimension( * ) iwork,
double precision sparse )

DLATM2

Purpose:
!>
!>    DLATM2 returns the (I,J) entry of a random matrix of dimension
!>    (M, N) described by the other parameters. It is called by the
!>    DLATMR routine in order to build random test matrices. No error
!>    checking on parameters is done, because this routine is called in
!>    a tight loop by DLATMR which has already checked the parameters.
!>
!>    Use of DLATM2 differs from SLATM3 in the order in which the random
!>    number generator is called to fill in random matrix entries.
!>    With DLATM2, the generator is called to fill in the pivoted matrix
!>    columnwise. With DLATM3, the generator is called to fill in the
!>    matrix columnwise, after which it is pivoted. Thus, DLATM3 can
!>    be used to construct random matrices which differ only in their
!>    order of rows and/or columns. DLATM2 is used to construct band
!>    matrices while avoiding calling the random number generator for
!>    entries outside the band (and therefore generating random numbers
!>
!>    The matrix whose (I,J) entry is returned is constructed as
!>    follows (this routine only computes one entry):
!>
!>      If I is outside (1..M) or J is outside (1..N), return zero
!>         (this is convenient for generating matrices in band format).
!>
!>      Generate a matrix A with random entries of distribution IDIST.
!>
!>      Set the diagonal to D.
!>
!>      Grade the matrix, if desired, from the left (by DL) and/or
!>         from the right (by DR or DL) as specified by IGRADE.
!>
!>      Permute, if desired, the rows and/or columns as specified by
!>         IPVTNG and IWORK.
!>
!>      Band the matrix to have lower bandwidth KL and upper
!>         bandwidth KU.
!>
!>      Set random entries to zero as specified by SPARSE.
!>
Parameters
[in]M
!>          M is INTEGER
!>           Number of rows of matrix. Not modified.
!>
[in]N
!>          N is INTEGER
!>           Number of columns of matrix. Not modified.
!>
[in]I
!>          I is INTEGER
!>           Row of entry to be returned. Not modified.
!>
[in]J
!>          J is INTEGER
!>           Column of entry to be returned. Not modified.
!>
[in]KL
!>          KL is INTEGER
!>           Lower bandwidth. Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           Upper bandwidth. Not modified.
!>
[in]IDIST
!>          IDIST is INTEGER
!>           On entry, IDIST specifies the type of distribution to be
!>           used to generate a random matrix .
!>           1 => UNIFORM( 0, 1 )
!>           2 => UNIFORM( -1, 1 )
!>           3 => NORMAL( 0, 1 )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array of dimension ( 4 )
!>           Seed for random number generator.
!>           Changed on exit.
!>
[in]D
!>          D is DOUBLE PRECISION array of dimension ( MIN( I , J ) )
!>           Diagonal entries of matrix. Not modified.
!>
[in]IGRADE
!>          IGRADE is INTEGER
!>           Specifies grading of matrix as follows:
!>           0  => no grading
!>           1  => matrix premultiplied by diag( DL )
!>           2  => matrix postmultiplied by diag( DR )
!>           3  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by diag( DR )
!>           4  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by inv( diag( DL ) )
!>           5  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by diag( DL )
!>           Not modified.
!>
[in]DL
!>          DL is DOUBLE PRECISION array ( I or J, as appropriate )
!>           Left scale factors for grading matrix.  Not modified.
!>
[in]DR
!>          DR is DOUBLE PRECISION array ( I or J, as appropriate )
!>           Right scale factors for grading matrix.  Not modified.
!>
[in]IPVTNG
!>          IPVTNG is INTEGER
!>           On entry specifies pivoting permutations as follows:
!>           0 => none.
!>           1 => row pivoting.
!>           2 => column pivoting.
!>           3 => full pivoting, i.e., on both sides.
!>           Not modified.
!>
[out]IWORK
!>          IWORK is INTEGER array ( I or J, as appropriate )
!>           This array specifies the permutation used. The
!>           row (or column) in position K was originally in
!>           position IWORK( K ).
!>           This differs from IWORK for DLATM3. Not modified.
!>
[in]SPARSE
!>          SPARSE is DOUBLE PRECISION between 0. and 1.
!>           On entry specifies the sparsity of the matrix
!>           if sparse matrix is to be generated.
!>           SPARSE should lie between 0 and 1.
!>           A uniform ( 0, 1 ) random number x is generated and
!>           compared to SPARSE; if x is larger the matrix entry
!>           is unchanged and if x is smaller the entry is set
!>           to zero. Thus on the average a fraction SPARSE of the
!>           entries will be set to zero.
!>           Not modified.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 206 of file dlatm2.f.

208*
209* -- LAPACK auxiliary routine --
210* -- LAPACK is a software package provided by Univ. of Tennessee, --
211* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212*
213* .. Scalar Arguments ..
214*
215 INTEGER I, IDIST, IGRADE, IPVTNG, J, KL, KU, M, N
216 DOUBLE PRECISION SPARSE
217* ..
218*
219* .. Array Arguments ..
220*
221 INTEGER ISEED( 4 ), IWORK( * )
222 DOUBLE PRECISION D( * ), DL( * ), DR( * )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228*
229 DOUBLE PRECISION ZERO
230 parameter( zero = 0.0d0 )
231* ..
232*
233* .. Local Scalars ..
234*
235 INTEGER ISUB, JSUB
236 DOUBLE PRECISION TEMP
237* ..
238*
239* .. External Functions ..
240*
241 DOUBLE PRECISION DLARAN, DLARND
242 EXTERNAL dlaran, dlarnd
243* ..
244*
245*-----------------------------------------------------------------------
246*
247* .. Executable Statements ..
248*
249*
250* Check for I and J in range
251*
252 IF( i.LT.1 .OR. i.GT.m .OR. j.LT.1 .OR. j.GT.n ) THEN
253 dlatm2 = zero
254 RETURN
255 END IF
256*
257* Check for banding
258*
259 IF( j.GT.i+ku .OR. j.LT.i-kl ) THEN
260 dlatm2 = zero
261 RETURN
262 END IF
263*
264* Check for sparsity
265*
266 IF( sparse.GT.zero ) THEN
267 IF( dlaran( iseed ).LT.sparse ) THEN
268 dlatm2 = zero
269 RETURN
270 END IF
271 END IF
272*
273* Compute subscripts depending on IPVTNG
274*
275 IF( ipvtng.EQ.0 ) THEN
276 isub = i
277 jsub = j
278 ELSE IF( ipvtng.EQ.1 ) THEN
279 isub = iwork( i )
280 jsub = j
281 ELSE IF( ipvtng.EQ.2 ) THEN
282 isub = i
283 jsub = iwork( j )
284 ELSE IF( ipvtng.EQ.3 ) THEN
285 isub = iwork( i )
286 jsub = iwork( j )
287 END IF
288*
289* Compute entry and grade it according to IGRADE
290*
291 IF( isub.EQ.jsub ) THEN
292 temp = d( isub )
293 ELSE
294 temp = dlarnd( idist, iseed )
295 END IF
296 IF( igrade.EQ.1 ) THEN
297 temp = temp*dl( isub )
298 ELSE IF( igrade.EQ.2 ) THEN
299 temp = temp*dr( jsub )
300 ELSE IF( igrade.EQ.3 ) THEN
301 temp = temp*dl( isub )*dr( jsub )
302 ELSE IF( igrade.EQ.4 .AND. isub.NE.jsub ) THEN
303 temp = temp*dl( isub ) / dl( jsub )
304 ELSE IF( igrade.EQ.5 ) THEN
305 temp = temp*dl( isub )*dl( jsub )
306 END IF
307 dlatm2 = temp
308 RETURN
309*
310* End of DLATM2
311*
dlatm2
double precision function dlatm2(m, n, i, j, kl, ku, idist, iseed, d, igrade, dl, dr, ipvtng, iwork, sparse)
DLATM2
Definition dlatm2.f:208

◆ dlatm3()

double precision function dlatm3 ( integer m,
integer n,
integer i,
integer j,
integer isub,
integer jsub,
integer kl,
integer ku,
integer idist,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) d,
integer igrade,
double precision, dimension( * ) dl,
double precision, dimension( * ) dr,
integer ipvtng,
integer, dimension( * ) iwork,
double precision sparse )

DLATM3

Purpose:
!>
!>    DLATM3 returns the (ISUB,JSUB) entry of a random matrix of
!>    dimension (M, N) described by the other parameters. (ISUB,JSUB)
!>    is the final position of the (I,J) entry after pivoting
!>    according to IPVTNG and IWORK. DLATM3 is called by the
!>    DLATMR routine in order to build random test matrices. No error
!>    checking on parameters is done, because this routine is called in
!>    a tight loop by DLATMR which has already checked the parameters.
!>
!>    Use of DLATM3 differs from SLATM2 in the order in which the random
!>    number generator is called to fill in random matrix entries.
!>    With DLATM2, the generator is called to fill in the pivoted matrix
!>    columnwise. With DLATM3, the generator is called to fill in the
!>    matrix columnwise, after which it is pivoted. Thus, DLATM3 can
!>    be used to construct random matrices which differ only in their
!>    order of rows and/or columns. DLATM2 is used to construct band
!>    matrices while avoiding calling the random number generator for
!>    entries outside the band (and therefore generating random numbers
!>    in different orders for different pivot orders).
!>
!>    The matrix whose (ISUB,JSUB) entry is returned is constructed as
!>    follows (this routine only computes one entry):
!>
!>      If ISUB is outside (1..M) or JSUB is outside (1..N), return zero
!>         (this is convenient for generating matrices in band format).
!>
!>      Generate a matrix A with random entries of distribution IDIST.
!>
!>      Set the diagonal to D.
!>
!>      Grade the matrix, if desired, from the left (by DL) and/or
!>         from the right (by DR or DL) as specified by IGRADE.
!>
!>      Permute, if desired, the rows and/or columns as specified by
!>         IPVTNG and IWORK.
!>
!>      Band the matrix to have lower bandwidth KL and upper
!>         bandwidth KU.
!>
!>      Set random entries to zero as specified by SPARSE.
!>
Parameters
[in]M
!>          M is INTEGER
!>           Number of rows of matrix. Not modified.
!>
[in]N
!>          N is INTEGER
!>           Number of columns of matrix. Not modified.
!>
[in]I
!>          I is INTEGER
!>           Row of unpivoted entry to be returned. Not modified.
!>
[in]J
!>          J is INTEGER
!>           Column of unpivoted entry to be returned. Not modified.
!>
[in,out]ISUB
!>          ISUB is INTEGER
!>           Row of pivoted entry to be returned. Changed on exit.
!>
[in,out]JSUB
!>          JSUB is INTEGER
!>           Column of pivoted entry to be returned. Changed on exit.
!>
[in]KL
!>          KL is INTEGER
!>           Lower bandwidth. Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           Upper bandwidth. Not modified.
!>
[in]IDIST
!>          IDIST is INTEGER
!>           On entry, IDIST specifies the type of distribution to be
!>           used to generate a random matrix .
!>           1 => UNIFORM( 0, 1 )
!>           2 => UNIFORM( -1, 1 )
!>           3 => NORMAL( 0, 1 )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array of dimension ( 4 )
!>           Seed for random number generator.
!>           Changed on exit.
!>
[in]D
!>          D is DOUBLE PRECISION array of dimension ( MIN( I , J ) )
!>           Diagonal entries of matrix. Not modified.
!>
[in]IGRADE
!>          IGRADE is INTEGER
!>           Specifies grading of matrix as follows:
!>           0  => no grading
!>           1  => matrix premultiplied by diag( DL )
!>           2  => matrix postmultiplied by diag( DR )
!>           3  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by diag( DR )
!>           4  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by inv( diag( DL ) )
!>           5  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by diag( DL )
!>           Not modified.
!>
[in]DL
!>          DL is DOUBLE PRECISION array ( I or J, as appropriate )
!>           Left scale factors for grading matrix.  Not modified.
!>
[in]DR
!>          DR is DOUBLE PRECISION array ( I or J, as appropriate )
!>           Right scale factors for grading matrix.  Not modified.
!>
[in]IPVTNG
!>          IPVTNG is INTEGER
!>           On entry specifies pivoting permutations as follows:
!>           0 => none.
!>           1 => row pivoting.
!>           2 => column pivoting.
!>           3 => full pivoting, i.e., on both sides.
!>           Not modified.
!>
[in]IWORK
!>          IWORK is INTEGER array ( I or J, as appropriate )
!>           This array specifies the permutation used. The
!>           row (or column) originally in position K is in
!>           position IWORK( K ) after pivoting.
!>           This differs from IWORK for DLATM2. Not modified.
!>
[in]SPARSE
!>          SPARSE is DOUBLE PRECISION between 0. and 1.
!>           On entry specifies the sparsity of the matrix
!>           if sparse matrix is to be generated.
!>           SPARSE should lie between 0 and 1.
!>           A uniform ( 0, 1 ) random number x is generated and
!>           compared to SPARSE; if x is larger the matrix entry
!>           is unchanged and if x is smaller the entry is set
!>           to zero. Thus on the average a fraction SPARSE of the
!>           entries will be set to zero.
!>           Not modified.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file dlatm3.f.

226*
227* -- LAPACK auxiliary routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232*
233 INTEGER I, IDIST, IGRADE, IPVTNG, ISUB, J, JSUB, KL,
234 $ KU, M, N
235 DOUBLE PRECISION SPARSE
236* ..
237*
238* .. Array Arguments ..
239*
240 INTEGER ISEED( 4 ), IWORK( * )
241 DOUBLE PRECISION D( * ), DL( * ), DR( * )
242* ..
243*
244* =====================================================================
245*
246* .. Parameters ..
247*
248 DOUBLE PRECISION ZERO
249 parameter( zero = 0.0d0 )
250* ..
251*
252* .. Local Scalars ..
253*
254 DOUBLE PRECISION TEMP
255* ..
256*
257* .. External Functions ..
258*
259 DOUBLE PRECISION DLARAN, DLARND
260 EXTERNAL dlaran, dlarnd
261* ..
262*
263*-----------------------------------------------------------------------
264*
265* .. Executable Statements ..
266*
267*
268* Check for I and J in range
269*
270 IF( i.LT.1 .OR. i.GT.m .OR. j.LT.1 .OR. j.GT.n ) THEN
271 isub = i
272 jsub = j
273 dlatm3 = zero
274 RETURN
275 END IF
276*
277* Compute subscripts depending on IPVTNG
278*
279 IF( ipvtng.EQ.0 ) THEN
280 isub = i
281 jsub = j
282 ELSE IF( ipvtng.EQ.1 ) THEN
283 isub = iwork( i )
284 jsub = j
285 ELSE IF( ipvtng.EQ.2 ) THEN
286 isub = i
287 jsub = iwork( j )
288 ELSE IF( ipvtng.EQ.3 ) THEN
289 isub = iwork( i )
290 jsub = iwork( j )
291 END IF
292*
293* Check for banding
294*
295 IF( jsub.GT.isub+ku .OR. jsub.LT.isub-kl ) THEN
296 dlatm3 = zero
297 RETURN
298 END IF
299*
300* Check for sparsity
301*
302 IF( sparse.GT.zero ) THEN
303 IF( dlaran( iseed ).LT.sparse ) THEN
304 dlatm3 = zero
305 RETURN
306 END IF
307 END IF
308*
309* Compute entry and grade it according to IGRADE
310*
311 IF( i.EQ.j ) THEN
312 temp = d( i )
313 ELSE
314 temp = dlarnd( idist, iseed )
315 END IF
316 IF( igrade.EQ.1 ) THEN
317 temp = temp*dl( i )
318 ELSE IF( igrade.EQ.2 ) THEN
319 temp = temp*dr( j )
320 ELSE IF( igrade.EQ.3 ) THEN
321 temp = temp*dl( i )*dr( j )
322 ELSE IF( igrade.EQ.4 .AND. i.NE.j ) THEN
323 temp = temp*dl( i ) / dl( j )
324 ELSE IF( igrade.EQ.5 ) THEN
325 temp = temp*dl( i )*dl( j )
326 END IF
327 dlatm3 = temp
328 RETURN
329*
330* End of DLATM3
331*
dlatm3
double precision function dlatm3(m, n, i, j, isub, jsub, kl, ku, idist, iseed, d, igrade, dl, dr, ipvtng, iwork, sparse)
DLATM3
Definition dlatm3.f:226

◆ dlatm5()

subroutine dlatm5 ( integer prtype,
integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldc, * ) c,
integer ldc,
double precision, dimension( ldd, * ) d,
integer ldd,
double precision, dimension( lde, * ) e,
integer lde,
double precision, dimension( ldf, * ) f,
integer ldf,
double precision, dimension( ldr, * ) r,
integer ldr,
double precision, dimension( ldl, * ) l,
integer ldl,
double precision alpha,
integer qblcka,
integer qblckb )

DLATM5

Purpose:
!>
!> DLATM5 generates matrices involved in the Generalized Sylvester
!> equation:
!>
!>     A * R - L * B = C
!>     D * R - L * E = F
!>
!> They also satisfy (the diagonalization condition)
!>
!>  [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] )
!>  [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] )
!>
!>
Parameters
[in]PRTYPE
!>          PRTYPE is INTEGER
!>           to a certain type of the matrices to generate
!>          (see further details).
!>
[in]M
!>          M is INTEGER
!>          Specifies the order of A and D and the number of rows in
!>          C, F,  R and L.
!>
[in]N
!>          N is INTEGER
!>          Specifies the order of B and E and the number of columns in
!>          C, F, R and L.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA, M).
!>          On exit A M-by-M is initialized according to PRTYPE.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.
!>
[out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N).
!>          On exit B N-by-N is initialized according to PRTYPE.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.
!>
[out]C
!>          C is DOUBLE PRECISION array, dimension (LDC, N).
!>          On exit C M-by-N is initialized according to PRTYPE.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of C.
!>
[out]D
!>          D is DOUBLE PRECISION array, dimension (LDD, M).
!>          On exit D M-by-M is initialized according to PRTYPE.
!>
[in]LDD
!>          LDD is INTEGER
!>          The leading dimension of D.
!>
[out]E
!>          E is DOUBLE PRECISION array, dimension (LDE, N).
!>          On exit E N-by-N is initialized according to PRTYPE.
!>
[in]LDE
!>          LDE is INTEGER
!>          The leading dimension of E.
!>
[out]F
!>          F is DOUBLE PRECISION array, dimension (LDF, N).
!>          On exit F M-by-N is initialized according to PRTYPE.
!>
[in]LDF
!>          LDF is INTEGER
!>          The leading dimension of F.
!>
[out]R
!>          R is DOUBLE PRECISION array, dimension (LDR, N).
!>          On exit R M-by-N is initialized according to PRTYPE.
!>
[in]LDR
!>          LDR is INTEGER
!>          The leading dimension of R.
!>
[out]L
!>          L is DOUBLE PRECISION array, dimension (LDL, N).
!>          On exit L M-by-N is initialized according to PRTYPE.
!>
[in]LDL
!>          LDL is INTEGER
!>          The leading dimension of L.
!>
[in]ALPHA
!>          ALPHA is DOUBLE PRECISION
!>          Parameter used in generating PRTYPE = 1 and 5 matrices.
!>
[in]QBLCKA
!>          QBLCKA is INTEGER
!>          When PRTYPE = 3, specifies the distance between 2-by-2
!>          blocks on the diagonal in A. Otherwise, QBLCKA is not
!>          referenced. QBLCKA > 1.
!>
[in]QBLCKB
!>          QBLCKB is INTEGER
!>          When PRTYPE = 3, specifies the distance between 2-by-2
!>          blocks on the diagonal in B. Otherwise, QBLCKB is not
!>          referenced. QBLCKB > 1.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
!>
!>             A : if (i == j) then A(i, j) = 1.0
!>                 if (j == i + 1) then A(i, j) = -1.0
!>                 else A(i, j) = 0.0,            i, j = 1...M
!>
!>             B : if (i == j) then B(i, j) = 1.0 - ALPHA
!>                 if (j == i + 1) then B(i, j) = 1.0
!>                 else B(i, j) = 0.0,            i, j = 1...N
!>
!>             D : if (i == j) then D(i, j) = 1.0
!>                 else D(i, j) = 0.0,            i, j = 1...M
!>
!>             E : if (i == j) then E(i, j) = 1.0
!>                 else E(i, j) = 0.0,            i, j = 1...N
!>
!>             L =  R are chosen from [-10...10],
!>                  which specifies the right hand sides (C, F).
!>
!>  PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
!>
!>             A : if (i <= j) then A(i, j) = [-1...1]
!>                 else A(i, j) = 0.0,             i, j = 1...M
!>
!>                 if (PRTYPE = 3) then
!>                    A(k + 1, k + 1) = A(k, k)
!>                    A(k + 1, k) = [-1...1]
!>                    sign(A(k, k + 1) = -(sin(A(k + 1, k))
!>                        k = 1, M - 1, QBLCKA
!>
!>             B : if (i <= j) then B(i, j) = [-1...1]
!>                 else B(i, j) = 0.0,            i, j = 1...N
!>
!>                 if (PRTYPE = 3) then
!>                    B(k + 1, k + 1) = B(k, k)
!>                    B(k + 1, k) = [-1...1]
!>                    sign(B(k, k + 1) = -(sign(B(k + 1, k))
!>                        k = 1, N - 1, QBLCKB
!>
!>             D : if (i <= j) then D(i, j) = [-1...1].
!>                 else D(i, j) = 0.0,            i, j = 1...M
!>
!>
!>             E : if (i <= j) then D(i, j) = [-1...1]
!>                 else E(i, j) = 0.0,            i, j = 1...N
!>
!>                 L, R are chosen from [-10...10],
!>                 which specifies the right hand sides (C, F).
!>
!>  PRTYPE = 4 Full
!>             A(i, j) = [-10...10]
!>             D(i, j) = [-1...1]    i,j = 1...M
!>             B(i, j) = [-10...10]
!>             E(i, j) = [-1...1]    i,j = 1...N
!>             R(i, j) = [-10...10]
!>             L(i, j) = [-1...1]    i = 1..M ,j = 1...N
!>
!>             L, R specifies the right hand sides (C, F).
!>
!>  PRTYPE = 5 special case common and/or close eigs.
!>

Definition at line 265 of file dlatm5.f.

268*
269* -- LAPACK computational routine --
270* -- LAPACK is a software package provided by Univ. of Tennessee, --
271* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
272*
273* .. Scalar Arguments ..
274 INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
275 $ PRTYPE, QBLCKA, QBLCKB
276 DOUBLE PRECISION ALPHA
277* ..
278* .. Array Arguments ..
279 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
280 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
281 $ L( LDL, * ), R( LDR, * )
282* ..
283*
284* =====================================================================
285*
286* .. Parameters ..
287 DOUBLE PRECISION ONE, ZERO, TWENTY, HALF, TWO
288 parameter( one = 1.0d+0, zero = 0.0d+0, twenty = 2.0d+1,
289 $ half = 0.5d+0, two = 2.0d+0 )
290* ..
291* .. Local Scalars ..
292 INTEGER I, J, K
293 DOUBLE PRECISION IMEPS, REEPS
294* ..
295* .. Intrinsic Functions ..
296 INTRINSIC dble, mod, sin
297* ..
298* .. External Subroutines ..
299 EXTERNAL dgemm
300* ..
301* .. Executable Statements ..
302*
303 IF( prtype.EQ.1 ) THEN
304 DO 20 i = 1, m
305 DO 10 j = 1, m
306 IF( i.EQ.j ) THEN
307 a( i, j ) = one
308 d( i, j ) = one
309 ELSE IF( i.EQ.j-1 ) THEN
310 a( i, j ) = -one
311 d( i, j ) = zero
312 ELSE
313 a( i, j ) = zero
314 d( i, j ) = zero
315 END IF
316 10 CONTINUE
317 20 CONTINUE
318*
319 DO 40 i = 1, n
320 DO 30 j = 1, n
321 IF( i.EQ.j ) THEN
322 b( i, j ) = one - alpha
323 e( i, j ) = one
324 ELSE IF( i.EQ.j-1 ) THEN
325 b( i, j ) = one
326 e( i, j ) = zero
327 ELSE
328 b( i, j ) = zero
329 e( i, j ) = zero
330 END IF
331 30 CONTINUE
332 40 CONTINUE
333*
334 DO 60 i = 1, m
335 DO 50 j = 1, n
336 r( i, j ) = ( half-sin( dble( i / j ) ) )*twenty
337 l( i, j ) = r( i, j )
338 50 CONTINUE
339 60 CONTINUE
340*
341 ELSE IF( prtype.EQ.2 .OR. prtype.EQ.3 ) THEN
342 DO 80 i = 1, m
343 DO 70 j = 1, m
344 IF( i.LE.j ) THEN
345 a( i, j ) = ( half-sin( dble( i ) ) )*two
346 d( i, j ) = ( half-sin( dble( i*j ) ) )*two
347 ELSE
348 a( i, j ) = zero
349 d( i, j ) = zero
350 END IF
351 70 CONTINUE
352 80 CONTINUE
353*
354 DO 100 i = 1, n
355 DO 90 j = 1, n
356 IF( i.LE.j ) THEN
357 b( i, j ) = ( half-sin( dble( i+j ) ) )*two
358 e( i, j ) = ( half-sin( dble( j ) ) )*two
359 ELSE
360 b( i, j ) = zero
361 e( i, j ) = zero
362 END IF
363 90 CONTINUE
364 100 CONTINUE
365*
366 DO 120 i = 1, m
367 DO 110 j = 1, n
368 r( i, j ) = ( half-sin( dble( i*j ) ) )*twenty
369 l( i, j ) = ( half-sin( dble( i+j ) ) )*twenty
370 110 CONTINUE
371 120 CONTINUE
372*
373 IF( prtype.EQ.3 ) THEN
374 IF( qblcka.LE.1 )
375 $ qblcka = 2
376 DO 130 k = 1, m - 1, qblcka
377 a( k+1, k+1 ) = a( k, k )
378 a( k+1, k ) = -sin( a( k, k+1 ) )
379 130 CONTINUE
380*
381 IF( qblckb.LE.1 )
382 $ qblckb = 2
383 DO 140 k = 1, n - 1, qblckb
384 b( k+1, k+1 ) = b( k, k )
385 b( k+1, k ) = -sin( b( k, k+1 ) )
386 140 CONTINUE
387 END IF
388*
389 ELSE IF( prtype.EQ.4 ) THEN
390 DO 160 i = 1, m
391 DO 150 j = 1, m
392 a( i, j ) = ( half-sin( dble( i*j ) ) )*twenty
393 d( i, j ) = ( half-sin( dble( i+j ) ) )*two
394 150 CONTINUE
395 160 CONTINUE
396*
397 DO 180 i = 1, n
398 DO 170 j = 1, n
399 b( i, j ) = ( half-sin( dble( i+j ) ) )*twenty
400 e( i, j ) = ( half-sin( dble( i*j ) ) )*two
401 170 CONTINUE
402 180 CONTINUE
403*
404 DO 200 i = 1, m
405 DO 190 j = 1, n
406 r( i, j ) = ( half-sin( dble( j / i ) ) )*twenty
407 l( i, j ) = ( half-sin( dble( i*j ) ) )*two
408 190 CONTINUE
409 200 CONTINUE
410*
411 ELSE IF( prtype.GE.5 ) THEN
412 reeps = half*two*twenty / alpha
413 imeps = ( half-two ) / alpha
414 DO 220 i = 1, m
415 DO 210 j = 1, n
416 r( i, j ) = ( half-sin( dble( i*j ) ) )*alpha / twenty
417 l( i, j ) = ( half-sin( dble( i+j ) ) )*alpha / twenty
418 210 CONTINUE
419 220 CONTINUE
420*
421 DO 230 i = 1, m
422 d( i, i ) = one
423 230 CONTINUE
424*
425 DO 240 i = 1, m
426 IF( i.LE.4 ) THEN
427 a( i, i ) = one
428 IF( i.GT.2 )
429 $ a( i, i ) = one + reeps
430 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
431 a( i, i+1 ) = imeps
432 ELSE IF( i.GT.1 ) THEN
433 a( i, i-1 ) = -imeps
434 END IF
435 ELSE IF( i.LE.8 ) THEN
436 IF( i.LE.6 ) THEN
437 a( i, i ) = reeps
438 ELSE
439 a( i, i ) = -reeps
440 END IF
441 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
442 a( i, i+1 ) = one
443 ELSE IF( i.GT.1 ) THEN
444 a( i, i-1 ) = -one
445 END IF
446 ELSE
447 a( i, i ) = one
448 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
449 a( i, i+1 ) = imeps*2
450 ELSE IF( i.GT.1 ) THEN
451 a( i, i-1 ) = -imeps*2
452 END IF
453 END IF
454 240 CONTINUE
455*
456 DO 250 i = 1, n
457 e( i, i ) = one
458 IF( i.LE.4 ) THEN
459 b( i, i ) = -one
460 IF( i.GT.2 )
461 $ b( i, i ) = one - reeps
462 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
463 b( i, i+1 ) = imeps
464 ELSE IF( i.GT.1 ) THEN
465 b( i, i-1 ) = -imeps
466 END IF
467 ELSE IF( i.LE.8 ) THEN
468 IF( i.LE.6 ) THEN
469 b( i, i ) = reeps
470 ELSE
471 b( i, i ) = -reeps
472 END IF
473 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
474 b( i, i+1 ) = one + imeps
475 ELSE IF( i.GT.1 ) THEN
476 b( i, i-1 ) = -one - imeps
477 END IF
478 ELSE
479 b( i, i ) = one - reeps
480 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
481 b( i, i+1 ) = imeps*2
482 ELSE IF( i.GT.1 ) THEN
483 b( i, i-1 ) = -imeps*2
484 END IF
485 END IF
486 250 CONTINUE
487 END IF
488*
489* Compute rhs (C, F)
490*
491 CALL dgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
492 CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
493 CALL dgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
494 CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
495*
496* End of DLATM5
497*
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2942
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187

◆ dlatm6()

subroutine dlatm6 ( integer type,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( lda, * ) b,
double precision, dimension( ldx, * ) x,
integer ldx,
double precision, dimension( ldy, * ) y,
integer ldy,
double precision alpha,
double precision beta,
double precision wx,
double precision wy,
double precision, dimension( * ) s,
double precision, dimension( * ) dif )

DLATM6

Purpose:
!>
!> DLATM6 generates test matrices for the generalized eigenvalue
!> problem, their corresponding right and left eigenvector matrices,
!> and also reciprocal condition numbers for all eigenvalues and
!> the reciprocal condition numbers of eigenvectors corresponding to
!> the 1th and 5th eigenvalues.
!>
!> Test Matrices
!> =============
!>
!> Two kinds of test matrix pairs
!>
!>       (A, B) = inverse(YH) * (Da, Db) * inverse(X)
!>
!> are used in the tests:
!>
!> Type 1:
!>    Da = 1+a   0    0    0    0    Db = 1   0   0   0   0
!>          0   2+a   0    0    0         0   1   0   0   0
!>          0    0   3+a   0    0         0   0   1   0   0
!>          0    0    0   4+a   0         0   0   0   1   0
!>          0    0    0    0   5+a ,      0   0   0   0   1 , and
!>
!> Type 2:
!>    Da =  1   -1    0    0    0    Db = 1   0   0   0   0
!>          1    1    0    0    0         0   1   0   0   0
!>          0    0    1    0    0         0   0   1   0   0
!>          0    0    0   1+a  1+b        0   0   0   1   0
!>          0    0    0  -1-b  1+a ,      0   0   0   0   1 .
!>
!> In both cases the same inverse(YH) and inverse(X) are used to compute
!> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
!>
!> YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x
!>         0    1   -y    y   -y         0   1   x  -x  -x
!>         0    0    1    0    0         0   0   1   0   0
!>         0    0    0    1    0         0   0   0   1   0
!>         0    0    0    0    1,        0   0   0   0   1 ,
!>
!> where a, b, x and y will have all values independently of each other.
!>
Parameters
[in]TYPE
!>          TYPE is INTEGER
!>          Specifies the problem type (see further details).
!>
[in]N
!>          N is INTEGER
!>          Size of the matrices A and B.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA, N).
!>          On exit A N-by-N is initialized according to TYPE.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A and of B.
!>
[out]B
!>          B is DOUBLE PRECISION array, dimension (LDA, N).
!>          On exit B N-by-N is initialized according to TYPE.
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (LDX, N).
!>          On exit X is the N-by-N matrix of right eigenvectors.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of X.
!>
[out]Y
!>          Y is DOUBLE PRECISION array, dimension (LDY, N).
!>          On exit Y is the N-by-N matrix of left eigenvectors.
!>
[in]LDY
!>          LDY is INTEGER
!>          The leading dimension of Y.
!>
[in]ALPHA
!>          ALPHA is DOUBLE PRECISION
!>
[in]BETA
!>          BETA is DOUBLE PRECISION
!>
!>          Weighting constants for matrix A.
!>
[in]WX
!>          WX is DOUBLE PRECISION
!>          Constant for right eigenvector matrix.
!>
[in]WY
!>          WY is DOUBLE PRECISION
!>          Constant for left eigenvector matrix.
!>
[out]S
!>          S is DOUBLE PRECISION array, dimension (N)
!>          S(i) is the reciprocal condition number for eigenvalue i.
!>
[out]DIF
!>          DIF is DOUBLE PRECISION array, dimension (N)
!>          DIF(i) is the reciprocal condition number for eigenvector i.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 174 of file dlatm6.f.

176*
177* -- LAPACK computational routine --
178* -- LAPACK is a software package provided by Univ. of Tennessee, --
179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181* .. Scalar Arguments ..
182 INTEGER LDA, LDX, LDY, N, TYPE
183 DOUBLE PRECISION ALPHA, BETA, WX, WY
184* ..
185* .. Array Arguments ..
186 DOUBLE PRECISION A( LDA, * ), B( LDA, * ), DIF( * ), S( * ),
187 $ X( LDX, * ), Y( LDY, * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 DOUBLE PRECISION ZERO, ONE, TWO, THREE
194 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
195 $ three = 3.0d+0 )
196* ..
197* .. Local Scalars ..
198 INTEGER I, INFO, J
199* ..
200* .. Local Arrays ..
201 DOUBLE PRECISION WORK( 100 ), Z( 12, 12 )
202* ..
203* .. Intrinsic Functions ..
204 INTRINSIC dble, sqrt
205* ..
206* .. External Subroutines ..
207 EXTERNAL dgesvd, dlacpy, dlakf2
208* ..
209* .. Executable Statements ..
210*
211* Generate test problem ...
212* (Da, Db) ...
213*
214 DO 20 i = 1, n
215 DO 10 j = 1, n
216*
217 IF( i.EQ.j ) THEN
218 a( i, i ) = dble( i ) + alpha
219 b( i, i ) = one
220 ELSE
221 a( i, j ) = zero
222 b( i, j ) = zero
223 END IF
224*
225 10 CONTINUE
226 20 CONTINUE
227*
228* Form X and Y
229*
230 CALL dlacpy( 'F', n, n, b, lda, y, ldy )
231 y( 3, 1 ) = -wy
232 y( 4, 1 ) = wy
233 y( 5, 1 ) = -wy
234 y( 3, 2 ) = -wy
235 y( 4, 2 ) = wy
236 y( 5, 2 ) = -wy
237*
238 CALL dlacpy( 'F', n, n, b, lda, x, ldx )
239 x( 1, 3 ) = -wx
240 x( 1, 4 ) = -wx
241 x( 1, 5 ) = wx
242 x( 2, 3 ) = wx
243 x( 2, 4 ) = -wx
244 x( 2, 5 ) = -wx
245*
246* Form (A, B)
247*
248 b( 1, 3 ) = wx + wy
249 b( 2, 3 ) = -wx + wy
250 b( 1, 4 ) = wx - wy
251 b( 2, 4 ) = wx - wy
252 b( 1, 5 ) = -wx + wy
253 b( 2, 5 ) = wx + wy
254 IF( type.EQ.1 ) THEN
255 a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
256 a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
257 a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
258 a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
259 a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
260 a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
261 ELSE IF( type.EQ.2 ) THEN
262 a( 1, 3 ) = two*wx + wy
263 a( 2, 3 ) = wy
264 a( 1, 4 ) = -wy*( two+alpha+beta )
265 a( 2, 4 ) = two*wx - wy*( two+alpha+beta )
266 a( 1, 5 ) = -two*wx + wy*( alpha-beta )
267 a( 2, 5 ) = wy*( alpha-beta )
268 a( 1, 1 ) = one
269 a( 1, 2 ) = -one
270 a( 2, 1 ) = one
271 a( 2, 2 ) = a( 1, 1 )
272 a( 3, 3 ) = one
273 a( 4, 4 ) = one + alpha
274 a( 4, 5 ) = one + beta
275 a( 5, 4 ) = -a( 4, 5 )
276 a( 5, 5 ) = a( 4, 4 )
277 END IF
278*
279* Compute condition numbers
280*
281 IF( type.EQ.1 ) THEN
282*
283 s( 1 ) = one / sqrt( ( one+three*wy*wy ) /
284 $ ( one+a( 1, 1 )*a( 1, 1 ) ) )
285 s( 2 ) = one / sqrt( ( one+three*wy*wy ) /
286 $ ( one+a( 2, 2 )*a( 2, 2 ) ) )
287 s( 3 ) = one / sqrt( ( one+two*wx*wx ) /
288 $ ( one+a( 3, 3 )*a( 3, 3 ) ) )
289 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
290 $ ( one+a( 4, 4 )*a( 4, 4 ) ) )
291 s( 5 ) = one / sqrt( ( one+two*wx*wx ) /
292 $ ( one+a( 5, 5 )*a( 5, 5 ) ) )
293*
294 CALL dlakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 12 )
295 CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
296 $ work( 10 ), 1, work( 11 ), 40, info )
297 dif( 1 ) = work( 8 )
298*
299 CALL dlakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 12 )
300 CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
301 $ work( 10 ), 1, work( 11 ), 40, info )
302 dif( 5 ) = work( 8 )
303*
304 ELSE IF( type.EQ.2 ) THEN
305*
306 s( 1 ) = one / sqrt( one / three+wy*wy )
307 s( 2 ) = s( 1 )
308 s( 3 ) = one / sqrt( one / two+wx*wx )
309 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
310 $ ( one+( one+alpha )*( one+alpha )+( one+beta )*( one+
311 $ beta ) ) )
312 s( 5 ) = s( 4 )
313*
314 CALL dlakf2( 2, 3, a, lda, a( 3, 3 ), b, b( 3, 3 ), z, 12 )
315 CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
316 $ work( 14 ), 1, work( 15 ), 60, info )
317 dif( 1 ) = work( 12 )
318*
319 CALL dlakf2( 3, 2, a, lda, a( 4, 4 ), b, b( 4, 4 ), z, 12 )
320 CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
321 $ work( 14 ), 1, work( 15 ), 60, info )
322 dif( 5 ) = work( 12 )
323*
324 END IF
325*
326 RETURN
327*
328* End of DLATM6
329*
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dgesvd
subroutine dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition dgesvd.f:211
dlakf2
subroutine dlakf2(m, n, a, lda, b, d, e, z, ldz)
DLAKF2
Definition dlakf2.f:105

◆ dlatm7()

subroutine dlatm7 ( integer mode,
double precision cond,
integer irsign,
integer idist,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) d,
integer n,
integer rank,
integer info )

DLATM7

Purpose:
!>
!>    DLATM7 computes the entries of D as specified by MODE
!>    COND and IRSIGN. IDIST and ISEED determine the generation
!>    of random numbers. DLATM7 is called by DLATMT to generate
!>    random test matrices.
!> 
!>  MODE   - INTEGER
!>           On entry describes how D is to be computed:
!>           MODE = 0 means do not change D.
!>
!>           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
!>           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1)) I=1:RANK
!>
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           Not modified.
!>
!>  COND   - DOUBLE PRECISION
!>           On entry, used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
!>  IRSIGN - INTEGER
!>           On entry, if MODE neither -6, 0 nor 6, determines sign of
!>           entries of D
!>           0 => leave entries of D unchanged
!>           1 => multiply each entry of D by 1 or -1 with probability .5
!>
!>  IDIST  - CHARACTER*1
!>           On entry, IDIST specifies the type of distribution to be
!>           used to generate a random matrix .
!>           1 => UNIFORM( 0, 1 )
!>           2 => UNIFORM( -1, 1 )
!>           3 => NORMAL( 0, 1 )
!>           Not modified.
!>
!>  ISEED  - INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. The random number generator uses a
!>           linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATM7
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
!>  D      - DOUBLE PRECISION array, dimension ( MIN( M , N ) )
!>           Array to be computed according to MODE, COND and IRSIGN.
!>           May be changed on exit if MODE is nonzero.
!>
!>  N      - INTEGER
!>           Number of entries of D. Not modified.
!>
!>  RANK   - INTEGER
!>           The rank of matrix to be generated for modes 1,2,3 only.
!>           D( RANK+1:N ) = 0.
!>           Not modified.
!>
!>  INFO   - INTEGER
!>            0  => normal termination
!>           -1  => if MODE not in range -6 to 6
!>           -2  => if MODE neither -6, 0 nor 6, and
!>                  IRSIGN neither 0 nor 1
!>           -3  => if MODE neither -6, 0 nor 6 and COND less than 1
!>           -4  => if MODE equals 6 or -6 and IDIST not in range 1 to 3
!>           -7  => if N negative
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 120 of file dlatm7.f.

122*
123* -- LAPACK computational routine --
124* -- LAPACK is a software package provided by Univ. of Tennessee, --
125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*
127* .. Scalar Arguments ..
128 DOUBLE PRECISION COND
129 INTEGER IDIST, INFO, IRSIGN, MODE, N, RANK
130* ..
131* .. Array Arguments ..
132 DOUBLE PRECISION D( * )
133 INTEGER ISEED( 4 )
134* ..
135*
136* =====================================================================
137*
138* .. Parameters ..
139 DOUBLE PRECISION ONE
140 parameter( one = 1.0d0 )
141 DOUBLE PRECISION ZERO
142 parameter( zero = 0.0d0 )
143 DOUBLE PRECISION HALF
144 parameter( half = 0.5d0 )
145* ..
146* .. Local Scalars ..
147 DOUBLE PRECISION ALPHA, TEMP
148 INTEGER I
149* ..
150* .. External Functions ..
151 DOUBLE PRECISION DLARAN
152 EXTERNAL dlaran
153* ..
154* .. External Subroutines ..
155 EXTERNAL dlarnv, xerbla
156* ..
157* .. Intrinsic Functions ..
158 INTRINSIC abs, dble, exp, log
159* ..
160* .. Executable Statements ..
161*
162* Decode and Test the input parameters. Initialize flags & seed.
163*
164 info = 0
165*
166* Quick return if possible
167*
168 IF( n.EQ.0 )
169 $ RETURN
170*
171* Set INFO if an error
172*
173 IF( mode.LT.-6 .OR. mode.GT.6 ) THEN
174 info = -1
175 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
176 $ ( irsign.NE.0 .AND. irsign.NE.1 ) ) THEN
177 info = -2
178 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
179 $ cond.LT.one ) THEN
180 info = -3
181 ELSE IF( ( mode.EQ.6 .OR. mode.EQ.-6 ) .AND.
182 $ ( idist.LT.1 .OR. idist.GT.3 ) ) THEN
183 info = -4
184 ELSE IF( n.LT.0 ) THEN
185 info = -7
186 END IF
187*
188 IF( info.NE.0 ) THEN
189 CALL xerbla( 'DLATM7', -info )
190 RETURN
191 END IF
192*
193* Compute D according to COND and MODE
194*
195 IF( mode.NE.0 ) THEN
196 GO TO ( 100, 130, 160, 190, 210, 230 )abs( mode )
197*
198* One large D value:
199*
200 100 CONTINUE
201 DO 110 i = 2, rank
202 d( i ) = one / cond
203 110 CONTINUE
204 DO 120 i = rank + 1, n
205 d( i ) = zero
206 120 CONTINUE
207 d( 1 ) = one
208 GO TO 240
209*
210* One small D value:
211*
212 130 CONTINUE
213 DO 140 i = 1, rank - 1
214 d( i ) = one
215 140 CONTINUE
216 DO 150 i = rank + 1, n
217 d( i ) = zero
218 150 CONTINUE
219 d( rank ) = one / cond
220 GO TO 240
221*
222* Exponentially distributed D values:
223*
224 160 CONTINUE
225 d( 1 ) = one
226 IF( n.GT.1 .AND. rank.GT.1 ) THEN
227 alpha = cond**( -one / dble( rank-1 ) )
228 DO 170 i = 2, rank
229 d( i ) = alpha**( i-1 )
230 170 CONTINUE
231 DO 180 i = rank + 1, n
232 d( i ) = zero
233 180 CONTINUE
234 END IF
235 GO TO 240
236*
237* Arithmetically distributed D values:
238*
239 190 CONTINUE
240 d( 1 ) = one
241 IF( n.GT.1 ) THEN
242 temp = one / cond
243 alpha = ( one-temp ) / dble( n-1 )
244 DO 200 i = 2, n
245 d( i ) = dble( n-i )*alpha + temp
246 200 CONTINUE
247 END IF
248 GO TO 240
249*
250* Randomly distributed D values on ( 1/COND , 1):
251*
252 210 CONTINUE
253 alpha = log( one / cond )
254 DO 220 i = 1, n
255 d( i ) = exp( alpha*dlaran( iseed ) )
256 220 CONTINUE
257 GO TO 240
258*
259* Randomly distributed D values from IDIST
260*
261 230 CONTINUE
262 CALL dlarnv( idist, iseed, n, d )
263*
264 240 CONTINUE
265*
266* If MODE neither -6 nor 0 nor 6, and IRSIGN = 1, assign
267* random signs to D
268*
269 IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
270 $ irsign.EQ.1 ) THEN
271 DO 250 i = 1, n
272 temp = dlaran( iseed )
273 IF( temp.GT.half )
274 $ d( i ) = -d( i )
275 250 CONTINUE
276 END IF
277*
278* Reverse if MODE < 0
279*
280 IF( mode.LT.0 ) THEN
281 DO 260 i = 1, n / 2
282 temp = d( i )
283 d( i ) = d( n+1-i )
284 d( n+1-i ) = temp
285 260 CONTINUE
286 END IF
287*
288 END IF
289*
290 RETURN
291*
292* End of DLATM7
293*

◆ dlatme()

subroutine dlatme ( integer n,
character dist,
integer, dimension( 4 ) iseed,
double precision, dimension( * ) d,
integer mode,
double precision cond,
double precision dmax,
character, dimension( * ) ei,
character rsign,
character upper,
character sim,
double precision, dimension( * ) ds,
integer modes,
double precision conds,
integer kl,
integer ku,
double precision anorm,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) work,
integer info )

DLATME

Purpose:
!>
!>    DLATME generates random non-symmetric square matrices with
!>    specified eigenvalues for testing LAPACK programs.
!>
!>    DLATME operates by applying the following sequence of
!>    operations:
!>
!>    1. Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX, and RSIGN
!>         as described below.
!>
!>    2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
!>         or MODE=5), certain pairs of adjacent elements of D are
!>         interpreted as the real and complex parts of a complex
!>         conjugate pair; A thus becomes block diagonal, with 1x1
!>         and 2x2 blocks.
!>
!>    3. If UPPER='T', the upper triangle of A is set to random values
!>         out of distribution DIST.
!>
!>    4. If SIM='T', A is multiplied on the left by a random matrix
!>         X, whose singular values are specified by DS, MODES, and
!>         CONDS, and on the right by X inverse.
!>
!>    5. If KL < N-1, the lower bandwidth is reduced to KL using
!>         Householder transformations.  If KU < N-1, the upper
!>         bandwidth is reduced to KU.
!>
!>    6. If ANORM is not negative, the matrix is scaled to have
!>         maximum-element-norm ANORM.
!>
!>    (Note: since the matrix cannot be reduced beyond Hessenberg form,
!>     no packing options are available.)
!>
Parameters
[in]N
!>          N is INTEGER
!>           The number of columns (or rows) of A. Not modified.
!>
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate the random eigen-/singular values, and for the
!>           upper triangle (see UPPER).
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATME
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension ( N )
!>           This array is used to specify the eigenvalues of A.  If
!>           MODE=0, then D is assumed to contain the eigenvalues (but
!>           see the description of EI), otherwise they will be
!>           computed according to MODE, COND, DMAX, and RSIGN and
!>           placed in D.
!>           Modified if MODE is nonzero.
!>
[in]MODE
!>          MODE is INTEGER
!>           On entry this describes how the eigenvalues are to
!>           be specified:
!>           MODE = 0 means use D (with EI) as input
!>           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
!>           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.  Each odd-even pair
!>                    of elements will be either used as two real
!>                    eigenvalues or as the real and imaginary part
!>                    of a complex conjugate pair of eigenvalues;
!>                    the choice of which is done is random, with
!>                    50-50 probability, for each pair.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is between 1 and 4, D has entries ranging
!>              from 1 to 1/COND, if between -1 and -4, D has entries
!>              ranging from 1/COND to 1,
!>           Not modified.
!>
[in]COND
!>          COND is DOUBLE PRECISION
!>           On entry, this is used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
[in]DMAX
!>          DMAX is DOUBLE PRECISION
!>           If MODE is neither -6, 0 nor 6, the contents of D, as
!>           computed according to MODE and COND, will be scaled by
!>           DMAX / max(abs(D(i))).  Note that DMAX need not be
!>           positive: if DMAX is negative (or zero), D will be
!>           scaled by a negative number (or zero).
!>           Not modified.
!>
[in]EI
!>          EI is CHARACTER*1 array, dimension ( N )
!>           If MODE is 0, and EI(1) is not ' ' (space character),
!>           this array specifies which elements of D (on input) are
!>           real eigenvalues and which are the real and imaginary parts
!>           of a complex conjugate pair of eigenvalues.  The elements
!>           of EI may then only have the values 'R' and 'I'.  If
!>           EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
!>           CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
!>           conjugate thereof.  If EI(j)=EI(j+1)='R', then the j-th
!>           eigenvalue is D(j) (i.e., real).  EI(1) may not be 'I',
!>           nor may two adjacent elements of EI both have the value 'I'.
!>           If MODE is not 0, then EI is ignored.  If MODE is 0 and
!>           EI(1)=' ', then the eigenvalues will all be real.
!>           Not modified.
!>
[in]RSIGN
!>          RSIGN is CHARACTER*1
!>           If MODE is not 0, 6, or -6, and RSIGN='T', then the
!>           elements of D, as computed according to MODE and COND, will
!>           be multiplied by a random sign (+1 or -1).  If RSIGN='F',
!>           they will not be.  RSIGN may only have the values 'T' or
!>           'F'.
!>           Not modified.
!>
[in]UPPER
!>          UPPER is CHARACTER*1
!>           If UPPER='T', then the elements of A above the diagonal
!>           (and above the 2x2 diagonal blocks, if A has complex
!>           eigenvalues) will be set to random numbers out of DIST.
!>           If UPPER='F', they will not.  UPPER may only have the
!>           values 'T' or 'F'.
!>           Not modified.
!>
[in]SIM
!>          SIM is CHARACTER*1
!>           If SIM='T', then A will be operated on by a , i.e., multiplied on the left by a matrix X and
!>           on the right by X inverse.  X = U S V, where U and V are
!>           random unitary matrices and S is a (diagonal) matrix of
!>           singular values specified by DS, MODES, and CONDS.  If
!>           SIM='F', then A will not be transformed.
!>           Not modified.
!>
[in,out]DS
!>          DS is DOUBLE PRECISION array, dimension ( N )
!>           This array is used to specify the singular values of X,
!>           in the same way that D specifies the eigenvalues of A.
!>           If MODE=0, the DS contains the singular values, which
!>           may not be zero.
!>           Modified if MODE is nonzero.
!>
[in]MODES
!>          MODES is INTEGER
!>
[in]CONDS
!>          CONDS is DOUBLE PRECISION
!>           Same as MODE and COND, but for specifying the diagonal
!>           of S.  MODES=-6 and +6 are not allowed (since they would
!>           result in randomly ill-conditioned eigenvalues.)
!>
[in]KL
!>          KL is INTEGER
!>           This specifies the lower bandwidth of the  matrix.  KL=1
!>           specifies upper Hessenberg form.  If KL is at least N-1,
!>           then A will have full lower bandwidth.  KL must be at
!>           least 1.
!>           Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           This specifies the upper bandwidth of the  matrix.  KU=1
!>           specifies lower Hessenberg form.  If KU is at least N-1,
!>           then A will have full upper bandwidth; if KU and KL
!>           are both at least N-1, then A will be dense.  Only one of
!>           KU and KL may be less than N-1.  KU must be at least 1.
!>           Not modified.
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>           If ANORM is not negative, then A will be scaled by a non-
!>           negative real number to make the maximum-element-norm of A
!>           to be ANORM.
!>           Not modified.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension ( LDA, N )
!>           On exit A is the desired test matrix.
!>           Modified.
!>
[in]LDA
!>          LDA is INTEGER
!>           LDA specifies the first dimension of A as declared in the
!>           calling program.  LDA must be at least N.
!>           Not modified.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension ( 3*N )
!>           Workspace.
!>           Modified.
!>
[out]INFO
!>          INFO is INTEGER
!>           Error code.  On exit, INFO will be set to one of the
!>           following values:
!>             0 => normal return
!>            -1 => N negative
!>            -2 => DIST illegal string
!>            -5 => MODE not in range -6 to 6
!>            -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>            -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
!>                  two adjacent elements of EI are 'I'.
!>            -9 => RSIGN is not 'T' or 'F'
!>           -10 => UPPER is not 'T' or 'F'
!>           -11 => SIM   is not 'T' or 'F'
!>           -12 => MODES=0 and DS has a zero singular value.
!>           -13 => MODES is not in the range -5 to 5.
!>           -14 => MODES is nonzero and CONDS is less than 1.
!>           -15 => KL is less than 1.
!>           -16 => KU is less than 1, or KL and KU are both less than
!>                  N-1.
!>           -19 => LDA is less than N.
!>            1  => Error return from DLATM1 (computing D)
!>            2  => Cannot scale to DMAX (max. eigenvalue is 0)
!>            3  => Error return from DLATM1 (computing DS)
!>            4  => Error return from DLARGE
!>            5  => Zero singular value from DLATM1.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 327 of file dlatme.f.

332*
333* -- LAPACK computational routine --
334* -- LAPACK is a software package provided by Univ. of Tennessee, --
335* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336*
337* .. Scalar Arguments ..
338 CHARACTER DIST, RSIGN, SIM, UPPER
339 INTEGER INFO, KL, KU, LDA, MODE, MODES, N
340 DOUBLE PRECISION ANORM, COND, CONDS, DMAX
341* ..
342* .. Array Arguments ..
343 CHARACTER EI( * )
344 INTEGER ISEED( 4 )
345 DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * )
346* ..
347*
348* =====================================================================
349*
350* .. Parameters ..
351 DOUBLE PRECISION ZERO
352 parameter( zero = 0.0d0 )
353 DOUBLE PRECISION ONE
354 parameter( one = 1.0d0 )
355 DOUBLE PRECISION HALF
356 parameter( half = 1.0d0 / 2.0d0 )
357* ..
358* .. Local Scalars ..
359 LOGICAL BADEI, BADS, USEEI
360 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
361 $ ISIM, IUPPER, J, JC, JCR, JR
362 DOUBLE PRECISION ALPHA, TAU, TEMP, XNORMS
363* ..
364* .. Local Arrays ..
365 DOUBLE PRECISION TEMPA( 1 )
366* ..
367* .. External Functions ..
368 LOGICAL LSAME
369 DOUBLE PRECISION DLANGE, DLARAN
370 EXTERNAL lsame, dlange, dlaran
371* ..
372* .. External Subroutines ..
373 EXTERNAL dcopy, dgemv, dger, dlarfg, dlarge, dlarnv,
374 $ dlaset, dlatm1, dscal, xerbla
375* ..
376* .. Intrinsic Functions ..
377 INTRINSIC abs, max, mod
378* ..
379* .. Executable Statements ..
380*
381* 1) Decode and Test the input parameters.
382* Initialize flags & seed.
383*
384 info = 0
385*
386* Quick return if possible
387*
388 IF( n.EQ.0 )
389 $ RETURN
390*
391* Decode DIST
392*
393 IF( lsame( dist, 'U' ) ) THEN
394 idist = 1
395 ELSE IF( lsame( dist, 'S' ) ) THEN
396 idist = 2
397 ELSE IF( lsame( dist, 'N' ) ) THEN
398 idist = 3
399 ELSE
400 idist = -1
401 END IF
402*
403* Check EI
404*
405 useei = .true.
406 badei = .false.
407 IF( lsame( ei( 1 ), ' ' ) .OR. mode.NE.0 ) THEN
408 useei = .false.
409 ELSE
410 IF( lsame( ei( 1 ), 'R' ) ) THEN
411 DO 10 j = 2, n
412 IF( lsame( ei( j ), 'I' ) ) THEN
413 IF( lsame( ei( j-1 ), 'I' ) )
414 $ badei = .true.
415 ELSE
416 IF( .NOT.lsame( ei( j ), 'R' ) )
417 $ badei = .true.
418 END IF
419 10 CONTINUE
420 ELSE
421 badei = .true.
422 END IF
423 END IF
424*
425* Decode RSIGN
426*
427 IF( lsame( rsign, 'T' ) ) THEN
428 irsign = 1
429 ELSE IF( lsame( rsign, 'F' ) ) THEN
430 irsign = 0
431 ELSE
432 irsign = -1
433 END IF
434*
435* Decode UPPER
436*
437 IF( lsame( upper, 'T' ) ) THEN
438 iupper = 1
439 ELSE IF( lsame( upper, 'F' ) ) THEN
440 iupper = 0
441 ELSE
442 iupper = -1
443 END IF
444*
445* Decode SIM
446*
447 IF( lsame( sim, 'T' ) ) THEN
448 isim = 1
449 ELSE IF( lsame( sim, 'F' ) ) THEN
450 isim = 0
451 ELSE
452 isim = -1
453 END IF
454*
455* Check DS, if MODES=0 and ISIM=1
456*
457 bads = .false.
458 IF( modes.EQ.0 .AND. isim.EQ.1 ) THEN
459 DO 20 j = 1, n
460 IF( ds( j ).EQ.zero )
461 $ bads = .true.
462 20 CONTINUE
463 END IF
464*
465* Set INFO if an error
466*
467 IF( n.LT.0 ) THEN
468 info = -1
469 ELSE IF( idist.EQ.-1 ) THEN
470 info = -2
471 ELSE IF( abs( mode ).GT.6 ) THEN
472 info = -5
473 ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
474 $ THEN
475 info = -6
476 ELSE IF( badei ) THEN
477 info = -8
478 ELSE IF( irsign.EQ.-1 ) THEN
479 info = -9
480 ELSE IF( iupper.EQ.-1 ) THEN
481 info = -10
482 ELSE IF( isim.EQ.-1 ) THEN
483 info = -11
484 ELSE IF( bads ) THEN
485 info = -12
486 ELSE IF( isim.EQ.1 .AND. abs( modes ).GT.5 ) THEN
487 info = -13
488 ELSE IF( isim.EQ.1 .AND. modes.NE.0 .AND. conds.LT.one ) THEN
489 info = -14
490 ELSE IF( kl.LT.1 ) THEN
491 info = -15
492 ELSE IF( ku.LT.1 .OR. ( ku.LT.n-1 .AND. kl.LT.n-1 ) ) THEN
493 info = -16
494 ELSE IF( lda.LT.max( 1, n ) ) THEN
495 info = -19
496 END IF
497*
498 IF( info.NE.0 ) THEN
499 CALL xerbla( 'DLATME', -info )
500 RETURN
501 END IF
502*
503* Initialize random number generator
504*
505 DO 30 i = 1, 4
506 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
507 30 CONTINUE
508*
509 IF( mod( iseed( 4 ), 2 ).NE.1 )
510 $ iseed( 4 ) = iseed( 4 ) + 1
511*
512* 2) Set up diagonal of A
513*
514* Compute D according to COND and MODE
515*
516 CALL dlatm1( mode, cond, irsign, idist, iseed, d, n, iinfo )
517 IF( iinfo.NE.0 ) THEN
518 info = 1
519 RETURN
520 END IF
521 IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
522*
523* Scale by DMAX
524*
525 temp = abs( d( 1 ) )
526 DO 40 i = 2, n
527 temp = max( temp, abs( d( i ) ) )
528 40 CONTINUE
529*
530 IF( temp.GT.zero ) THEN
531 alpha = dmax / temp
532 ELSE IF( dmax.NE.zero ) THEN
533 info = 2
534 RETURN
535 ELSE
536 alpha = zero
537 END IF
538*
539 CALL dscal( n, alpha, d, 1 )
540*
541 END IF
542*
543 CALL dlaset( 'Full', n, n, zero, zero, a, lda )
544 CALL dcopy( n, d, 1, a, lda+1 )
545*
546* Set up complex conjugate pairs
547*
548 IF( mode.EQ.0 ) THEN
549 IF( useei ) THEN
550 DO 50 j = 2, n
551 IF( lsame( ei( j ), 'I' ) ) THEN
552 a( j-1, j ) = a( j, j )
553 a( j, j-1 ) = -a( j, j )
554 a( j, j ) = a( j-1, j-1 )
555 END IF
556 50 CONTINUE
557 END IF
558*
559 ELSE IF( abs( mode ).EQ.5 ) THEN
560*
561 DO 60 j = 2, n, 2
562 IF( dlaran( iseed ).GT.half ) THEN
563 a( j-1, j ) = a( j, j )
564 a( j, j-1 ) = -a( j, j )
565 a( j, j ) = a( j-1, j-1 )
566 END IF
567 60 CONTINUE
568 END IF
569*
570* 3) If UPPER='T', set upper triangle of A to random numbers.
571* (but don't modify the corners of 2x2 blocks.)
572*
573 IF( iupper.NE.0 ) THEN
574 DO 70 jc = 2, n
575 IF( a( jc-1, jc ).NE.zero ) THEN
576 jr = jc - 2
577 ELSE
578 jr = jc - 1
579 END IF
580 CALL dlarnv( idist, iseed, jr, a( 1, jc ) )
581 70 CONTINUE
582 END IF
583*
584* 4) If SIM='T', apply similarity transformation.
585*
586* -1
587* Transform is X A X , where X = U S V, thus
588*
589* it is U S V A V' (1/S) U'
590*
591 IF( isim.NE.0 ) THEN
592*
593* Compute S (singular values of the eigenvector matrix)
594* according to CONDS and MODES
595*
596 CALL dlatm1( modes, conds, 0, 0, iseed, ds, n, iinfo )
597 IF( iinfo.NE.0 ) THEN
598 info = 3
599 RETURN
600 END IF
601*
602* Multiply by V and V'
603*
604 CALL dlarge( n, a, lda, iseed, work, iinfo )
605 IF( iinfo.NE.0 ) THEN
606 info = 4
607 RETURN
608 END IF
609*
610* Multiply by S and (1/S)
611*
612 DO 80 j = 1, n
613 CALL dscal( n, ds( j ), a( j, 1 ), lda )
614 IF( ds( j ).NE.zero ) THEN
615 CALL dscal( n, one / ds( j ), a( 1, j ), 1 )
616 ELSE
617 info = 5
618 RETURN
619 END IF
620 80 CONTINUE
621*
622* Multiply by U and U'
623*
624 CALL dlarge( n, a, lda, iseed, work, iinfo )
625 IF( iinfo.NE.0 ) THEN
626 info = 4
627 RETURN
628 END IF
629 END IF
630*
631* 5) Reduce the bandwidth.
632*
633 IF( kl.LT.n-1 ) THEN
634*
635* Reduce bandwidth -- kill column
636*
637 DO 90 jcr = kl + 1, n - 1
638 ic = jcr - kl
639 irows = n + 1 - jcr
640 icols = n + kl - jcr
641*
642 CALL dcopy( irows, a( jcr, ic ), 1, work, 1 )
643 xnorms = work( 1 )
644 CALL dlarfg( irows, xnorms, work( 2 ), 1, tau )
645 work( 1 ) = one
646*
647 CALL dgemv( 'T', irows, icols, one, a( jcr, ic+1 ), lda,
648 $ work, 1, zero, work( irows+1 ), 1 )
649 CALL dger( irows, icols, -tau, work, 1, work( irows+1 ), 1,
650 $ a( jcr, ic+1 ), lda )
651*
652 CALL dgemv( 'N', n, irows, one, a( 1, jcr ), lda, work, 1,
653 $ zero, work( irows+1 ), 1 )
654 CALL dger( n, irows, -tau, work( irows+1 ), 1, work, 1,
655 $ a( 1, jcr ), lda )
656*
657 a( jcr, ic ) = xnorms
658 CALL dlaset( 'Full', irows-1, 1, zero, zero, a( jcr+1, ic ),
659 $ lda )
660 90 CONTINUE
661 ELSE IF( ku.LT.n-1 ) THEN
662*
663* Reduce upper bandwidth -- kill a row at a time.
664*
665 DO 100 jcr = ku + 1, n - 1
666 ir = jcr - ku
667 irows = n + ku - jcr
668 icols = n + 1 - jcr
669*
670 CALL dcopy( icols, a( ir, jcr ), lda, work, 1 )
671 xnorms = work( 1 )
672 CALL dlarfg( icols, xnorms, work( 2 ), 1, tau )
673 work( 1 ) = one
674*
675 CALL dgemv( 'N', irows, icols, one, a( ir+1, jcr ), lda,
676 $ work, 1, zero, work( icols+1 ), 1 )
677 CALL dger( irows, icols, -tau, work( icols+1 ), 1, work, 1,
678 $ a( ir+1, jcr ), lda )
679*
680 CALL dgemv( 'C', icols, n, one, a( jcr, 1 ), lda, work, 1,
681 $ zero, work( icols+1 ), 1 )
682 CALL dger( icols, n, -tau, work, 1, work( icols+1 ), 1,
683 $ a( jcr, 1 ), lda )
684*
685 a( ir, jcr ) = xnorms
686 CALL dlaset( 'Full', 1, icols-1, zero, zero, a( ir, jcr+1 ),
687 $ lda )
688 100 CONTINUE
689 END IF
690*
691* Scale the matrix to have norm ANORM
692*
693 IF( anorm.GE.zero ) THEN
694 temp = dlange( 'M', n, n, a, lda, tempa )
695 IF( temp.GT.zero ) THEN
696 alpha = anorm / temp
697 DO 110 j = 1, n
698 CALL dscal( n, alpha, a( 1, j ), 1 )
699 110 CONTINUE
700 END IF
701 END IF
702*
703 RETURN
704*
705* End of DLATME
706*
dlange
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 ...
Definition dlange.f:114
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlatm1
subroutine dlatm1(mode, cond, irsign, idist, iseed, d, n, info)
DLATM1
Definition dlatm1.f:135
dlarge
subroutine dlarge(n, a, lda, iseed, work, info)
DLARGE
Definition dlarge.f:87
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339

◆ dlatmr()

subroutine dlatmr ( integer m,
integer n,
character dist,
integer, dimension( 4 ) iseed,
character sym,
double precision, dimension( * ) d,
integer mode,
double precision cond,
double precision dmax,
character rsign,
character grade,
double precision, dimension( * ) dl,
integer model,
double precision condl,
double precision, dimension( * ) dr,
integer moder,
double precision condr,
character pivtng,
integer, dimension( * ) ipivot,
integer kl,
integer ku,
double precision sparse,
double precision anorm,
character pack,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) iwork,
integer info )

DLATMR

Purpose:
!>
!>    DLATMR generates random matrices of various types for testing
!>    LAPACK programs.
!>
!>    DLATMR operates by applying the following sequence of
!>    operations:
!>
!>      Generate a matrix A with random entries of distribution DIST
!>         which is symmetric if SYM='S', and nonsymmetric
!>         if SYM='N'.
!>
!>      Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX and RSIGN
!>         as described below.
!>
!>      Grade the matrix, if desired, from the left and/or right
!>         as specified by GRADE. The inputs DL, MODEL, CONDL, DR,
!>         MODER and CONDR also determine the grading as described
!>         below.
!>
!>      Permute, if desired, the rows and/or columns as specified by
!>         PIVTNG and IPIVOT.
!>
!>      Set random entries to zero, if desired, to get a random sparse
!>         matrix as specified by SPARSE.
!>
!>      Make A a band matrix, if desired, by zeroing out the matrix
!>         outside a band of lower bandwidth KL and upper bandwidth KU.
!>
!>      Scale A, if desired, to have maximum entry ANORM.
!>
!>      Pack the matrix if desired. Options specified by PACK are:
!>         no packing
!>         zero out upper half (if symmetric)
!>         zero out lower half (if symmetric)
!>         store the upper half columnwise (if symmetric or
!>             square upper triangular)
!>         store the lower half columnwise (if symmetric or
!>             square lower triangular)
!>             same as upper half rowwise if symmetric
!>         store the lower triangle in banded format (if symmetric)
!>         store the upper triangle in banded format (if symmetric)
!>         store the entire matrix in banded format
!>
!>    Note: If two calls to DLATMR differ only in the PACK parameter,
!>          they will generate mathematically equivalent matrices.
!>
!>          If two calls to DLATMR both have full bandwidth (KL = M-1
!>          and KU = N-1), and differ only in the PIVTNG and PACK
!>          parameters, then the matrices generated will differ only
!>          in the order of the rows and/or columns, and otherwise
!>          contain the same data. This consistency cannot be and
!>          is not maintained with less than full bandwidth.
!>
Parameters
[in]M
!>          M is INTEGER
!>           Number of rows of A. Not modified.
!>
[in]N
!>          N is INTEGER
!>           Number of columns of A. Not modified.
!>
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate a random matrix .
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATMR
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
[in]SYM
!>          SYM is CHARACTER*1
!>           If SYM='S' or 'H', generated matrix is symmetric.
!>           If SYM='N', generated matrix is nonsymmetric.
!>           Not modified.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (min(M,N))
!>           On entry this array specifies the diagonal entries
!>           of the diagonal of A.  D may either be specified
!>           on entry, or set according to MODE and COND as described
!>           below. May be changed on exit if MODE is nonzero.
!>
[in]MODE
!>          MODE is INTEGER
!>           On entry describes how D is to be used:
!>           MODE = 0 means use D as input
!>           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
!>           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           Not modified.
!>
[in]COND
!>          COND is DOUBLE PRECISION
!>           On entry, used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
[in]DMAX
!>          DMAX is DOUBLE PRECISION
!>           If MODE neither -6, 0 nor 6, the diagonal is scaled by
!>           DMAX / max(abs(D(i))), so that maximum absolute entry
!>           of diagonal is abs(DMAX). If DMAX is negative (or zero),
!>           diagonal will be scaled by a negative number (or zero).
!>
[in]RSIGN
!>          RSIGN is CHARACTER*1
!>           If MODE neither -6, 0 nor 6, specifies sign of diagonal
!>           as follows:
!>           'T' => diagonal entries are multiplied by 1 or -1
!>                  with probability .5
!>           'F' => diagonal unchanged
!>           Not modified.
!>
[in]GRADE
!>          GRADE is CHARACTER*1
!>           Specifies grading of matrix as follows:
!>           'N'  => no grading
!>           'L'  => matrix premultiplied by diag( DL )
!>                   (only if matrix nonsymmetric)
!>           'R'  => matrix postmultiplied by diag( DR )
!>                   (only if matrix nonsymmetric)
!>           'B'  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by diag( DR )
!>                   (only if matrix nonsymmetric)
!>           'S' or 'H'  => matrix premultiplied by diag( DL ) and
!>                          postmultiplied by diag( DL )
!>                          ('S' for symmetric, or 'H' for Hermitian)
!>           'E'  => matrix premultiplied by diag( DL ) and
!>                         postmultiplied by inv( diag( DL ) )
!>                         ( 'E' for eigenvalue invariance)
!>                   (only if matrix nonsymmetric)
!>                   Note: if GRADE='E', then M must equal N.
!>           Not modified.
!>
[in,out]DL
!>          DL is DOUBLE PRECISION array, dimension (M)
!>           If MODEL=0, then on entry this array specifies the diagonal
!>           entries of a diagonal matrix used as described under GRADE
!>           above. If MODEL is not zero, then DL will be set according
!>           to MODEL and CONDL, analogous to the way D is set according
!>           to MODE and COND (except there is no DMAX parameter for DL).
!>           If GRADE='E', then DL cannot have zero entries.
!>           Not referenced if GRADE = 'N' or 'R'. Changed on exit.
!>
[in]MODEL
!>          MODEL is INTEGER
!>           This specifies how the diagonal array DL is to be computed,
!>           just as MODE specifies how D is to be computed.
!>           Not modified.
!>
[in]CONDL
!>          CONDL is DOUBLE PRECISION
!>           When MODEL is not zero, this specifies the condition number
!>           of the computed DL.  Not modified.
!>
[in,out]DR
!>          DR is DOUBLE PRECISION array, dimension (N)
!>           If MODER=0, then on entry this array specifies the diagonal
!>           entries of a diagonal matrix used as described under GRADE
!>           above. If MODER is not zero, then DR will be set according
!>           to MODER and CONDR, analogous to the way D is set according
!>           to MODE and COND (except there is no DMAX parameter for DR).
!>           Not referenced if GRADE = 'N', 'L', 'H', 'S' or 'E'.
!>           Changed on exit.
!>
[in]MODER
!>          MODER is INTEGER
!>           This specifies how the diagonal array DR is to be computed,
!>           just as MODE specifies how D is to be computed.
!>           Not modified.
!>
[in]CONDR
!>          CONDR is DOUBLE PRECISION
!>           When MODER is not zero, this specifies the condition number
!>           of the computed DR.  Not modified.
!>
[in]PIVTNG
!>          PIVTNG is CHARACTER*1
!>           On entry specifies pivoting permutations as follows:
!>           'N' or ' ' => none.
!>           'L' => left or row pivoting (matrix must be nonsymmetric).
!>           'R' => right or column pivoting (matrix must be
!>                  nonsymmetric).
!>           'B' or 'F' => both or full pivoting, i.e., on both sides.
!>                         In this case, M must equal N
!>
!>           If two calls to DLATMR both have full bandwidth (KL = M-1
!>           and KU = N-1), and differ only in the PIVTNG and PACK
!>           parameters, then the matrices generated will differ only
!>           in the order of the rows and/or columns, and otherwise
!>           contain the same data. This consistency cannot be
!>           maintained with less than full bandwidth.
!>
[in]IPIVOT
!>          IPIVOT is INTEGER array, dimension (N or M)
!>           This array specifies the permutation used.  After the
!>           basic matrix is generated, the rows, columns, or both
!>           are permuted.   If, say, row pivoting is selected, DLATMR
!>           starts with the *last* row and interchanges the M-th and
!>           IPIVOT(M)-th rows, then moves to the next-to-last row,
!>           interchanging the (M-1)-th and the IPIVOT(M-1)-th rows,
!>           and so on.  In terms of , the permutation is
!>           (1 IPIVOT(1)) (2 IPIVOT(2)) ... (M IPIVOT(M))
!>           where the rightmost cycle is applied first.  This is the
!>           *inverse* of the effect of pivoting in LINPACK.  The idea
!>           is that factoring (with pivoting) an identity matrix
!>           which has been inverse-pivoted in this way should
!>           result in a pivot vector identical to IPIVOT.
!>           Not referenced if PIVTNG = 'N'. Not modified.
!>
[in]KL
!>          KL is INTEGER
!>           On entry specifies the lower bandwidth of the  matrix. For
!>           example, KL=0 implies upper triangular, KL=1 implies upper
!>           Hessenberg, and KL at least M-1 implies the matrix is not
!>           banded. Must equal KU if matrix is symmetric.
!>           Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           On entry specifies the upper bandwidth of the  matrix. For
!>           example, KU=0 implies lower triangular, KU=1 implies lower
!>           Hessenberg, and KU at least N-1 implies the matrix is not
!>           banded. Must equal KL if matrix is symmetric.
!>           Not modified.
!>
[in]SPARSE
!>          SPARSE is DOUBLE PRECISION
!>           On entry specifies the sparsity of the matrix if a sparse
!>           matrix is to be generated. SPARSE should lie between
!>           0 and 1. To generate a sparse matrix, for each matrix entry
!>           a uniform ( 0, 1 ) random number x is generated and
!>           compared to SPARSE; if x is larger the matrix entry
!>           is unchanged and if x is smaller the entry is set
!>           to zero. Thus on the average a fraction SPARSE of the
!>           entries will be set to zero.
!>           Not modified.
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>           On entry specifies maximum entry of output matrix
!>           (output matrix will by multiplied by a constant so that
!>           its largest absolute entry equal ANORM)
!>           if ANORM is nonnegative. If ANORM is negative no scaling
!>           is done. Not modified.
!>
[in]PACK
!>          PACK is CHARACTER*1
!>           On entry specifies packing of matrix as follows:
!>           'N' => no packing
!>           'U' => zero out all subdiagonal entries (if symmetric)
!>           'L' => zero out all superdiagonal entries (if symmetric)
!>           'C' => store the upper triangle columnwise
!>                  (only if matrix symmetric or square upper triangular)
!>           'R' => store the lower triangle columnwise
!>                  (only if matrix symmetric or square lower triangular)
!>                  (same as upper half rowwise if symmetric)
!>           'B' => store the lower triangle in band storage scheme
!>                  (only if matrix symmetric)
!>           'Q' => store the upper triangle in band storage scheme
!>                  (only if matrix symmetric)
!>           'Z' => store the entire matrix in band storage scheme
!>                      (pivoting can be provided for by using this
!>                      option to store A in the trailing rows of
!>                      the allocated storage)
!>
!>           Using these options, the various LAPACK packed and banded
!>           storage schemes can be obtained:
!>           GB               - use 'Z'
!>           PB, SB or TB     - use 'B' or 'Q'
!>           PP, SP or TP     - use 'C' or 'R'
!>
!>           If two calls to DLATMR differ only in the PACK parameter,
!>           they will generate mathematically equivalent matrices.
!>           Not modified.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>           On exit A is the desired test matrix. Only those
!>           entries of A which are significant on output
!>           will be referenced (even if A is in packed or band
!>           storage format). The 'unoccupied corners' of A in
!>           band format will be zeroed out.
!>
[in]LDA
!>          LDA is INTEGER
!>           on entry LDA specifies the first dimension of A as
!>           declared in the calling program.
!>           If PACK='N', 'U' or 'L', LDA must be at least max ( 1, M ).
!>           If PACK='C' or 'R', LDA must be at least 1.
!>           If PACK='B', or 'Q', LDA must be MIN ( KU+1, N )
!>           If PACK='Z', LDA must be at least KUU+KLL+1, where
!>           KUU = MIN ( KU, N-1 ) and KLL = MIN ( KL, M-1 )
!>           Not modified.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension ( N or M)
!>           Workspace. Not referenced if PIVTNG = 'N'. Changed on exit.
!>
[out]INFO
!>          INFO is INTEGER
!>           Error parameter on exit:
!>             0 => normal return
!>            -1 => M negative or unequal to N and SYM='S' or 'H'
!>            -2 => N negative
!>            -3 => DIST illegal string
!>            -5 => SYM illegal string
!>            -7 => MODE not in range -6 to 6
!>            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>           -10 => MODE neither -6, 0 nor 6 and RSIGN illegal string
!>           -11 => GRADE illegal string, or GRADE='E' and
!>                  M not equal to N, or GRADE='L', 'R', 'B' or 'E' and
!>                  SYM = 'S' or 'H'
!>           -12 => GRADE = 'E' and DL contains zero
!>           -13 => MODEL not in range -6 to 6 and GRADE= 'L', 'B', 'H',
!>                  'S' or 'E'
!>           -14 => CONDL less than 1.0, GRADE='L', 'B', 'H', 'S' or 'E',
!>                  and MODEL neither -6, 0 nor 6
!>           -16 => MODER not in range -6 to 6 and GRADE= 'R' or 'B'
!>           -17 => CONDR less than 1.0, GRADE='R' or 'B', and
!>                  MODER neither -6, 0 nor 6
!>           -18 => PIVTNG illegal string, or PIVTNG='B' or 'F' and
!>                  M not equal to N, or PIVTNG='L' or 'R' and SYM='S'
!>                  or 'H'
!>           -19 => IPIVOT contains out of range number and
!>                  PIVTNG not equal to 'N'
!>           -20 => KL negative
!>           -21 => KU negative, or SYM='S' or 'H' and KU not equal to KL
!>           -22 => SPARSE not in range 0. to 1.
!>           -24 => PACK illegal string, or PACK='U', 'L', 'B' or 'Q'
!>                  and SYM='N', or PACK='C' and SYM='N' and either KL
!>                  not equal to 0 or N not equal to M, or PACK='R' and
!>                  SYM='N', and either KU not equal to 0 or N not equal
!>                  to M
!>           -26 => LDA too small
!>             1 => Error return from DLATM1 (computing D)
!>             2 => Cannot scale diagonal to DMAX (max. entry is 0)
!>             3 => Error return from DLATM1 (computing DL)
!>             4 => Error return from DLATM1 (computing DR)
!>             5 => ANORM is positive, but matrix constructed prior to
!>                  attempting to scale it to have norm ANORM, is zero
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 467 of file dlatmr.f.

471*
472* -- LAPACK computational routine --
473* -- LAPACK is a software package provided by Univ. of Tennessee, --
474* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
475*
476* .. Scalar Arguments ..
477 CHARACTER DIST, GRADE, PACK, PIVTNG, RSIGN, SYM
478 INTEGER INFO, KL, KU, LDA, M, MODE, MODEL, MODER, N
479 DOUBLE PRECISION ANORM, COND, CONDL, CONDR, DMAX, SPARSE
480* ..
481* .. Array Arguments ..
482 INTEGER IPIVOT( * ), ISEED( 4 ), IWORK( * )
483 DOUBLE PRECISION A( LDA, * ), D( * ), DL( * ), DR( * )
484* ..
485*
486* =====================================================================
487*
488* .. Parameters ..
489 DOUBLE PRECISION ZERO
490 parameter( zero = 0.0d0 )
491 DOUBLE PRECISION ONE
492 parameter( one = 1.0d0 )
493* ..
494* .. Local Scalars ..
495 LOGICAL BADPVT, DZERO, FULBND
496 INTEGER I, IDIST, IGRADE, IISUB, IPACK, IPVTNG, IRSIGN,
497 $ ISUB, ISYM, J, JJSUB, JSUB, K, KLL, KUU, MNMIN,
498 $ MNSUB, MXSUB, NPVTS
499 DOUBLE PRECISION ALPHA, ONORM, TEMP
500* ..
501* .. Local Arrays ..
502 DOUBLE PRECISION TEMPA( 1 )
503* ..
504* .. External Functions ..
505 LOGICAL LSAME
506 DOUBLE PRECISION DLANGB, DLANGE, DLANSB, DLANSP, DLANSY, DLATM2,
507 $ DLATM3
508 EXTERNAL lsame, dlangb, dlange, dlansb, dlansp, dlansy,
509 $ dlatm2, dlatm3
510* ..
511* .. External Subroutines ..
512 EXTERNAL dlatm1, dscal, xerbla
513* ..
514* .. Intrinsic Functions ..
515 INTRINSIC abs, max, min, mod
516* ..
517* .. Executable Statements ..
518*
519* 1) Decode and Test the input parameters.
520* Initialize flags & seed.
521*
522 info = 0
523*
524* Quick return if possible
525*
526 IF( m.EQ.0 .OR. n.EQ.0 )
527 $ RETURN
528*
529* Decode DIST
530*
531 IF( lsame( dist, 'U' ) ) THEN
532 idist = 1
533 ELSE IF( lsame( dist, 'S' ) ) THEN
534 idist = 2
535 ELSE IF( lsame( dist, 'N' ) ) THEN
536 idist = 3
537 ELSE
538 idist = -1
539 END IF
540*
541* Decode SYM
542*
543 IF( lsame( sym, 'S' ) ) THEN
544 isym = 0
545 ELSE IF( lsame( sym, 'N' ) ) THEN
546 isym = 1
547 ELSE IF( lsame( sym, 'H' ) ) THEN
548 isym = 0
549 ELSE
550 isym = -1
551 END IF
552*
553* Decode RSIGN
554*
555 IF( lsame( rsign, 'F' ) ) THEN
556 irsign = 0
557 ELSE IF( lsame( rsign, 'T' ) ) THEN
558 irsign = 1
559 ELSE
560 irsign = -1
561 END IF
562*
563* Decode PIVTNG
564*
565 IF( lsame( pivtng, 'N' ) ) THEN
566 ipvtng = 0
567 ELSE IF( lsame( pivtng, ' ' ) ) THEN
568 ipvtng = 0
569 ELSE IF( lsame( pivtng, 'L' ) ) THEN
570 ipvtng = 1
571 npvts = m
572 ELSE IF( lsame( pivtng, 'R' ) ) THEN
573 ipvtng = 2
574 npvts = n
575 ELSE IF( lsame( pivtng, 'B' ) ) THEN
576 ipvtng = 3
577 npvts = min( n, m )
578 ELSE IF( lsame( pivtng, 'F' ) ) THEN
579 ipvtng = 3
580 npvts = min( n, m )
581 ELSE
582 ipvtng = -1
583 END IF
584*
585* Decode GRADE
586*
587 IF( lsame( grade, 'N' ) ) THEN
588 igrade = 0
589 ELSE IF( lsame( grade, 'L' ) ) THEN
590 igrade = 1
591 ELSE IF( lsame( grade, 'R' ) ) THEN
592 igrade = 2
593 ELSE IF( lsame( grade, 'B' ) ) THEN
594 igrade = 3
595 ELSE IF( lsame( grade, 'E' ) ) THEN
596 igrade = 4
597 ELSE IF( lsame( grade, 'H' ) .OR. lsame( grade, 'S' ) ) THEN
598 igrade = 5
599 ELSE
600 igrade = -1
601 END IF
602*
603* Decode PACK
604*
605 IF( lsame( pack, 'N' ) ) THEN
606 ipack = 0
607 ELSE IF( lsame( pack, 'U' ) ) THEN
608 ipack = 1
609 ELSE IF( lsame( pack, 'L' ) ) THEN
610 ipack = 2
611 ELSE IF( lsame( pack, 'C' ) ) THEN
612 ipack = 3
613 ELSE IF( lsame( pack, 'R' ) ) THEN
614 ipack = 4
615 ELSE IF( lsame( pack, 'B' ) ) THEN
616 ipack = 5
617 ELSE IF( lsame( pack, 'Q' ) ) THEN
618 ipack = 6
619 ELSE IF( lsame( pack, 'Z' ) ) THEN
620 ipack = 7
621 ELSE
622 ipack = -1
623 END IF
624*
625* Set certain internal parameters
626*
627 mnmin = min( m, n )
628 kll = min( kl, m-1 )
629 kuu = min( ku, n-1 )
630*
631* If inv(DL) is used, check to see if DL has a zero entry.
632*
633 dzero = .false.
634 IF( igrade.EQ.4 .AND. model.EQ.0 ) THEN
635 DO 10 i = 1, m
636 IF( dl( i ).EQ.zero )
637 $ dzero = .true.
638 10 CONTINUE
639 END IF
640*
641* Check values in IPIVOT
642*
643 badpvt = .false.
644 IF( ipvtng.GT.0 ) THEN
645 DO 20 j = 1, npvts
646 IF( ipivot( j ).LE.0 .OR. ipivot( j ).GT.npvts )
647 $ badpvt = .true.
648 20 CONTINUE
649 END IF
650*
651* Set INFO if an error
652*
653 IF( m.LT.0 ) THEN
654 info = -1
655 ELSE IF( m.NE.n .AND. isym.EQ.0 ) THEN
656 info = -1
657 ELSE IF( n.LT.0 ) THEN
658 info = -2
659 ELSE IF( idist.EQ.-1 ) THEN
660 info = -3
661 ELSE IF( isym.EQ.-1 ) THEN
662 info = -5
663 ELSE IF( mode.LT.-6 .OR. mode.GT.6 ) THEN
664 info = -7
665 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
666 $ cond.LT.one ) THEN
667 info = -8
668 ELSE IF( ( mode.NE.-6 .AND. mode.NE.0 .AND. mode.NE.6 ) .AND.
669 $ irsign.EQ.-1 ) THEN
670 info = -10
671 ELSE IF( igrade.EQ.-1 .OR. ( igrade.EQ.4 .AND. m.NE.n ) .OR.
672 $ ( ( igrade.GE.1 .AND. igrade.LE.4 ) .AND. isym.EQ.0 ) )
673 $ THEN
674 info = -11
675 ELSE IF( igrade.EQ.4 .AND. dzero ) THEN
676 info = -12
677 ELSE IF( ( igrade.EQ.1 .OR. igrade.EQ.3 .OR. igrade.EQ.4 .OR.
678 $ igrade.EQ.5 ) .AND. ( model.LT.-6 .OR. model.GT.6 ) )
679 $ THEN
680 info = -13
681 ELSE IF( ( igrade.EQ.1 .OR. igrade.EQ.3 .OR. igrade.EQ.4 .OR.
682 $ igrade.EQ.5 ) .AND. ( model.NE.-6 .AND. model.NE.0 .AND.
683 $ model.NE.6 ) .AND. condl.LT.one ) THEN
684 info = -14
685 ELSE IF( ( igrade.EQ.2 .OR. igrade.EQ.3 ) .AND.
686 $ ( moder.LT.-6 .OR. moder.GT.6 ) ) THEN
687 info = -16
688 ELSE IF( ( igrade.EQ.2 .OR. igrade.EQ.3 ) .AND.
689 $ ( moder.NE.-6 .AND. moder.NE.0 .AND. moder.NE.6 ) .AND.
690 $ condr.LT.one ) THEN
691 info = -17
692 ELSE IF( ipvtng.EQ.-1 .OR. ( ipvtng.EQ.3 .AND. m.NE.n ) .OR.
693 $ ( ( ipvtng.EQ.1 .OR. ipvtng.EQ.2 ) .AND. isym.EQ.0 ) )
694 $ THEN
695 info = -18
696 ELSE IF( ipvtng.NE.0 .AND. badpvt ) THEN
697 info = -19
698 ELSE IF( kl.LT.0 ) THEN
699 info = -20
700 ELSE IF( ku.LT.0 .OR. ( isym.EQ.0 .AND. kl.NE.ku ) ) THEN
701 info = -21
702 ELSE IF( sparse.LT.zero .OR. sparse.GT.one ) THEN
703 info = -22
704 ELSE IF( ipack.EQ.-1 .OR. ( ( ipack.EQ.1 .OR. ipack.EQ.2 .OR.
705 $ ipack.EQ.5 .OR. ipack.EQ.6 ) .AND. isym.EQ.1 ) .OR.
706 $ ( ipack.EQ.3 .AND. isym.EQ.1 .AND. ( kl.NE.0 .OR. m.NE.
707 $ n ) ) .OR. ( ipack.EQ.4 .AND. isym.EQ.1 .AND. ( ku.NE.
708 $ 0 .OR. m.NE.n ) ) ) THEN
709 info = -24
710 ELSE IF( ( ( ipack.EQ.0 .OR. ipack.EQ.1 .OR. ipack.EQ.2 ) .AND.
711 $ lda.LT.max( 1, m ) ) .OR. ( ( ipack.EQ.3 .OR. ipack.EQ.
712 $ 4 ) .AND. lda.LT.1 ) .OR. ( ( ipack.EQ.5 .OR. ipack.EQ.
713 $ 6 ) .AND. lda.LT.kuu+1 ) .OR.
714 $ ( ipack.EQ.7 .AND. lda.LT.kll+kuu+1 ) ) THEN
715 info = -26
716 END IF
717*
718 IF( info.NE.0 ) THEN
719 CALL xerbla( 'DLATMR', -info )
720 RETURN
721 END IF
722*
723* Decide if we can pivot consistently
724*
725 fulbnd = .false.
726 IF( kuu.EQ.n-1 .AND. kll.EQ.m-1 )
727 $ fulbnd = .true.
728*
729* Initialize random number generator
730*
731 DO 30 i = 1, 4
732 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
733 30 CONTINUE
734*
735 iseed( 4 ) = 2*( iseed( 4 ) / 2 ) + 1
736*
737* 2) Set up D, DL, and DR, if indicated.
738*
739* Compute D according to COND and MODE
740*
741 CALL dlatm1( mode, cond, irsign, idist, iseed, d, mnmin, info )
742 IF( info.NE.0 ) THEN
743 info = 1
744 RETURN
745 END IF
746 IF( mode.NE.0 .AND. mode.NE.-6 .AND. mode.NE.6 ) THEN
747*
748* Scale by DMAX
749*
750 temp = abs( d( 1 ) )
751 DO 40 i = 2, mnmin
752 temp = max( temp, abs( d( i ) ) )
753 40 CONTINUE
754 IF( temp.EQ.zero .AND. dmax.NE.zero ) THEN
755 info = 2
756 RETURN
757 END IF
758 IF( temp.NE.zero ) THEN
759 alpha = dmax / temp
760 ELSE
761 alpha = one
762 END IF
763 DO 50 i = 1, mnmin
764 d( i ) = alpha*d( i )
765 50 CONTINUE
766*
767 END IF
768*
769* Compute DL if grading set
770*
771 IF( igrade.EQ.1 .OR. igrade.EQ.3 .OR. igrade.EQ.4 .OR. igrade.EQ.
772 $ 5 ) THEN
773 CALL dlatm1( model, condl, 0, idist, iseed, dl, m, info )
774 IF( info.NE.0 ) THEN
775 info = 3
776 RETURN
777 END IF
778 END IF
779*
780* Compute DR if grading set
781*
782 IF( igrade.EQ.2 .OR. igrade.EQ.3 ) THEN
783 CALL dlatm1( moder, condr, 0, idist, iseed, dr, n, info )
784 IF( info.NE.0 ) THEN
785 info = 4
786 RETURN
787 END IF
788 END IF
789*
790* 3) Generate IWORK if pivoting
791*
792 IF( ipvtng.GT.0 ) THEN
793 DO 60 i = 1, npvts
794 iwork( i ) = i
795 60 CONTINUE
796 IF( fulbnd ) THEN
797 DO 70 i = 1, npvts
798 k = ipivot( i )
799 j = iwork( i )
800 iwork( i ) = iwork( k )
801 iwork( k ) = j
802 70 CONTINUE
803 ELSE
804 DO 80 i = npvts, 1, -1
805 k = ipivot( i )
806 j = iwork( i )
807 iwork( i ) = iwork( k )
808 iwork( k ) = j
809 80 CONTINUE
810 END IF
811 END IF
812*
813* 4) Generate matrices for each kind of PACKing
814* Always sweep matrix columnwise (if symmetric, upper
815* half only) so that matrix generated does not depend
816* on PACK
817*
818 IF( fulbnd ) THEN
819*
820* Use DLATM3 so matrices generated with differing PIVOTing only
821* differ only in the order of their rows and/or columns.
822*
823 IF( ipack.EQ.0 ) THEN
824 IF( isym.EQ.0 ) THEN
825 DO 100 j = 1, n
826 DO 90 i = 1, j
827 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku,
828 $ idist, iseed, d, igrade, dl, dr, ipvtng,
829 $ iwork, sparse )
830 a( isub, jsub ) = temp
831 a( jsub, isub ) = temp
832 90 CONTINUE
833 100 CONTINUE
834 ELSE IF( isym.EQ.1 ) THEN
835 DO 120 j = 1, n
836 DO 110 i = 1, m
837 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku,
838 $ idist, iseed, d, igrade, dl, dr, ipvtng,
839 $ iwork, sparse )
840 a( isub, jsub ) = temp
841 110 CONTINUE
842 120 CONTINUE
843 END IF
844*
845 ELSE IF( ipack.EQ.1 ) THEN
846*
847 DO 140 j = 1, n
848 DO 130 i = 1, j
849 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku, idist,
850 $ iseed, d, igrade, dl, dr, ipvtng, iwork,
851 $ sparse )
852 mnsub = min( isub, jsub )
853 mxsub = max( isub, jsub )
854 a( mnsub, mxsub ) = temp
855 IF( mnsub.NE.mxsub )
856 $ a( mxsub, mnsub ) = zero
857 130 CONTINUE
858 140 CONTINUE
859*
860 ELSE IF( ipack.EQ.2 ) THEN
861*
862 DO 160 j = 1, n
863 DO 150 i = 1, j
864 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku, idist,
865 $ iseed, d, igrade, dl, dr, ipvtng, iwork,
866 $ sparse )
867 mnsub = min( isub, jsub )
868 mxsub = max( isub, jsub )
869 a( mxsub, mnsub ) = temp
870 IF( mnsub.NE.mxsub )
871 $ a( mnsub, mxsub ) = zero
872 150 CONTINUE
873 160 CONTINUE
874*
875 ELSE IF( ipack.EQ.3 ) THEN
876*
877 DO 180 j = 1, n
878 DO 170 i = 1, j
879 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku, idist,
880 $ iseed, d, igrade, dl, dr, ipvtng, iwork,
881 $ sparse )
882*
883* Compute K = location of (ISUB,JSUB) entry in packed
884* array
885*
886 mnsub = min( isub, jsub )
887 mxsub = max( isub, jsub )
888 k = mxsub*( mxsub-1 ) / 2 + mnsub
889*
890* Convert K to (IISUB,JJSUB) location
891*
892 jjsub = ( k-1 ) / lda + 1
893 iisub = k - lda*( jjsub-1 )
894*
895 a( iisub, jjsub ) = temp
896 170 CONTINUE
897 180 CONTINUE
898*
899 ELSE IF( ipack.EQ.4 ) THEN
900*
901 DO 200 j = 1, n
902 DO 190 i = 1, j
903 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku, idist,
904 $ iseed, d, igrade, dl, dr, ipvtng, iwork,
905 $ sparse )
906*
907* Compute K = location of (I,J) entry in packed array
908*
909 mnsub = min( isub, jsub )
910 mxsub = max( isub, jsub )
911 IF( mnsub.EQ.1 ) THEN
912 k = mxsub
913 ELSE
914 k = n*( n+1 ) / 2 - ( n-mnsub+1 )*( n-mnsub+2 ) /
915 $ 2 + mxsub - mnsub + 1
916 END IF
917*
918* Convert K to (IISUB,JJSUB) location
919*
920 jjsub = ( k-1 ) / lda + 1
921 iisub = k - lda*( jjsub-1 )
922*
923 a( iisub, jjsub ) = temp
924 190 CONTINUE
925 200 CONTINUE
926*
927 ELSE IF( ipack.EQ.5 ) THEN
928*
929 DO 220 j = 1, n
930 DO 210 i = j - kuu, j
931 IF( i.LT.1 ) THEN
932 a( j-i+1, i+n ) = zero
933 ELSE
934 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku,
935 $ idist, iseed, d, igrade, dl, dr, ipvtng,
936 $ iwork, sparse )
937 mnsub = min( isub, jsub )
938 mxsub = max( isub, jsub )
939 a( mxsub-mnsub+1, mnsub ) = temp
940 END IF
941 210 CONTINUE
942 220 CONTINUE
943*
944 ELSE IF( ipack.EQ.6 ) THEN
945*
946 DO 240 j = 1, n
947 DO 230 i = j - kuu, j
948 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku, idist,
949 $ iseed, d, igrade, dl, dr, ipvtng, iwork,
950 $ sparse )
951 mnsub = min( isub, jsub )
952 mxsub = max( isub, jsub )
953 a( mnsub-mxsub+kuu+1, mxsub ) = temp
954 230 CONTINUE
955 240 CONTINUE
956*
957 ELSE IF( ipack.EQ.7 ) THEN
958*
959 IF( isym.EQ.0 ) THEN
960 DO 260 j = 1, n
961 DO 250 i = j - kuu, j
962 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku,
963 $ idist, iseed, d, igrade, dl, dr, ipvtng,
964 $ iwork, sparse )
965 mnsub = min( isub, jsub )
966 mxsub = max( isub, jsub )
967 a( mnsub-mxsub+kuu+1, mxsub ) = temp
968 IF( i.LT.1 )
969 $ a( j-i+1+kuu, i+n ) = zero
970 IF( i.GE.1 .AND. mnsub.NE.mxsub )
971 $ a( mxsub-mnsub+1+kuu, mnsub ) = temp
972 250 CONTINUE
973 260 CONTINUE
974 ELSE IF( isym.EQ.1 ) THEN
975 DO 280 j = 1, n
976 DO 270 i = j - kuu, j + kll
977 temp = dlatm3( m, n, i, j, isub, jsub, kl, ku,
978 $ idist, iseed, d, igrade, dl, dr, ipvtng,
979 $ iwork, sparse )
980 a( isub-jsub+kuu+1, jsub ) = temp
981 270 CONTINUE
982 280 CONTINUE
983 END IF
984*
985 END IF
986*
987 ELSE
988*
989* Use DLATM2
990*
991 IF( ipack.EQ.0 ) THEN
992 IF( isym.EQ.0 ) THEN
993 DO 300 j = 1, n
994 DO 290 i = 1, j
995 a( i, j ) = dlatm2( m, n, i, j, kl, ku, idist,
996 $ iseed, d, igrade, dl, dr, ipvtng,
997 $ iwork, sparse )
998 a( j, i ) = a( i, j )
999 290 CONTINUE
1000 300 CONTINUE
1001 ELSE IF( isym.EQ.1 ) THEN
1002 DO 320 j = 1, n
1003 DO 310 i = 1, m
1004 a( i, j ) = dlatm2( m, n, i, j, kl, ku, idist,
1005 $ iseed, d, igrade, dl, dr, ipvtng,
1006 $ iwork, sparse )
1007 310 CONTINUE
1008 320 CONTINUE
1009 END IF
1010*
1011 ELSE IF( ipack.EQ.1 ) THEN
1012*
1013 DO 340 j = 1, n
1014 DO 330 i = 1, j
1015 a( i, j ) = dlatm2( m, n, i, j, kl, ku, idist, iseed,
1016 $ d, igrade, dl, dr, ipvtng, iwork, sparse )
1017 IF( i.NE.j )
1018 $ a( j, i ) = zero
1019 330 CONTINUE
1020 340 CONTINUE
1021*
1022 ELSE IF( ipack.EQ.2 ) THEN
1023*
1024 DO 360 j = 1, n
1025 DO 350 i = 1, j
1026 a( j, i ) = dlatm2( m, n, i, j, kl, ku, idist, iseed,
1027 $ d, igrade, dl, dr, ipvtng, iwork, sparse )
1028 IF( i.NE.j )
1029 $ a( i, j ) = zero
1030 350 CONTINUE
1031 360 CONTINUE
1032*
1033 ELSE IF( ipack.EQ.3 ) THEN
1034*
1035 isub = 0
1036 jsub = 1
1037 DO 380 j = 1, n
1038 DO 370 i = 1, j
1039 isub = isub + 1
1040 IF( isub.GT.lda ) THEN
1041 isub = 1
1042 jsub = jsub + 1
1043 END IF
1044 a( isub, jsub ) = dlatm2( m, n, i, j, kl, ku, idist,
1045 $ iseed, d, igrade, dl, dr, ipvtng,
1046 $ iwork, sparse )
1047 370 CONTINUE
1048 380 CONTINUE
1049*
1050 ELSE IF( ipack.EQ.4 ) THEN
1051*
1052 IF( isym.EQ.0 ) THEN
1053 DO 400 j = 1, n
1054 DO 390 i = 1, j
1055*
1056* Compute K = location of (I,J) entry in packed array
1057*
1058 IF( i.EQ.1 ) THEN
1059 k = j
1060 ELSE
1061 k = n*( n+1 ) / 2 - ( n-i+1 )*( n-i+2 ) / 2 +
1062 $ j - i + 1
1063 END IF
1064*
1065* Convert K to (ISUB,JSUB) location
1066*
1067 jsub = ( k-1 ) / lda + 1
1068 isub = k - lda*( jsub-1 )
1069*
1070 a( isub, jsub ) = dlatm2( m, n, i, j, kl, ku,
1071 $ idist, iseed, d, igrade, dl, dr,
1072 $ ipvtng, iwork, sparse )
1073 390 CONTINUE
1074 400 CONTINUE
1075 ELSE
1076 isub = 0
1077 jsub = 1
1078 DO 420 j = 1, n
1079 DO 410 i = j, m
1080 isub = isub + 1
1081 IF( isub.GT.lda ) THEN
1082 isub = 1
1083 jsub = jsub + 1
1084 END IF
1085 a( isub, jsub ) = dlatm2( m, n, i, j, kl, ku,
1086 $ idist, iseed, d, igrade, dl, dr,
1087 $ ipvtng, iwork, sparse )
1088 410 CONTINUE
1089 420 CONTINUE
1090 END IF
1091*
1092 ELSE IF( ipack.EQ.5 ) THEN
1093*
1094 DO 440 j = 1, n
1095 DO 430 i = j - kuu, j
1096 IF( i.LT.1 ) THEN
1097 a( j-i+1, i+n ) = zero
1098 ELSE
1099 a( j-i+1, i ) = dlatm2( m, n, i, j, kl, ku, idist,
1100 $ iseed, d, igrade, dl, dr, ipvtng,
1101 $ iwork, sparse )
1102 END IF
1103 430 CONTINUE
1104 440 CONTINUE
1105*
1106 ELSE IF( ipack.EQ.6 ) THEN
1107*
1108 DO 460 j = 1, n
1109 DO 450 i = j - kuu, j
1110 a( i-j+kuu+1, j ) = dlatm2( m, n, i, j, kl, ku, idist,
1111 $ iseed, d, igrade, dl, dr, ipvtng,
1112 $ iwork, sparse )
1113 450 CONTINUE
1114 460 CONTINUE
1115*
1116 ELSE IF( ipack.EQ.7 ) THEN
1117*
1118 IF( isym.EQ.0 ) THEN
1119 DO 480 j = 1, n
1120 DO 470 i = j - kuu, j
1121 a( i-j+kuu+1, j ) = dlatm2( m, n, i, j, kl, ku,
1122 $ idist, iseed, d, igrade, dl,
1123 $ dr, ipvtng, iwork, sparse )
1124 IF( i.LT.1 )
1125 $ a( j-i+1+kuu, i+n ) = zero
1126 IF( i.GE.1 .AND. i.NE.j )
1127 $ a( j-i+1+kuu, i ) = a( i-j+kuu+1, j )
1128 470 CONTINUE
1129 480 CONTINUE
1130 ELSE IF( isym.EQ.1 ) THEN
1131 DO 500 j = 1, n
1132 DO 490 i = j - kuu, j + kll
1133 a( i-j+kuu+1, j ) = dlatm2( m, n, i, j, kl, ku,
1134 $ idist, iseed, d, igrade, dl,
1135 $ dr, ipvtng, iwork, sparse )
1136 490 CONTINUE
1137 500 CONTINUE
1138 END IF
1139*
1140 END IF
1141*
1142 END IF
1143*
1144* 5) Scaling the norm
1145*
1146 IF( ipack.EQ.0 ) THEN
1147 onorm = dlange( 'M', m, n, a, lda, tempa )
1148 ELSE IF( ipack.EQ.1 ) THEN
1149 onorm = dlansy( 'M', 'U', n, a, lda, tempa )
1150 ELSE IF( ipack.EQ.2 ) THEN
1151 onorm = dlansy( 'M', 'L', n, a, lda, tempa )
1152 ELSE IF( ipack.EQ.3 ) THEN
1153 onorm = dlansp( 'M', 'U', n, a, tempa )
1154 ELSE IF( ipack.EQ.4 ) THEN
1155 onorm = dlansp( 'M', 'L', n, a, tempa )
1156 ELSE IF( ipack.EQ.5 ) THEN
1157 onorm = dlansb( 'M', 'L', n, kll, a, lda, tempa )
1158 ELSE IF( ipack.EQ.6 ) THEN
1159 onorm = dlansb( 'M', 'U', n, kuu, a, lda, tempa )
1160 ELSE IF( ipack.EQ.7 ) THEN
1161 onorm = dlangb( 'M', n, kll, kuu, a, lda, tempa )
1162 END IF
1163*
1164 IF( anorm.GE.zero ) THEN
1165*
1166 IF( anorm.GT.zero .AND. onorm.EQ.zero ) THEN
1167*
1168* Desired scaling impossible
1169*
1170 info = 5
1171 RETURN
1172*
1173 ELSE IF( ( anorm.GT.one .AND. onorm.LT.one ) .OR.
1174 $ ( anorm.LT.one .AND. onorm.GT.one ) ) THEN
1175*
1176* Scale carefully to avoid over / underflow
1177*
1178 IF( ipack.LE.2 ) THEN
1179 DO 510 j = 1, n
1180 CALL dscal( m, one / onorm, a( 1, j ), 1 )
1181 CALL dscal( m, anorm, a( 1, j ), 1 )
1182 510 CONTINUE
1183*
1184 ELSE IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1185*
1186 CALL dscal( n*( n+1 ) / 2, one / onorm, a, 1 )
1187 CALL dscal( n*( n+1 ) / 2, anorm, a, 1 )
1188*
1189 ELSE IF( ipack.GE.5 ) THEN
1190*
1191 DO 520 j = 1, n
1192 CALL dscal( kll+kuu+1, one / onorm, a( 1, j ), 1 )
1193 CALL dscal( kll+kuu+1, anorm, a( 1, j ), 1 )
1194 520 CONTINUE
1195*
1196 END IF
1197*
1198 ELSE
1199*
1200* Scale straightforwardly
1201*
1202 IF( ipack.LE.2 ) THEN
1203 DO 530 j = 1, n
1204 CALL dscal( m, anorm / onorm, a( 1, j ), 1 )
1205 530 CONTINUE
1206*
1207 ELSE IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1208*
1209 CALL dscal( n*( n+1 ) / 2, anorm / onorm, a, 1 )
1210*
1211 ELSE IF( ipack.GE.5 ) THEN
1212*
1213 DO 540 j = 1, n
1214 CALL dscal( kll+kuu+1, anorm / onorm, a( 1, j ), 1 )
1215 540 CONTINUE
1216 END IF
1217*
1218 END IF
1219*
1220 END IF
1221*
1222* End of DLATMR
1223*
dlangb
double precision function dlangb(norm, n, kl, ku, ab, ldab, work)
DLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlangb.f:124
dlansp
double precision function dlansp(norm, uplo, n, ap, work)
DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansp.f:114
dlansb
double precision function dlansb(norm, uplo, n, k, ab, ldab, work)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansb.f:129
dlansy
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,...
Definition dlansy.f:122

◆ dlatms()

subroutine dlatms ( integer m,
integer n,
character dist,
integer, dimension( 4 ) iseed,
character sym,
double precision, dimension( * ) d,
integer mode,
double precision cond,
double precision dmax,
integer kl,
integer ku,
character pack,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) work,
integer info )

DLATMS

Purpose:
!>
!>    DLATMS generates random matrices with specified singular values
!>    (or symmetric/hermitian with specified eigenvalues)
!>    for testing LAPACK programs.
!>
!>    DLATMS operates by applying the following sequence of
!>    operations:
!>
!>      Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX, and SYM
!>         as described below.
!>
!>      Generate a matrix with the appropriate band structure, by one
!>         of two methods:
!>
!>      Method A:
!>          Generate a dense M x N matrix by multiplying D on the left
!>              and the right by random unitary matrices, then:
!>
!>          Reduce the bandwidth according to KL and KU, using
!>          Householder transformations.
!>
!>      Method B:
!>          Convert the bandwidth-0 (i.e., diagonal) matrix to a
!>              bandwidth-1 matrix using Givens rotations, 
!>              out-of-band elements back, much as in QR; then
!>              convert the bandwidth-1 to a bandwidth-2 matrix, etc.
!>              Note that for reasonably small bandwidths (relative to
!>              M and N) this requires less storage, as a dense matrix
!>              is not generated.  Also, for symmetric matrices, only
!>              one triangle is generated.
!>
!>      Method A is chosen if the bandwidth is a large fraction of the
!>          order of the matrix, and LDA is at least M (so a dense
!>          matrix can be stored.)  Method B is chosen if the bandwidth
!>          is small (< 1/2 N for symmetric, < .3 N+M for
!>          non-symmetric), or LDA is less than M and not less than the
!>          bandwidth.
!>
!>      Pack the matrix if desired. Options specified by PACK are:
!>         no packing
!>         zero out upper half (if symmetric)
!>         zero out lower half (if symmetric)
!>         store the upper half columnwise (if symmetric or upper
!>               triangular)
!>         store the lower half columnwise (if symmetric or lower
!>               triangular)
!>         store the lower triangle in banded format (if symmetric
!>               or lower triangular)
!>         store the upper triangle in banded format (if symmetric
!>               or upper triangular)
!>         store the entire matrix in banded format
!>      If Method B is chosen, and band format is specified, then the
!>         matrix will be generated in the band format, so no repacking
!>         will be necessary.
!>
Parameters
[in]M
!>          M is INTEGER
!>           The number of rows of A. Not modified.
!>
[in]N
!>          N is INTEGER
!>           The number of columns of A. Not modified.
!>
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate the random eigen-/singular values.
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATMS
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
[in]SYM
!>          SYM is CHARACTER*1
!>           If SYM='S' or 'H', the generated matrix is symmetric, with
!>             eigenvalues specified by D, COND, MODE, and DMAX; they
!>             may be positive, negative, or zero.
!>           If SYM='P', the generated matrix is symmetric, with
!>             eigenvalues (= singular values) specified by D, COND,
!>             MODE, and DMAX; they will not be negative.
!>           If SYM='N', the generated matrix is nonsymmetric, with
!>             singular values specified by D, COND, MODE, and DMAX;
!>             they will not be negative.
!>           Not modified.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension ( MIN( M , N ) )
!>           This array is used to specify the singular values or
!>           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
!>           assumed to contain the singular/eigenvalues, otherwise
!>           they will be computed according to MODE, COND, and DMAX,
!>           and placed in D.
!>           Modified if MODE is nonzero.
!>
[in]MODE
!>          MODE is INTEGER
!>           On entry this describes how the singular/eigenvalues are to
!>           be specified:
!>           MODE = 0 means use D as input
!>           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
!>           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
!>              the elements of D will also be multiplied by a random
!>              sign (i.e., +1 or -1.)
!>           Not modified.
!>
[in]COND
!>          COND is DOUBLE PRECISION
!>           On entry, this is used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
[in]DMAX
!>          DMAX is DOUBLE PRECISION
!>           If MODE is neither -6, 0 nor 6, the contents of D, as
!>           computed according to MODE and COND, will be scaled by
!>           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
!>           singular value (which is to say the norm) will be abs(DMAX).
!>           Note that DMAX need not be positive: if DMAX is negative
!>           (or zero), D will be scaled by a negative number (or zero).
!>           Not modified.
!>
[in]KL
!>          KL is INTEGER
!>           This specifies the lower bandwidth of the  matrix. For
!>           example, KL=0 implies upper triangular, KL=1 implies upper
!>           Hessenberg, and KL being at least M-1 means that the matrix
!>           has full lower bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           This specifies the upper bandwidth of the  matrix. For
!>           example, KU=0 implies lower triangular, KU=1 implies lower
!>           Hessenberg, and KU being at least N-1 means that the matrix
!>           has full upper bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!>
[in]PACK
!>          PACK is CHARACTER*1
!>           This specifies packing of matrix as follows:
!>           'N' => no packing
!>           'U' => zero out all subdiagonal entries (if symmetric)
!>           'L' => zero out all superdiagonal entries (if symmetric)
!>           'C' => store the upper triangle columnwise
!>                  (only if the matrix is symmetric or upper triangular)
!>           'R' => store the lower triangle columnwise
!>                  (only if the matrix is symmetric or lower triangular)
!>           'B' => store the lower triangle in band storage scheme
!>                  (only if matrix symmetric or lower triangular)
!>           'Q' => store the upper triangle in band storage scheme
!>                  (only if matrix symmetric or upper triangular)
!>           'Z' => store the entire matrix in band storage scheme
!>                      (pivoting can be provided for by using this
!>                      option to store A in the trailing rows of
!>                      the allocated storage)
!>
!>           Using these options, the various LAPACK packed and banded
!>           storage schemes can be obtained:
!>           GB               - use 'Z'
!>           PB, SB or TB     - use 'B' or 'Q'
!>           PP, SP or TP     - use 'C' or 'R'
!>
!>           If two calls to DLATMS differ only in the PACK parameter,
!>           they will generate mathematically equivalent matrices.
!>           Not modified.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension ( LDA, N )
!>           On exit A is the desired test matrix.  A is first generated
!>           in full (unpacked) form, and then packed, if so specified
!>           by PACK.  Thus, the first M elements of the first N
!>           columns will always be modified.  If PACK specifies a
!>           packed or banded storage scheme, all LDA elements of the
!>           first N columns will be modified; the elements of the
!>           array which do not correspond to elements of the generated
!>           matrix are set to zero.
!>           Modified.
!>
[in]LDA
!>          LDA is INTEGER
!>           LDA specifies the first dimension of A as declared in the
!>           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
!>           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
!>           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
!>           If PACK='Z', LDA must be large enough to hold the packed
!>           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
!>           Not modified.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
!>           Workspace.
!>           Modified.
!>
[out]INFO
!>          INFO is INTEGER
!>           Error code.  On exit, INFO will be set to one of the
!>           following values:
!>             0 => normal return
!>            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
!>            -2 => N negative
!>            -3 => DIST illegal string
!>            -5 => SYM illegal string
!>            -7 => MODE not in range -6 to 6
!>            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>           -10 => KL negative
!>           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
!>           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
!>                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
!>                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
!>                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
!>                  N.
!>           -14 => LDA is less than M, or PACK='Z' and LDA is less than
!>                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
!>            1  => Error return from DLATM1
!>            2  => Cannot scale to DMAX (max. sing. value is 0)
!>            3  => Error return from DLAGGE or SLAGSY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 319 of file dlatms.f.

321*
322* -- LAPACK computational routine --
323* -- LAPACK is a software package provided by Univ. of Tennessee, --
324* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
325*
326* .. Scalar Arguments ..
327 CHARACTER DIST, PACK, SYM
328 INTEGER INFO, KL, KU, LDA, M, MODE, N
329 DOUBLE PRECISION COND, DMAX
330* ..
331* .. Array Arguments ..
332 INTEGER ISEED( 4 )
333 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
334* ..
335*
336* =====================================================================
337*
338* .. Parameters ..
339 DOUBLE PRECISION ZERO
340 parameter( zero = 0.0d0 )
341 DOUBLE PRECISION ONE
342 parameter( one = 1.0d0 )
343 DOUBLE PRECISION TWOPI
344 parameter( twopi = 6.28318530717958647692528676655900576839d+0 )
345* ..
346* .. Local Scalars ..
347 LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN
348 INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
349 $ IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
350 $ IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
351 $ JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
352 $ UUB
353 DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
354* ..
355* .. External Functions ..
356 LOGICAL LSAME
357 DOUBLE PRECISION DLARND
358 EXTERNAL lsame, dlarnd
359* ..
360* .. External Subroutines ..
361 EXTERNAL dcopy, dlagge, dlagsy, dlarot, dlartg, dlaset,
362 $ dlatm1, dscal, xerbla
363* ..
364* .. Intrinsic Functions ..
365 INTRINSIC abs, cos, dble, max, min, mod, sin
366* ..
367* .. Executable Statements ..
368*
369* 1) Decode and Test the input parameters.
370* Initialize flags & seed.
371*
372 info = 0
373*
374* Quick return if possible
375*
376 IF( m.EQ.0 .OR. n.EQ.0 )
377 $ RETURN
378*
379* Decode DIST
380*
381 IF( lsame( dist, 'U' ) ) THEN
382 idist = 1
383 ELSE IF( lsame( dist, 'S' ) ) THEN
384 idist = 2
385 ELSE IF( lsame( dist, 'N' ) ) THEN
386 idist = 3
387 ELSE
388 idist = -1
389 END IF
390*
391* Decode SYM
392*
393 IF( lsame( sym, 'N' ) ) THEN
394 isym = 1
395 irsign = 0
396 ELSE IF( lsame( sym, 'P' ) ) THEN
397 isym = 2
398 irsign = 0
399 ELSE IF( lsame( sym, 'S' ) ) THEN
400 isym = 2
401 irsign = 1
402 ELSE IF( lsame( sym, 'H' ) ) THEN
403 isym = 2
404 irsign = 1
405 ELSE
406 isym = -1
407 END IF
408*
409* Decode PACK
410*
411 isympk = 0
412 IF( lsame( pack, 'N' ) ) THEN
413 ipack = 0
414 ELSE IF( lsame( pack, 'U' ) ) THEN
415 ipack = 1
416 isympk = 1
417 ELSE IF( lsame( pack, 'L' ) ) THEN
418 ipack = 2
419 isympk = 1
420 ELSE IF( lsame( pack, 'C' ) ) THEN
421 ipack = 3
422 isympk = 2
423 ELSE IF( lsame( pack, 'R' ) ) THEN
424 ipack = 4
425 isympk = 3
426 ELSE IF( lsame( pack, 'B' ) ) THEN
427 ipack = 5
428 isympk = 3
429 ELSE IF( lsame( pack, 'Q' ) ) THEN
430 ipack = 6
431 isympk = 2
432 ELSE IF( lsame( pack, 'Z' ) ) THEN
433 ipack = 7
434 ELSE
435 ipack = -1
436 END IF
437*
438* Set certain internal parameters
439*
440 mnmin = min( m, n )
441 llb = min( kl, m-1 )
442 uub = min( ku, n-1 )
443 mr = min( m, n+llb )
444 nc = min( n, m+uub )
445*
446 IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
447 minlda = uub + 1
448 ELSE IF( ipack.EQ.7 ) THEN
449 minlda = llb + uub + 1
450 ELSE
451 minlda = m
452 END IF
453*
454* Use Givens rotation method if bandwidth small enough,
455* or if LDA is too small to store the matrix unpacked.
456*
457 givens = .false.
458 IF( isym.EQ.1 ) THEN
459 IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
460 $ givens = .true.
461 ELSE
462 IF( 2*llb.LT.m )
463 $ givens = .true.
464 END IF
465 IF( lda.LT.m .AND. lda.GE.minlda )
466 $ givens = .true.
467*
468* Set INFO if an error
469*
470 IF( m.LT.0 ) THEN
471 info = -1
472 ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
473 info = -1
474 ELSE IF( n.LT.0 ) THEN
475 info = -2
476 ELSE IF( idist.EQ.-1 ) THEN
477 info = -3
478 ELSE IF( isym.EQ.-1 ) THEN
479 info = -5
480 ELSE IF( abs( mode ).GT.6 ) THEN
481 info = -7
482 ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
483 $ THEN
484 info = -8
485 ELSE IF( kl.LT.0 ) THEN
486 info = -10
487 ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
488 info = -11
489 ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
490 $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
491 $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
492 $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
493 info = -12
494 ELSE IF( lda.LT.max( 1, minlda ) ) THEN
495 info = -14
496 END IF
497*
498 IF( info.NE.0 ) THEN
499 CALL xerbla( 'DLATMS', -info )
500 RETURN
501 END IF
502*
503* Initialize random number generator
504*
505 DO 10 i = 1, 4
506 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
507 10 CONTINUE
508*
509 IF( mod( iseed( 4 ), 2 ).NE.1 )
510 $ iseed( 4 ) = iseed( 4 ) + 1
511*
512* 2) Set up D if indicated.
513*
514* Compute D according to COND and MODE
515*
516 CALL dlatm1( mode, cond, irsign, idist, iseed, d, mnmin, iinfo )
517 IF( iinfo.NE.0 ) THEN
518 info = 1
519 RETURN
520 END IF
521*
522* Choose Top-Down if D is (apparently) increasing,
523* Bottom-Up if D is (apparently) decreasing.
524*
525 IF( abs( d( 1 ) ).LE.abs( d( mnmin ) ) ) THEN
526 topdwn = .true.
527 ELSE
528 topdwn = .false.
529 END IF
530*
531 IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
532*
533* Scale by DMAX
534*
535 temp = abs( d( 1 ) )
536 DO 20 i = 2, mnmin
537 temp = max( temp, abs( d( i ) ) )
538 20 CONTINUE
539*
540 IF( temp.GT.zero ) THEN
541 alpha = dmax / temp
542 ELSE
543 info = 2
544 RETURN
545 END IF
546*
547 CALL dscal( mnmin, alpha, d, 1 )
548*
549 END IF
550*
551* 3) Generate Banded Matrix using Givens rotations.
552* Also the special case of UUB=LLB=0
553*
554* Compute Addressing constants to cover all
555* storage formats. Whether GE, SY, GB, or SB,
556* upper or lower triangle or both,
557* the (i,j)-th element is in
558* A( i - ISKEW*j + IOFFST, j )
559*
560 IF( ipack.GT.4 ) THEN
561 ilda = lda - 1
562 iskew = 1
563 IF( ipack.GT.5 ) THEN
564 ioffst = uub + 1
565 ELSE
566 ioffst = 1
567 END IF
568 ELSE
569 ilda = lda
570 iskew = 0
571 ioffst = 0
572 END IF
573*
574* IPACKG is the format that the matrix is generated in. If this is
575* different from IPACK, then the matrix must be repacked at the
576* end. It also signals how to compute the norm, for scaling.
577*
578 ipackg = 0
579 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
580*
581* Diagonal Matrix -- We are done, unless it
582* is to be stored SP/PP/TP (PACK='R' or 'C')
583*
584 IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
585 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
586 IF( ipack.LE.2 .OR. ipack.GE.5 )
587 $ ipackg = ipack
588*
589 ELSE IF( givens ) THEN
590*
591* Check whether to use Givens rotations,
592* Householder transformations, or nothing.
593*
594 IF( isym.EQ.1 ) THEN
595*
596* Non-symmetric -- A = U D V
597*
598 IF( ipack.GT.4 ) THEN
599 ipackg = ipack
600 ELSE
601 ipackg = 0
602 END IF
603*
604 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
605*
606 IF( topdwn ) THEN
607 jkl = 0
608 DO 50 jku = 1, uub
609*
610* Transform from bandwidth JKL, JKU-1 to JKL, JKU
611*
612* Last row actually rotated is M
613* Last column actually rotated is MIN( M+JKU, N )
614*
615 DO 40 jr = 1, min( m+jku, n ) + jkl - 1
616 extra = zero
617 angle = twopi*dlarnd( 1, iseed )
618 c = cos( angle )
619 s = sin( angle )
620 icol = max( 1, jr-jkl )
621 IF( jr.LT.m ) THEN
622 il = min( n, jr+jku ) + 1 - icol
623 CALL dlarot( .true., jr.GT.jkl, .false., il, c,
624 $ s, a( jr-iskew*icol+ioffst, icol ),
625 $ ilda, extra, dummy )
626 END IF
627*
628* Chase "EXTRA" back up
629*
630 ir = jr
631 ic = icol
632 DO 30 jch = jr - jkl, 1, -jkl - jku
633 IF( ir.LT.m ) THEN
634 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
635 $ ic+1 ), extra, c, s, dummy )
636 END IF
637 irow = max( 1, jch-jku )
638 il = ir + 2 - irow
639 temp = zero
640 iltemp = jch.GT.jku
641 CALL dlarot( .false., iltemp, .true., il, c, -s,
642 $ a( irow-iskew*ic+ioffst, ic ),
643 $ ilda, temp, extra )
644 IF( iltemp ) THEN
645 CALL dlartg( a( irow+1-iskew*( ic+1 )+ioffst,
646 $ ic+1 ), temp, c, s, dummy )
647 icol = max( 1, jch-jku-jkl )
648 il = ic + 2 - icol
649 extra = zero
650 CALL dlarot( .true., jch.GT.jku+jkl, .true.,
651 $ il, c, -s, a( irow-iskew*icol+
652 $ ioffst, icol ), ilda, extra,
653 $ temp )
654 ic = icol
655 ir = irow
656 END IF
657 30 CONTINUE
658 40 CONTINUE
659 50 CONTINUE
660*
661 jku = uub
662 DO 80 jkl = 1, llb
663*
664* Transform from bandwidth JKL-1, JKU to JKL, JKU
665*
666 DO 70 jc = 1, min( n+jkl, m ) + jku - 1
667 extra = zero
668 angle = twopi*dlarnd( 1, iseed )
669 c = cos( angle )
670 s = sin( angle )
671 irow = max( 1, jc-jku )
672 IF( jc.LT.n ) THEN
673 il = min( m, jc+jkl ) + 1 - irow
674 CALL dlarot( .false., jc.GT.jku, .false., il, c,
675 $ s, a( irow-iskew*jc+ioffst, jc ),
676 $ ilda, extra, dummy )
677 END IF
678*
679* Chase "EXTRA" back up
680*
681 ic = jc
682 ir = irow
683 DO 60 jch = jc - jku, 1, -jkl - jku
684 IF( ic.LT.n ) THEN
685 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
686 $ ic+1 ), extra, c, s, dummy )
687 END IF
688 icol = max( 1, jch-jkl )
689 il = ic + 2 - icol
690 temp = zero
691 iltemp = jch.GT.jkl
692 CALL dlarot( .true., iltemp, .true., il, c, -s,
693 $ a( ir-iskew*icol+ioffst, icol ),
694 $ ilda, temp, extra )
695 IF( iltemp ) THEN
696 CALL dlartg( a( ir+1-iskew*( icol+1 )+ioffst,
697 $ icol+1 ), temp, c, s, dummy )
698 irow = max( 1, jch-jkl-jku )
699 il = ir + 2 - irow
700 extra = zero
701 CALL dlarot( .false., jch.GT.jkl+jku, .true.,
702 $ il, c, -s, a( irow-iskew*icol+
703 $ ioffst, icol ), ilda, extra,
704 $ temp )
705 ic = icol
706 ir = irow
707 END IF
708 60 CONTINUE
709 70 CONTINUE
710 80 CONTINUE
711*
712 ELSE
713*
714* Bottom-Up -- Start at the bottom right.
715*
716 jkl = 0
717 DO 110 jku = 1, uub
718*
719* Transform from bandwidth JKL, JKU-1 to JKL, JKU
720*
721* First row actually rotated is M
722* First column actually rotated is MIN( M+JKU, N )
723*
724 iendch = min( m, n+jkl ) - 1
725 DO 100 jc = min( m+jku, n ) - 1, 1 - jkl, -1
726 extra = zero
727 angle = twopi*dlarnd( 1, iseed )
728 c = cos( angle )
729 s = sin( angle )
730 irow = max( 1, jc-jku+1 )
731 IF( jc.GT.0 ) THEN
732 il = min( m, jc+jkl+1 ) + 1 - irow
733 CALL dlarot( .false., .false., jc+jkl.LT.m, il,
734 $ c, s, a( irow-iskew*jc+ioffst,
735 $ jc ), ilda, dummy, extra )
736 END IF
737*
738* Chase "EXTRA" back down
739*
740 ic = jc
741 DO 90 jch = jc + jkl, iendch, jkl + jku
742 ilextr = ic.GT.0
743 IF( ilextr ) THEN
744 CALL dlartg( a( jch-iskew*ic+ioffst, ic ),
745 $ extra, c, s, dummy )
746 END IF
747 ic = max( 1, ic )
748 icol = min( n-1, jch+jku )
749 iltemp = jch + jku.LT.n
750 temp = zero
751 CALL dlarot( .true., ilextr, iltemp, icol+2-ic,
752 $ c, s, a( jch-iskew*ic+ioffst, ic ),
753 $ ilda, extra, temp )
754 IF( iltemp ) THEN
755 CALL dlartg( a( jch-iskew*icol+ioffst,
756 $ icol ), temp, c, s, dummy )
757 il = min( iendch, jch+jkl+jku ) + 2 - jch
758 extra = zero
759 CALL dlarot( .false., .true.,
760 $ jch+jkl+jku.LE.iendch, il, c, s,
761 $ a( jch-iskew*icol+ioffst,
762 $ icol ), ilda, temp, extra )
763 ic = icol
764 END IF
765 90 CONTINUE
766 100 CONTINUE
767 110 CONTINUE
768*
769 jku = uub
770 DO 140 jkl = 1, llb
771*
772* Transform from bandwidth JKL-1, JKU to JKL, JKU
773*
774* First row actually rotated is MIN( N+JKL, M )
775* First column actually rotated is N
776*
777 iendch = min( n, m+jku ) - 1
778 DO 130 jr = min( n+jkl, m ) - 1, 1 - jku, -1
779 extra = zero
780 angle = twopi*dlarnd( 1, iseed )
781 c = cos( angle )
782 s = sin( angle )
783 icol = max( 1, jr-jkl+1 )
784 IF( jr.GT.0 ) THEN
785 il = min( n, jr+jku+1 ) + 1 - icol
786 CALL dlarot( .true., .false., jr+jku.LT.n, il,
787 $ c, s, a( jr-iskew*icol+ioffst,
788 $ icol ), ilda, dummy, extra )
789 END IF
790*
791* Chase "EXTRA" back down
792*
793 ir = jr
794 DO 120 jch = jr + jku, iendch, jkl + jku
795 ilextr = ir.GT.0
796 IF( ilextr ) THEN
797 CALL dlartg( a( ir-iskew*jch+ioffst, jch ),
798 $ extra, c, s, dummy )
799 END IF
800 ir = max( 1, ir )
801 irow = min( m-1, jch+jkl )
802 iltemp = jch + jkl.LT.m
803 temp = zero
804 CALL dlarot( .false., ilextr, iltemp, irow+2-ir,
805 $ c, s, a( ir-iskew*jch+ioffst,
806 $ jch ), ilda, extra, temp )
807 IF( iltemp ) THEN
808 CALL dlartg( a( irow-iskew*jch+ioffst, jch ),
809 $ temp, c, s, dummy )
810 il = min( iendch, jch+jkl+jku ) + 2 - jch
811 extra = zero
812 CALL dlarot( .true., .true.,
813 $ jch+jkl+jku.LE.iendch, il, c, s,
814 $ a( irow-iskew*jch+ioffst, jch ),
815 $ ilda, temp, extra )
816 ir = irow
817 END IF
818 120 CONTINUE
819 130 CONTINUE
820 140 CONTINUE
821 END IF
822*
823 ELSE
824*
825* Symmetric -- A = U D U'
826*
827 ipackg = ipack
828 ioffg = ioffst
829*
830 IF( topdwn ) THEN
831*
832* Top-Down -- Generate Upper triangle only
833*
834 IF( ipack.GE.5 ) THEN
835 ipackg = 6
836 ioffg = uub + 1
837 ELSE
838 ipackg = 1
839 END IF
840 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
841*
842 DO 170 k = 1, uub
843 DO 160 jc = 1, n - 1
844 irow = max( 1, jc-k )
845 il = min( jc+1, k+2 )
846 extra = zero
847 temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
848 angle = twopi*dlarnd( 1, iseed )
849 c = cos( angle )
850 s = sin( angle )
851 CALL dlarot( .false., jc.GT.k, .true., il, c, s,
852 $ a( irow-iskew*jc+ioffg, jc ), ilda,
853 $ extra, temp )
854 CALL dlarot( .true., .true., .false.,
855 $ min( k, n-jc )+1, c, s,
856 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
857 $ temp, dummy )
858*
859* Chase EXTRA back up the matrix
860*
861 icol = jc
862 DO 150 jch = jc - k, 1, -k
863 CALL dlartg( a( jch+1-iskew*( icol+1 )+ioffg,
864 $ icol+1 ), extra, c, s, dummy )
865 temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
866 CALL dlarot( .true., .true., .true., k+2, c, -s,
867 $ a( ( 1-iskew )*jch+ioffg, jch ),
868 $ ilda, temp, extra )
869 irow = max( 1, jch-k )
870 il = min( jch+1, k+2 )
871 extra = zero
872 CALL dlarot( .false., jch.GT.k, .true., il, c,
873 $ -s, a( irow-iskew*jch+ioffg, jch ),
874 $ ilda, extra, temp )
875 icol = jch
876 150 CONTINUE
877 160 CONTINUE
878 170 CONTINUE
879*
880* If we need lower triangle, copy from upper. Note that
881* the order of copying is chosen to work for 'q' -> 'b'
882*
883 IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
884 DO 190 jc = 1, n
885 irow = ioffst - iskew*jc
886 DO 180 jr = jc, min( n, jc+uub )
887 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
888 180 CONTINUE
889 190 CONTINUE
890 IF( ipack.EQ.5 ) THEN
891 DO 210 jc = n - uub + 1, n
892 DO 200 jr = n + 2 - jc, uub + 1
893 a( jr, jc ) = zero
894 200 CONTINUE
895 210 CONTINUE
896 END IF
897 IF( ipackg.EQ.6 ) THEN
898 ipackg = ipack
899 ELSE
900 ipackg = 0
901 END IF
902 END IF
903 ELSE
904*
905* Bottom-Up -- Generate Lower triangle only
906*
907 IF( ipack.GE.5 ) THEN
908 ipackg = 5
909 IF( ipack.EQ.6 )
910 $ ioffg = 1
911 ELSE
912 ipackg = 2
913 END IF
914 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
915*
916 DO 240 k = 1, uub
917 DO 230 jc = n - 1, 1, -1
918 il = min( n+1-jc, k+2 )
919 extra = zero
920 temp = a( 1+( 1-iskew )*jc+ioffg, jc )
921 angle = twopi*dlarnd( 1, iseed )
922 c = cos( angle )
923 s = -sin( angle )
924 CALL dlarot( .false., .true., n-jc.GT.k, il, c, s,
925 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
926 $ temp, extra )
927 icol = max( 1, jc-k+1 )
928 CALL dlarot( .true., .false., .true., jc+2-icol, c,
929 $ s, a( jc-iskew*icol+ioffg, icol ),
930 $ ilda, dummy, temp )
931*
932* Chase EXTRA back down the matrix
933*
934 icol = jc
935 DO 220 jch = jc + k, n - 1, k
936 CALL dlartg( a( jch-iskew*icol+ioffg, icol ),
937 $ extra, c, s, dummy )
938 temp = a( 1+( 1-iskew )*jch+ioffg, jch )
939 CALL dlarot( .true., .true., .true., k+2, c, s,
940 $ a( jch-iskew*icol+ioffg, icol ),
941 $ ilda, extra, temp )
942 il = min( n+1-jch, k+2 )
943 extra = zero
944 CALL dlarot( .false., .true., n-jch.GT.k, il, c,
945 $ s, a( ( 1-iskew )*jch+ioffg, jch ),
946 $ ilda, temp, extra )
947 icol = jch
948 220 CONTINUE
949 230 CONTINUE
950 240 CONTINUE
951*
952* If we need upper triangle, copy from lower. Note that
953* the order of copying is chosen to work for 'b' -> 'q'
954*
955 IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
956 DO 260 jc = n, 1, -1
957 irow = ioffst - iskew*jc
958 DO 250 jr = jc, max( 1, jc-uub ), -1
959 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
960 250 CONTINUE
961 260 CONTINUE
962 IF( ipack.EQ.6 ) THEN
963 DO 280 jc = 1, uub
964 DO 270 jr = 1, uub + 1 - jc
965 a( jr, jc ) = zero
966 270 CONTINUE
967 280 CONTINUE
968 END IF
969 IF( ipackg.EQ.5 ) THEN
970 ipackg = ipack
971 ELSE
972 ipackg = 0
973 END IF
974 END IF
975 END IF
976 END IF
977*
978 ELSE
979*
980* 4) Generate Banded Matrix by first
981* Rotating by random Unitary matrices,
982* then reducing the bandwidth using Householder
983* transformations.
984*
985* Note: we should get here only if LDA .ge. N
986*
987 IF( isym.EQ.1 ) THEN
988*
989* Non-symmetric -- A = U D V
990*
991 CALL dlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
992 $ iinfo )
993 ELSE
994*
995* Symmetric -- A = U D U'
996*
997 CALL dlagsy( m, llb, d, a, lda, iseed, work, iinfo )
998*
999 END IF
1000 IF( iinfo.NE.0 ) THEN
1001 info = 3
1002 RETURN
1003 END IF
1004 END IF
1005*
1006* 5) Pack the matrix
1007*
1008 IF( ipack.NE.ipackg ) THEN
1009 IF( ipack.EQ.1 ) THEN
1010*
1011* 'U' -- Upper triangular, not packed
1012*
1013 DO 300 j = 1, m
1014 DO 290 i = j + 1, m
1015 a( i, j ) = zero
1016 290 CONTINUE
1017 300 CONTINUE
1018*
1019 ELSE IF( ipack.EQ.2 ) THEN
1020*
1021* 'L' -- Lower triangular, not packed
1022*
1023 DO 320 j = 2, m
1024 DO 310 i = 1, j - 1
1025 a( i, j ) = zero
1026 310 CONTINUE
1027 320 CONTINUE
1028*
1029 ELSE IF( ipack.EQ.3 ) THEN
1030*
1031* 'C' -- Upper triangle packed Columnwise.
1032*
1033 icol = 1
1034 irow = 0
1035 DO 340 j = 1, m
1036 DO 330 i = 1, j
1037 irow = irow + 1
1038 IF( irow.GT.lda ) THEN
1039 irow = 1
1040 icol = icol + 1
1041 END IF
1042 a( irow, icol ) = a( i, j )
1043 330 CONTINUE
1044 340 CONTINUE
1045*
1046 ELSE IF( ipack.EQ.4 ) THEN
1047*
1048* 'R' -- Lower triangle packed Columnwise.
1049*
1050 icol = 1
1051 irow = 0
1052 DO 360 j = 1, m
1053 DO 350 i = j, m
1054 irow = irow + 1
1055 IF( irow.GT.lda ) THEN
1056 irow = 1
1057 icol = icol + 1
1058 END IF
1059 a( irow, icol ) = a( i, j )
1060 350 CONTINUE
1061 360 CONTINUE
1062*
1063 ELSE IF( ipack.GE.5 ) THEN
1064*
1065* 'B' -- The lower triangle is packed as a band matrix.
1066* 'Q' -- The upper triangle is packed as a band matrix.
1067* 'Z' -- The whole matrix is packed as a band matrix.
1068*
1069 IF( ipack.EQ.5 )
1070 $ uub = 0
1071 IF( ipack.EQ.6 )
1072 $ llb = 0
1073*
1074 DO 380 j = 1, uub
1075 DO 370 i = min( j+llb, m ), 1, -1
1076 a( i-j+uub+1, j ) = a( i, j )
1077 370 CONTINUE
1078 380 CONTINUE
1079*
1080 DO 400 j = uub + 2, n
1081 DO 390 i = j - uub, min( j+llb, m )
1082 a( i-j+uub+1, j ) = a( i, j )
1083 390 CONTINUE
1084 400 CONTINUE
1085 END IF
1086*
1087* If packed, zero out extraneous elements.
1088*
1089* Symmetric/Triangular Packed --
1090* zero out everything after A(IROW,ICOL)
1091*
1092 IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1093 DO 420 jc = icol, m
1094 DO 410 jr = irow + 1, lda
1095 a( jr, jc ) = zero
1096 410 CONTINUE
1097 irow = 0
1098 420 CONTINUE
1099*
1100 ELSE IF( ipack.GE.5 ) THEN
1101*
1102* Packed Band --
1103* 1st row is now in A( UUB+2-j, j), zero above it
1104* m-th row is now in A( M+UUB-j,j), zero below it
1105* last non-zero diagonal is now in A( UUB+LLB+1,j ),
1106* zero below it, too.
1107*
1108 ir1 = uub + llb + 2
1109 ir2 = uub + m + 2
1110 DO 450 jc = 1, n
1111 DO 430 jr = 1, uub + 1 - jc
1112 a( jr, jc ) = zero
1113 430 CONTINUE
1114 DO 440 jr = max( 1, min( ir1, ir2-jc ) ), lda
1115 a( jr, jc ) = zero
1116 440 CONTINUE
1117 450 CONTINUE
1118 END IF
1119 END IF
1120*
1121 RETURN
1122*
1123* End of DLATMS
1124*
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:113
dlagsy
subroutine dlagsy(n, k, d, a, lda, iseed, work, info)
DLAGSY
Definition dlagsy.f:101
dlagge
subroutine dlagge(m, n, kl, ku, d, a, lda, iseed, work, info)
DLAGGE
Definition dlagge.f:113
dlarot
subroutine dlarot(lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
DLAROT
Definition dlarot.f:226

◆ dlatmt()

subroutine dlatmt ( integer m,
integer n,
character dist,
integer, dimension( 4 ) iseed,
character sym,
double precision, dimension( * ) d,
integer mode,
double precision cond,
double precision dmax,
integer rank,
integer kl,
integer ku,
character pack,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) work,
integer info )

DLATMT

Purpose:
!>
!>    DLATMT generates random matrices with specified singular values
!>    (or symmetric/hermitian with specified eigenvalues)
!>    for testing LAPACK programs.
!>
!>    DLATMT operates by applying the following sequence of
!>    operations:
!>
!>      Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX, and SYM
!>         as described below.
!>
!>      Generate a matrix with the appropriate band structure, by one
!>         of two methods:
!>
!>      Method A:
!>          Generate a dense M x N matrix by multiplying D on the left
!>              and the right by random unitary matrices, then:
!>
!>          Reduce the bandwidth according to KL and KU, using
!>          Householder transformations.
!>
!>      Method B:
!>          Convert the bandwidth-0 (i.e., diagonal) matrix to a
!>              bandwidth-1 matrix using Givens rotations, 
!>              out-of-band elements back, much as in QR; then
!>              convert the bandwidth-1 to a bandwidth-2 matrix, etc.
!>              Note that for reasonably small bandwidths (relative to
!>              M and N) this requires less storage, as a dense matrix
!>              is not generated.  Also, for symmetric matrices, only
!>              one triangle is generated.
!>
!>      Method A is chosen if the bandwidth is a large fraction of the
!>          order of the matrix, and LDA is at least M (so a dense
!>          matrix can be stored.)  Method B is chosen if the bandwidth
!>          is small (< 1/2 N for symmetric, < .3 N+M for
!>          non-symmetric), or LDA is less than M and not less than the
!>          bandwidth.
!>
!>      Pack the matrix if desired. Options specified by PACK are:
!>         no packing
!>         zero out upper half (if symmetric)
!>         zero out lower half (if symmetric)
!>         store the upper half columnwise (if symmetric or upper
!>               triangular)
!>         store the lower half columnwise (if symmetric or lower
!>               triangular)
!>         store the lower triangle in banded format (if symmetric
!>               or lower triangular)
!>         store the upper triangle in banded format (if symmetric
!>               or upper triangular)
!>         store the entire matrix in banded format
!>      If Method B is chosen, and band format is specified, then the
!>         matrix will be generated in the band format, so no repacking
!>         will be necessary.
!>
Parameters
[in]M
!>          M is INTEGER
!>           The number of rows of A. Not modified.
!>
[in]N
!>          N is INTEGER
!>           The number of columns of A. Not modified.
!>
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate the random eigen-/singular values.
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!>
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to DLATMT
!>           to continue the same random number sequence.
!>           Changed on exit.
!>
[in]SYM
!>          SYM is CHARACTER*1
!>           If SYM='S' or 'H', the generated matrix is symmetric, with
!>             eigenvalues specified by D, COND, MODE, and DMAX; they
!>             may be positive, negative, or zero.
!>           If SYM='P', the generated matrix is symmetric, with
!>             eigenvalues (= singular values) specified by D, COND,
!>             MODE, and DMAX; they will not be negative.
!>           If SYM='N', the generated matrix is nonsymmetric, with
!>             singular values specified by D, COND, MODE, and DMAX;
!>             they will not be negative.
!>           Not modified.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension ( MIN( M , N ) )
!>           This array is used to specify the singular values or
!>           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
!>           assumed to contain the singular/eigenvalues, otherwise
!>           they will be computed according to MODE, COND, and DMAX,
!>           and placed in D.
!>           Modified if MODE is nonzero.
!>
[in]MODE
!>          MODE is INTEGER
!>           On entry this describes how the singular/eigenvalues are to
!>           be specified:
!>           MODE = 0 means use D as input
!>
!>           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
!>           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
!>
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
!>              the elements of D will also be multiplied by a random
!>              sign (i.e., +1 or -1.)
!>           Not modified.
!>
[in]COND
!>          COND is DOUBLE PRECISION
!>           On entry, this is used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!>
[in]DMAX
!>          DMAX is DOUBLE PRECISION
!>           If MODE is neither -6, 0 nor 6, the contents of D, as
!>           computed according to MODE and COND, will be scaled by
!>           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
!>           singular value (which is to say the norm) will be abs(DMAX).
!>           Note that DMAX need not be positive: if DMAX is negative
!>           (or zero), D will be scaled by a negative number (or zero).
!>           Not modified.
!>
[in]RANK
!>          RANK is INTEGER
!>           The rank of matrix to be generated for modes 1,2,3 only.
!>           D( RANK+1:N ) = 0.
!>           Not modified.
!>
[in]KL
!>          KL is INTEGER
!>           This specifies the lower bandwidth of the  matrix. For
!>           example, KL=0 implies upper triangular, KL=1 implies upper
!>           Hessenberg, and KL being at least M-1 means that the matrix
!>           has full lower bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!>
[in]KU
!>          KU is INTEGER
!>           This specifies the upper bandwidth of the  matrix. For
!>           example, KU=0 implies lower triangular, KU=1 implies lower
!>           Hessenberg, and KU being at least N-1 means that the matrix
!>           has full upper bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!>
[in]PACK
!>          PACK is CHARACTER*1
!>           This specifies packing of matrix as follows:
!>           'N' => no packing
!>           'U' => zero out all subdiagonal entries (if symmetric)
!>           'L' => zero out all superdiagonal entries (if symmetric)
!>           'C' => store the upper triangle columnwise
!>                  (only if the matrix is symmetric or upper triangular)
!>           'R' => store the lower triangle columnwise
!>                  (only if the matrix is symmetric or lower triangular)
!>           'B' => store the lower triangle in band storage scheme
!>                  (only if matrix symmetric or lower triangular)
!>           'Q' => store the upper triangle in band storage scheme
!>                  (only if matrix symmetric or upper triangular)
!>           'Z' => store the entire matrix in band storage scheme
!>                      (pivoting can be provided for by using this
!>                      option to store A in the trailing rows of
!>                      the allocated storage)
!>
!>           Using these options, the various LAPACK packed and banded
!>           storage schemes can be obtained:
!>           GB               - use 'Z'
!>           PB, SB or TB     - use 'B' or 'Q'
!>           PP, SP or TP     - use 'C' or 'R'
!>
!>           If two calls to DLATMT differ only in the PACK parameter,
!>           they will generate mathematically equivalent matrices.
!>           Not modified.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension ( LDA, N )
!>           On exit A is the desired test matrix.  A is first generated
!>           in full (unpacked) form, and then packed, if so specified
!>           by PACK.  Thus, the first M elements of the first N
!>           columns will always be modified.  If PACK specifies a
!>           packed or banded storage scheme, all LDA elements of the
!>           first N columns will be modified; the elements of the
!>           array which do not correspond to elements of the generated
!>           matrix are set to zero.
!>           Modified.
!>
[in]LDA
!>          LDA is INTEGER
!>           LDA specifies the first dimension of A as declared in the
!>           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
!>           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
!>           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
!>           If PACK='Z', LDA must be large enough to hold the packed
!>           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
!>           Not modified.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
!>           Workspace.
!>           Modified.
!>
[out]INFO
!>          INFO is INTEGER
!>           Error code.  On exit, INFO will be set to one of the
!>           following values:
!>             0 => normal return
!>            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
!>            -2 => N negative
!>            -3 => DIST illegal string
!>            -5 => SYM illegal string
!>            -7 => MODE not in range -6 to 6
!>            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>           -10 => KL negative
!>           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
!>           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
!>                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
!>                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
!>                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
!>                  N.
!>           -14 => LDA is less than M, or PACK='Z' and LDA is less than
!>                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
!>            1  => Error return from DLATM7
!>            2  => Cannot scale to DMAX (max. sing. value is 0)
!>            3  => Error return from DLAGGE or DLAGSY
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 329 of file dlatmt.f.

331*
332* -- LAPACK computational routine --
333* -- LAPACK is a software package provided by Univ. of Tennessee, --
334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335*
336* .. Scalar Arguments ..
337 DOUBLE PRECISION COND, DMAX
338 INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK
339 CHARACTER DIST, PACK, SYM
340* ..
341* .. Array Arguments ..
342 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
343 INTEGER ISEED( 4 )
344* ..
345*
346* =====================================================================
347*
348* .. Parameters ..
349 DOUBLE PRECISION ZERO
350 parameter( zero = 0.0d0 )
351 DOUBLE PRECISION ONE
352 parameter( one = 1.0d0 )
353 DOUBLE PRECISION TWOPI
354 parameter( twopi = 6.28318530717958647692528676655900576839d+0 )
355* ..
356* .. Local Scalars ..
357 DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
358 INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
359 $ IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
360 $ IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
361 $ JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
362 $ UUB
363 LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN
364* ..
365* .. External Functions ..
366 DOUBLE PRECISION DLARND
367 LOGICAL LSAME
368 EXTERNAL dlarnd, lsame
369* ..
370* .. External Subroutines ..
371 EXTERNAL dlatm7, dcopy, dlagge, dlagsy, dlarot,
372 $ dlartg, dlaset, dscal, xerbla
373* ..
374* .. Intrinsic Functions ..
375 INTRINSIC abs, cos, dble, max, min, mod, sin
376* ..
377* .. Executable Statements ..
378*
379* 1) Decode and Test the input parameters.
380* Initialize flags & seed.
381*
382 info = 0
383*
384* Quick return if possible
385*
386 IF( m.EQ.0 .OR. n.EQ.0 )
387 $ RETURN
388*
389* Decode DIST
390*
391 IF( lsame( dist, 'U' ) ) THEN
392 idist = 1
393 ELSE IF( lsame( dist, 'S' ) ) THEN
394 idist = 2
395 ELSE IF( lsame( dist, 'N' ) ) THEN
396 idist = 3
397 ELSE
398 idist = -1
399 END IF
400*
401* Decode SYM
402*
403 IF( lsame( sym, 'N' ) ) THEN
404 isym = 1
405 irsign = 0
406 ELSE IF( lsame( sym, 'P' ) ) THEN
407 isym = 2
408 irsign = 0
409 ELSE IF( lsame( sym, 'S' ) ) THEN
410 isym = 2
411 irsign = 1
412 ELSE IF( lsame( sym, 'H' ) ) THEN
413 isym = 2
414 irsign = 1
415 ELSE
416 isym = -1
417 END IF
418*
419* Decode PACK
420*
421 isympk = 0
422 IF( lsame( pack, 'N' ) ) THEN
423 ipack = 0
424 ELSE IF( lsame( pack, 'U' ) ) THEN
425 ipack = 1
426 isympk = 1
427 ELSE IF( lsame( pack, 'L' ) ) THEN
428 ipack = 2
429 isympk = 1
430 ELSE IF( lsame( pack, 'C' ) ) THEN
431 ipack = 3
432 isympk = 2
433 ELSE IF( lsame( pack, 'R' ) ) THEN
434 ipack = 4
435 isympk = 3
436 ELSE IF( lsame( pack, 'B' ) ) THEN
437 ipack = 5
438 isympk = 3
439 ELSE IF( lsame( pack, 'Q' ) ) THEN
440 ipack = 6
441 isympk = 2
442 ELSE IF( lsame( pack, 'Z' ) ) THEN
443 ipack = 7
444 ELSE
445 ipack = -1
446 END IF
447*
448* Set certain internal parameters
449*
450 mnmin = min( m, n )
451 llb = min( kl, m-1 )
452 uub = min( ku, n-1 )
453 mr = min( m, n+llb )
454 nc = min( n, m+uub )
455*
456 IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
457 minlda = uub + 1
458 ELSE IF( ipack.EQ.7 ) THEN
459 minlda = llb + uub + 1
460 ELSE
461 minlda = m
462 END IF
463*
464* Use Givens rotation method if bandwidth small enough,
465* or if LDA is too small to store the matrix unpacked.
466*
467 givens = .false.
468 IF( isym.EQ.1 ) THEN
469 IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
470 $ givens = .true.
471 ELSE
472 IF( 2*llb.LT.m )
473 $ givens = .true.
474 END IF
475 IF( lda.LT.m .AND. lda.GE.minlda )
476 $ givens = .true.
477*
478* Set INFO if an error
479*
480 IF( m.LT.0 ) THEN
481 info = -1
482 ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
483 info = -1
484 ELSE IF( n.LT.0 ) THEN
485 info = -2
486 ELSE IF( idist.EQ.-1 ) THEN
487 info = -3
488 ELSE IF( isym.EQ.-1 ) THEN
489 info = -5
490 ELSE IF( abs( mode ).GT.6 ) THEN
491 info = -7
492 ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
493 $ THEN
494 info = -8
495 ELSE IF( kl.LT.0 ) THEN
496 info = -10
497 ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
498 info = -11
499 ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
500 $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
501 $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
502 $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
503 info = -12
504 ELSE IF( lda.LT.max( 1, minlda ) ) THEN
505 info = -14
506 END IF
507*
508 IF( info.NE.0 ) THEN
509 CALL xerbla( 'DLATMT', -info )
510 RETURN
511 END IF
512*
513* Initialize random number generator
514*
515 DO 100 i = 1, 4
516 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
517 100 CONTINUE
518*
519 IF( mod( iseed( 4 ), 2 ).NE.1 )
520 $ iseed( 4 ) = iseed( 4 ) + 1
521*
522* 2) Set up D if indicated.
523*
524* Compute D according to COND and MODE
525*
526 CALL dlatm7( mode, cond, irsign, idist, iseed, d, mnmin, rank,
527 $ iinfo )
528 IF( iinfo.NE.0 ) THEN
529 info = 1
530 RETURN
531 END IF
532*
533* Choose Top-Down if D is (apparently) increasing,
534* Bottom-Up if D is (apparently) decreasing.
535*
536 IF( abs( d( 1 ) ).LE.abs( d( rank ) ) ) THEN
537 topdwn = .true.
538 ELSE
539 topdwn = .false.
540 END IF
541*
542 IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
543*
544* Scale by DMAX
545*
546 temp = abs( d( 1 ) )
547 DO 110 i = 2, rank
548 temp = max( temp, abs( d( i ) ) )
549 110 CONTINUE
550*
551 IF( temp.GT.zero ) THEN
552 alpha = dmax / temp
553 ELSE
554 info = 2
555 RETURN
556 END IF
557*
558 CALL dscal( rank, alpha, d, 1 )
559*
560 END IF
561*
562* 3) Generate Banded Matrix using Givens rotations.
563* Also the special case of UUB=LLB=0
564*
565* Compute Addressing constants to cover all
566* storage formats. Whether GE, SY, GB, or SB,
567* upper or lower triangle or both,
568* the (i,j)-th element is in
569* A( i - ISKEW*j + IOFFST, j )
570*
571 IF( ipack.GT.4 ) THEN
572 ilda = lda - 1
573 iskew = 1
574 IF( ipack.GT.5 ) THEN
575 ioffst = uub + 1
576 ELSE
577 ioffst = 1
578 END IF
579 ELSE
580 ilda = lda
581 iskew = 0
582 ioffst = 0
583 END IF
584*
585* IPACKG is the format that the matrix is generated in. If this is
586* different from IPACK, then the matrix must be repacked at the
587* end. It also signals how to compute the norm, for scaling.
588*
589 ipackg = 0
590 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
591*
592* Diagonal Matrix -- We are done, unless it
593* is to be stored SP/PP/TP (PACK='R' or 'C')
594*
595 IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
596 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
597 IF( ipack.LE.2 .OR. ipack.GE.5 )
598 $ ipackg = ipack
599*
600 ELSE IF( givens ) THEN
601*
602* Check whether to use Givens rotations,
603* Householder transformations, or nothing.
604*
605 IF( isym.EQ.1 ) THEN
606*
607* Non-symmetric -- A = U D V
608*
609 IF( ipack.GT.4 ) THEN
610 ipackg = ipack
611 ELSE
612 ipackg = 0
613 END IF
614*
615 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
616*
617 IF( topdwn ) THEN
618 jkl = 0
619 DO 140 jku = 1, uub
620*
621* Transform from bandwidth JKL, JKU-1 to JKL, JKU
622*
623* Last row actually rotated is M
624* Last column actually rotated is MIN( M+JKU, N )
625*
626 DO 130 jr = 1, min( m+jku, n ) + jkl - 1
627 extra = zero
628 angle = twopi*dlarnd( 1, iseed )
629 c = cos( angle )
630 s = sin( angle )
631 icol = max( 1, jr-jkl )
632 IF( jr.LT.m ) THEN
633 il = min( n, jr+jku ) + 1 - icol
634 CALL dlarot( .true., jr.GT.jkl, .false., il, c,
635 $ s, a( jr-iskew*icol+ioffst, icol ),
636 $ ilda, extra, dummy )
637 END IF
638*
639* Chase "EXTRA" back up
640*
641 ir = jr
642 ic = icol
643 DO 120 jch = jr - jkl, 1, -jkl - jku
644 IF( ir.LT.m ) THEN
645 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
646 $ ic+1 ), extra, c, s, dummy )
647 END IF
648 irow = max( 1, jch-jku )
649 il = ir + 2 - irow
650 temp = zero
651 iltemp = jch.GT.jku
652 CALL dlarot( .false., iltemp, .true., il, c, -s,
653 $ a( irow-iskew*ic+ioffst, ic ),
654 $ ilda, temp, extra )
655 IF( iltemp ) THEN
656 CALL dlartg( a( irow+1-iskew*( ic+1 )+ioffst,
657 $ ic+1 ), temp, c, s, dummy )
658 icol = max( 1, jch-jku-jkl )
659 il = ic + 2 - icol
660 extra = zero
661 CALL dlarot( .true., jch.GT.jku+jkl, .true.,
662 $ il, c, -s, a( irow-iskew*icol+
663 $ ioffst, icol ), ilda, extra,
664 $ temp )
665 ic = icol
666 ir = irow
667 END IF
668 120 CONTINUE
669 130 CONTINUE
670 140 CONTINUE
671*
672 jku = uub
673 DO 170 jkl = 1, llb
674*
675* Transform from bandwidth JKL-1, JKU to JKL, JKU
676*
677 DO 160 jc = 1, min( n+jkl, m ) + jku - 1
678 extra = zero
679 angle = twopi*dlarnd( 1, iseed )
680 c = cos( angle )
681 s = sin( angle )
682 irow = max( 1, jc-jku )
683 IF( jc.LT.n ) THEN
684 il = min( m, jc+jkl ) + 1 - irow
685 CALL dlarot( .false., jc.GT.jku, .false., il, c,
686 $ s, a( irow-iskew*jc+ioffst, jc ),
687 $ ilda, extra, dummy )
688 END IF
689*
690* Chase "EXTRA" back up
691*
692 ic = jc
693 ir = irow
694 DO 150 jch = jc - jku, 1, -jkl - jku
695 IF( ic.LT.n ) THEN
696 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
697 $ ic+1 ), extra, c, s, dummy )
698 END IF
699 icol = max( 1, jch-jkl )
700 il = ic + 2 - icol
701 temp = zero
702 iltemp = jch.GT.jkl
703 CALL dlarot( .true., iltemp, .true., il, c, -s,
704 $ a( ir-iskew*icol+ioffst, icol ),
705 $ ilda, temp, extra )
706 IF( iltemp ) THEN
707 CALL dlartg( a( ir+1-iskew*( icol+1 )+ioffst,
708 $ icol+1 ), temp, c, s, dummy )
709 irow = max( 1, jch-jkl-jku )
710 il = ir + 2 - irow
711 extra = zero
712 CALL dlarot( .false., jch.GT.jkl+jku, .true.,
713 $ il, c, -s, a( irow-iskew*icol+
714 $ ioffst, icol ), ilda, extra,
715 $ temp )
716 ic = icol
717 ir = irow
718 END IF
719 150 CONTINUE
720 160 CONTINUE
721 170 CONTINUE
722*
723 ELSE
724*
725* Bottom-Up -- Start at the bottom right.
726*
727 jkl = 0
728 DO 200 jku = 1, uub
729*
730* Transform from bandwidth JKL, JKU-1 to JKL, JKU
731*
732* First row actually rotated is M
733* First column actually rotated is MIN( M+JKU, N )
734*
735 iendch = min( m, n+jkl ) - 1
736 DO 190 jc = min( m+jku, n ) - 1, 1 - jkl, -1
737 extra = zero
738 angle = twopi*dlarnd( 1, iseed )
739 c = cos( angle )
740 s = sin( angle )
741 irow = max( 1, jc-jku+1 )
742 IF( jc.GT.0 ) THEN
743 il = min( m, jc+jkl+1 ) + 1 - irow
744 CALL dlarot( .false., .false., jc+jkl.LT.m, il,
745 $ c, s, a( irow-iskew*jc+ioffst,
746 $ jc ), ilda, dummy, extra )
747 END IF
748*
749* Chase "EXTRA" back down
750*
751 ic = jc
752 DO 180 jch = jc + jkl, iendch, jkl + jku
753 ilextr = ic.GT.0
754 IF( ilextr ) THEN
755 CALL dlartg( a( jch-iskew*ic+ioffst, ic ),
756 $ extra, c, s, dummy )
757 END IF
758 ic = max( 1, ic )
759 icol = min( n-1, jch+jku )
760 iltemp = jch + jku.LT.n
761 temp = zero
762 CALL dlarot( .true., ilextr, iltemp, icol+2-ic,
763 $ c, s, a( jch-iskew*ic+ioffst, ic ),
764 $ ilda, extra, temp )
765 IF( iltemp ) THEN
766 CALL dlartg( a( jch-iskew*icol+ioffst,
767 $ icol ), temp, c, s, dummy )
768 il = min( iendch, jch+jkl+jku ) + 2 - jch
769 extra = zero
770 CALL dlarot( .false., .true.,
771 $ jch+jkl+jku.LE.iendch, il, c, s,
772 $ a( jch-iskew*icol+ioffst,
773 $ icol ), ilda, temp, extra )
774 ic = icol
775 END IF
776 180 CONTINUE
777 190 CONTINUE
778 200 CONTINUE
779*
780 jku = uub
781 DO 230 jkl = 1, llb
782*
783* Transform from bandwidth JKL-1, JKU to JKL, JKU
784*
785* First row actually rotated is MIN( N+JKL, M )
786* First column actually rotated is N
787*
788 iendch = min( n, m+jku ) - 1
789 DO 220 jr = min( n+jkl, m ) - 1, 1 - jku, -1
790 extra = zero
791 angle = twopi*dlarnd( 1, iseed )
792 c = cos( angle )
793 s = sin( angle )
794 icol = max( 1, jr-jkl+1 )
795 IF( jr.GT.0 ) THEN
796 il = min( n, jr+jku+1 ) + 1 - icol
797 CALL dlarot( .true., .false., jr+jku.LT.n, il,
798 $ c, s, a( jr-iskew*icol+ioffst,
799 $ icol ), ilda, dummy, extra )
800 END IF
801*
802* Chase "EXTRA" back down
803*
804 ir = jr
805 DO 210 jch = jr + jku, iendch, jkl + jku
806 ilextr = ir.GT.0
807 IF( ilextr ) THEN
808 CALL dlartg( a( ir-iskew*jch+ioffst, jch ),
809 $ extra, c, s, dummy )
810 END IF
811 ir = max( 1, ir )
812 irow = min( m-1, jch+jkl )
813 iltemp = jch + jkl.LT.m
814 temp = zero
815 CALL dlarot( .false., ilextr, iltemp, irow+2-ir,
816 $ c, s, a( ir-iskew*jch+ioffst,
817 $ jch ), ilda, extra, temp )
818 IF( iltemp ) THEN
819 CALL dlartg( a( irow-iskew*jch+ioffst, jch ),
820 $ temp, c, s, dummy )
821 il = min( iendch, jch+jkl+jku ) + 2 - jch
822 extra = zero
823 CALL dlarot( .true., .true.,
824 $ jch+jkl+jku.LE.iendch, il, c, s,
825 $ a( irow-iskew*jch+ioffst, jch ),
826 $ ilda, temp, extra )
827 ir = irow
828 END IF
829 210 CONTINUE
830 220 CONTINUE
831 230 CONTINUE
832 END IF
833*
834 ELSE
835*
836* Symmetric -- A = U D U'
837*
838 ipackg = ipack
839 ioffg = ioffst
840*
841 IF( topdwn ) THEN
842*
843* Top-Down -- Generate Upper triangle only
844*
845 IF( ipack.GE.5 ) THEN
846 ipackg = 6
847 ioffg = uub + 1
848 ELSE
849 ipackg = 1
850 END IF
851 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
852*
853 DO 260 k = 1, uub
854 DO 250 jc = 1, n - 1
855 irow = max( 1, jc-k )
856 il = min( jc+1, k+2 )
857 extra = zero
858 temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
859 angle = twopi*dlarnd( 1, iseed )
860 c = cos( angle )
861 s = sin( angle )
862 CALL dlarot( .false., jc.GT.k, .true., il, c, s,
863 $ a( irow-iskew*jc+ioffg, jc ), ilda,
864 $ extra, temp )
865 CALL dlarot( .true., .true., .false.,
866 $ min( k, n-jc )+1, c, s,
867 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
868 $ temp, dummy )
869*
870* Chase EXTRA back up the matrix
871*
872 icol = jc
873 DO 240 jch = jc - k, 1, -k
874 CALL dlartg( a( jch+1-iskew*( icol+1 )+ioffg,
875 $ icol+1 ), extra, c, s, dummy )
876 temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
877 CALL dlarot( .true., .true., .true., k+2, c, -s,
878 $ a( ( 1-iskew )*jch+ioffg, jch ),
879 $ ilda, temp, extra )
880 irow = max( 1, jch-k )
881 il = min( jch+1, k+2 )
882 extra = zero
883 CALL dlarot( .false., jch.GT.k, .true., il, c,
884 $ -s, a( irow-iskew*jch+ioffg, jch ),
885 $ ilda, extra, temp )
886 icol = jch
887 240 CONTINUE
888 250 CONTINUE
889 260 CONTINUE
890*
891* If we need lower triangle, copy from upper. Note that
892* the order of copying is chosen to work for 'q' -> 'b'
893*
894 IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
895 DO 280 jc = 1, n
896 irow = ioffst - iskew*jc
897 DO 270 jr = jc, min( n, jc+uub )
898 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
899 270 CONTINUE
900 280 CONTINUE
901 IF( ipack.EQ.5 ) THEN
902 DO 300 jc = n - uub + 1, n
903 DO 290 jr = n + 2 - jc, uub + 1
904 a( jr, jc ) = zero
905 290 CONTINUE
906 300 CONTINUE
907 END IF
908 IF( ipackg.EQ.6 ) THEN
909 ipackg = ipack
910 ELSE
911 ipackg = 0
912 END IF
913 END IF
914 ELSE
915*
916* Bottom-Up -- Generate Lower triangle only
917*
918 IF( ipack.GE.5 ) THEN
919 ipackg = 5
920 IF( ipack.EQ.6 )
921 $ ioffg = 1
922 ELSE
923 ipackg = 2
924 END IF
925 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
926*
927 DO 330 k = 1, uub
928 DO 320 jc = n - 1, 1, -1
929 il = min( n+1-jc, k+2 )
930 extra = zero
931 temp = a( 1+( 1-iskew )*jc+ioffg, jc )
932 angle = twopi*dlarnd( 1, iseed )
933 c = cos( angle )
934 s = -sin( angle )
935 CALL dlarot( .false., .true., n-jc.GT.k, il, c, s,
936 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
937 $ temp, extra )
938 icol = max( 1, jc-k+1 )
939 CALL dlarot( .true., .false., .true., jc+2-icol, c,
940 $ s, a( jc-iskew*icol+ioffg, icol ),
941 $ ilda, dummy, temp )
942*
943* Chase EXTRA back down the matrix
944*
945 icol = jc
946 DO 310 jch = jc + k, n - 1, k
947 CALL dlartg( a( jch-iskew*icol+ioffg, icol ),
948 $ extra, c, s, dummy )
949 temp = a( 1+( 1-iskew )*jch+ioffg, jch )
950 CALL dlarot( .true., .true., .true., k+2, c, s,
951 $ a( jch-iskew*icol+ioffg, icol ),
952 $ ilda, extra, temp )
953 il = min( n+1-jch, k+2 )
954 extra = zero
955 CALL dlarot( .false., .true., n-jch.GT.k, il, c,
956 $ s, a( ( 1-iskew )*jch+ioffg, jch ),
957 $ ilda, temp, extra )
958 icol = jch
959 310 CONTINUE
960 320 CONTINUE
961 330 CONTINUE
962*
963* If we need upper triangle, copy from lower. Note that
964* the order of copying is chosen to work for 'b' -> 'q'
965*
966 IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
967 DO 350 jc = n, 1, -1
968 irow = ioffst - iskew*jc
969 DO 340 jr = jc, max( 1, jc-uub ), -1
970 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
971 340 CONTINUE
972 350 CONTINUE
973 IF( ipack.EQ.6 ) THEN
974 DO 370 jc = 1, uub
975 DO 360 jr = 1, uub + 1 - jc
976 a( jr, jc ) = zero
977 360 CONTINUE
978 370 CONTINUE
979 END IF
980 IF( ipackg.EQ.5 ) THEN
981 ipackg = ipack
982 ELSE
983 ipackg = 0
984 END IF
985 END IF
986 END IF
987 END IF
988*
989 ELSE
990*
991* 4) Generate Banded Matrix by first
992* Rotating by random Unitary matrices,
993* then reducing the bandwidth using Householder
994* transformations.
995*
996* Note: we should get here only if LDA .ge. N
997*
998 IF( isym.EQ.1 ) THEN
999*
1000* Non-symmetric -- A = U D V
1001*
1002 CALL dlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1003 $ iinfo )
1004 ELSE
1005*
1006* Symmetric -- A = U D U'
1007*
1008 CALL dlagsy( m, llb, d, a, lda, iseed, work, iinfo )
1009*
1010 END IF
1011 IF( iinfo.NE.0 ) THEN
1012 info = 3
1013 RETURN
1014 END IF
1015 END IF
1016*
1017* 5) Pack the matrix
1018*
1019 IF( ipack.NE.ipackg ) THEN
1020 IF( ipack.EQ.1 ) THEN
1021*
1022* 'U' -- Upper triangular, not packed
1023*
1024 DO 390 j = 1, m
1025 DO 380 i = j + 1, m
1026 a( i, j ) = zero
1027 380 CONTINUE
1028 390 CONTINUE
1029*
1030 ELSE IF( ipack.EQ.2 ) THEN
1031*
1032* 'L' -- Lower triangular, not packed
1033*
1034 DO 410 j = 2, m
1035 DO 400 i = 1, j - 1
1036 a( i, j ) = zero
1037 400 CONTINUE
1038 410 CONTINUE
1039*
1040 ELSE IF( ipack.EQ.3 ) THEN
1041*
1042* 'C' -- Upper triangle packed Columnwise.
1043*
1044 icol = 1
1045 irow = 0
1046 DO 430 j = 1, m
1047 DO 420 i = 1, j
1048 irow = irow + 1
1049 IF( irow.GT.lda ) THEN
1050 irow = 1
1051 icol = icol + 1
1052 END IF
1053 a( irow, icol ) = a( i, j )
1054 420 CONTINUE
1055 430 CONTINUE
1056*
1057 ELSE IF( ipack.EQ.4 ) THEN
1058*
1059* 'R' -- Lower triangle packed Columnwise.
1060*
1061 icol = 1
1062 irow = 0
1063 DO 450 j = 1, m
1064 DO 440 i = j, m
1065 irow = irow + 1
1066 IF( irow.GT.lda ) THEN
1067 irow = 1
1068 icol = icol + 1
1069 END IF
1070 a( irow, icol ) = a( i, j )
1071 440 CONTINUE
1072 450 CONTINUE
1073*
1074 ELSE IF( ipack.GE.5 ) THEN
1075*
1076* 'B' -- The lower triangle is packed as a band matrix.
1077* 'Q' -- The upper triangle is packed as a band matrix.
1078* 'Z' -- The whole matrix is packed as a band matrix.
1079*
1080 IF( ipack.EQ.5 )
1081 $ uub = 0
1082 IF( ipack.EQ.6 )
1083 $ llb = 0
1084*
1085 DO 470 j = 1, uub
1086 DO 460 i = min( j+llb, m ), 1, -1
1087 a( i-j+uub+1, j ) = a( i, j )
1088 460 CONTINUE
1089 470 CONTINUE
1090*
1091 DO 490 j = uub + 2, n
1092 DO 480 i = j - uub, min( j+llb, m )
1093 a( i-j+uub+1, j ) = a( i, j )
1094 480 CONTINUE
1095 490 CONTINUE
1096 END IF
1097*
1098* If packed, zero out extraneous elements.
1099*
1100* Symmetric/Triangular Packed --
1101* zero out everything after A(IROW,ICOL)
1102*
1103 IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1104 DO 510 jc = icol, m
1105 DO 500 jr = irow + 1, lda
1106 a( jr, jc ) = zero
1107 500 CONTINUE
1108 irow = 0
1109 510 CONTINUE
1110*
1111 ELSE IF( ipack.GE.5 ) THEN
1112*
1113* Packed Band --
1114* 1st row is now in A( UUB+2-j, j), zero above it
1115* m-th row is now in A( M+UUB-j,j), zero below it
1116* last non-zero diagonal is now in A( UUB+LLB+1,j ),
1117* zero below it, too.
1118*
1119 ir1 = uub + llb + 2
1120 ir2 = uub + m + 2
1121 DO 540 jc = 1, n
1122 DO 520 jr = 1, uub + 1 - jc
1123 a( jr, jc ) = zero
1124 520 CONTINUE
1125 DO 530 jr = max( 1, min( ir1, ir2-jc ) ), lda
1126 a( jr, jc ) = zero
1127 530 CONTINUE
1128 540 CONTINUE
1129 END IF
1130 END IF
1131*
1132 RETURN
1133*
1134* End of DLATMT
1135*
dlatm7
subroutine dlatm7(mode, cond, irsign, idist, iseed, d, n, rank, info)
DLATM7
Definition dlatm7.f:122