LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for double:

Functions

subroutine dlagge (M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
 DLAGGE More...
 
subroutine dlagsy (N, K, D, A, LDA, ISEED, WORK, INFO)
 DLAGSY More...
 
subroutine dlahilb (N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
 DLAHILB More...
 
subroutine dlakf2 (M, N, A, LDA, B, D, E, Z, LDZ)
 DLAKF2 More...
 
subroutine dlarge (N, A, LDA, ISEED, WORK, INFO)
 DLARGE More...
 
double precision function dlarnd (IDIST, ISEED)
 DLARND More...
 
subroutine dlaror (SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
 DLAROR More...
 
subroutine dlarot (LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
 DLAROT More...
 
subroutine dlatm1 (MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
 DLATM1 More...
 
double precision function dlatm2 (M, N, I, J, KL, KU, IDIST, ISEED, D, IGRADE, DL, DR, IPVTNG, IWORK, SPARSE)
 DLATM2 More...
 
double precision function dlatm3 (M, N, I, J, ISUB, JSUB, KL, KU, IDIST, ISEED, D, IGRADE, DL, DR, IPVTNG, IWORK, SPARSE)
 DLATM3 More...
 
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 More...
 
subroutine dlatm6 (TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
 DLATM6 More...
 
subroutine dlatm7 (MODE, COND, IRSIGN, IDIST, ISEED, D, N, RANK, INFO)
 DLATM7 More...
 
subroutine dlatme (N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
 DLATME More...
 
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 More...
 
subroutine dlatms (M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
 DLATMS More...
 
subroutine dlatmt (M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RANK, KL, KU, PACK, A, LDA, WORK, INFO)
 DLATMT More...
 

Detailed Description

This is the group of double LAPACK TESTING MATGEN routines.

Function Documentation

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.
Date
November 2015

Definition at line 115 of file dlagge.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 103 of file dlagsy.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2015

Definition at line 126 of file dlahilb.f.

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

Here is the call graph for this function:

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.
Date
November 2011

Definition at line 107 of file dlakf2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 89 of file dlarge.f.

89 *
90 * -- LAPACK auxiliary routine (version 3.4.0) --
91 * -- LAPACK is a software package provided by Univ. of Tennessee, --
92 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
93 * November 2011
94 *
95 * .. Scalar Arguments ..
96  INTEGER info, lda, n
97 * ..
98 * .. Array Arguments ..
99  INTEGER iseed( 4 )
100  DOUBLE PRECISION a( lda, * ), work( * )
101 * ..
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106  DOUBLE PRECISION zero, one
107  parameter( zero = 0.0d+0, one = 1.0d+0 )
108 * ..
109 * .. Local Scalars ..
110  INTEGER i
111  DOUBLE PRECISION tau, wa, wb, wn
112 * ..
113 * .. External Subroutines ..
114  EXTERNAL dgemv, dger, dlarnv, dscal, xerbla
115 * ..
116 * .. Intrinsic Functions ..
117  INTRINSIC max, sign
118 * ..
119 * .. External Functions ..
120  DOUBLE PRECISION dnrm2
121  EXTERNAL dnrm2
122 * ..
123 * .. Executable Statements ..
124 *
125 * Test the input arguments
126 *
127  info = 0
128  IF( n.LT.0 ) THEN
129  info = -1
130  ELSE IF( lda.LT.max( 1, n ) ) THEN
131  info = -3
132  END IF
133  IF( info.LT.0 ) THEN
134  CALL xerbla( 'DLARGE', -info )
135  RETURN
136  END IF
137 *
138 * pre- and post-multiply A by random orthogonal matrix
139 *
140  DO 10 i = n, 1, -1
141 *
142 * generate random reflection
143 *
144  CALL dlarnv( 3, iseed, n-i+1, work )
145  wn = dnrm2( n-i+1, work, 1 )
146  wa = sign( wn, work( 1 ) )
147  IF( wn.EQ.zero ) THEN
148  tau = zero
149  ELSE
150  wb = work( 1 ) + wa
151  CALL dscal( n-i, one / wb, work( 2 ), 1 )
152  work( 1 ) = one
153  tau = wb / wa
154  END IF
155 *
156 * multiply A(i:n,1:n) by random reflection from the left
157 *
158  CALL dgemv( 'Transpose', n-i+1, n, one, a( i, 1 ), lda, work,
159  $ 1, zero, work( n+1 ), 1 )
160  CALL dger( n-i+1, n, -tau, work, 1, work( n+1 ), 1, a( i, 1 ),
161  $ lda )
162 *
163 * multiply A(1:n,i:n) by random reflection from the right
164 *
165  CALL dgemv( 'No transpose', n, n-i+1, one, a( 1, i ), lda,
166  $ work, 1, zero, work( n+1 ), 1 )
167  CALL dger( n, n-i+1, -tau, work( n+1 ), 1, work, 1, a( 1, i ),
168  $ lda )
169  10 CONTINUE
170  RETURN
171 *
172 * End of DLARGE
173 *
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011
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 75 of file dlarnd.f.

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

Definition at line 148 of file dlaror.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "p" and "q":

    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 "effective" 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.
Date
November 2011

Definition at line 228 of file dlarot.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2015

Definition at line 137 of file dlatm1.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 paramters. 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 matix 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.
Date
November 2011

Definition at line 210 of file dlatm2.f.

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

Definition at line 228 of file dlatm3.f.

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

Here is the caller graph for this function:

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
          "Points" to a certian type of the matrices to generate
          (see futher 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.
Date
November 2011
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 270 of file dlatm5.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 futher 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.
Date
November 2011

Definition at line 178 of file dlatm6.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
November 2011

Definition at line 124 of file dlatm7.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "similarity
           transform", 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.
Date
November 2011

Definition at line 334 of file dlatme.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 "2-cycles", 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]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]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]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, N-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.
Date
November 2011

Definition at line 473 of file dlatmr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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, "chasing"
              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.
Date
November 2011

Definition at line 323 of file dlatms.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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, "chasing"
              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.
Date
November 2011

Definition at line 333 of file dlatmt.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: