|
kusano |
7d535a |
/* -- translated by f2c (version 19940927).
|
|
kusano |
7d535a |
You must link the resulting object file with the libraries:
|
|
kusano |
7d535a |
-lf2c -lm (in that order)
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
#include "f2c.h"
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Table of constant values */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
static integer c__1 = 1;
|
|
kusano |
7d535a |
static real c_b23 = 0.f;
|
|
kusano |
7d535a |
static integer c__0 = 0;
|
|
kusano |
7d535a |
static real c_b39 = 1.f;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Subroutine */ int slatme_(integer *n, char *dist, integer *iseed, real *d,
|
|
kusano |
7d535a |
integer *mode, real *cond, real *dmax__, char *ei, char *rsign, char *
|
|
kusano |
7d535a |
upper, char *sim, real *ds, integer *modes, real *conds, integer *kl,
|
|
kusano |
7d535a |
integer *ku, real *anorm, real *a, integer *lda, real *work, integer *
|
|
kusano |
7d535a |
info)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
/* System generated locals */
|
|
kusano |
7d535a |
integer a_dim1, a_offset, i__1, i__2;
|
|
kusano |
7d535a |
real r__1, r__2, r__3;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Local variables */
|
|
kusano |
7d535a |
static logical bads;
|
|
kusano |
7d535a |
extern /* Subroutine */ int sger_(integer *, integer *, real *, real *,
|
|
kusano |
7d535a |
integer *, real *, integer *, real *, integer *);
|
|
kusano |
7d535a |
static integer isim;
|
|
kusano |
7d535a |
static real temp;
|
|
kusano |
7d535a |
static logical badei;
|
|
kusano |
7d535a |
static integer i, j;
|
|
kusano |
7d535a |
static real alpha;
|
|
kusano |
7d535a |
extern logical lsame_(char *, char *);
|
|
kusano |
7d535a |
static integer iinfo;
|
|
kusano |
7d535a |
extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
|
|
kusano |
7d535a |
static real tempa[1];
|
|
kusano |
7d535a |
static integer icols;
|
|
kusano |
7d535a |
static logical useei;
|
|
kusano |
7d535a |
static integer idist;
|
|
kusano |
7d535a |
extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *,
|
|
kusano |
7d535a |
real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *);
|
|
kusano |
7d535a |
static integer irows;
|
|
kusano |
7d535a |
extern /* Subroutine */ int slatm1_(integer *, real *, integer *, integer
|
|
kusano |
7d535a |
*, integer *, real *, integer *, integer *);
|
|
kusano |
7d535a |
static integer ic, jc, ir, jr;
|
|
kusano |
7d535a |
extern doublereal slange_(char *, integer *, integer *, real *, integer *,
|
|
kusano |
7d535a |
real *);
|
|
kusano |
7d535a |
extern /* Subroutine */ int slarge_(integer *, real *, integer *, integer
|
|
kusano |
7d535a |
*, real *, integer *), slarfg_(integer *, real *, real *, integer
|
|
kusano |
7d535a |
*, real *), xerbla_(char *, integer *);
|
|
kusano |
7d535a |
extern doublereal slaran_(integer *);
|
|
kusano |
7d535a |
static integer irsign;
|
|
kusano |
7d535a |
extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *,
|
|
kusano |
7d535a |
real *, real *, integer *);
|
|
kusano |
7d535a |
static integer iupper;
|
|
kusano |
7d535a |
extern /* Subroutine */ int slarnv_(integer *, integer *, integer *, real
|
|
kusano |
7d535a |
*);
|
|
kusano |
7d535a |
static real xnorms;
|
|
kusano |
7d535a |
static integer jcr;
|
|
kusano |
7d535a |
static real tau;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* -- LAPACK test routine (version 2.0) --
|
|
kusano |
7d535a |
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
kusano |
7d535a |
Courant Institute, Argonne National Lab, and Rice University
|
|
kusano |
7d535a |
September 30, 1994
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Purpose
|
|
kusano |
7d535a |
=======
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SLATME generates random non-symmetric square matrices with
|
|
kusano |
7d535a |
specified eigenvalues for testing LAPACK programs.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SLATME operates by applying the following sequence of
|
|
kusano |
7d535a |
operations:
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
1. Set the diagonal to D, where D may be input or
|
|
kusano |
7d535a |
computed according to MODE, COND, DMAX, and RSIGN
|
|
kusano |
7d535a |
as described below.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
|
|
kusano |
7d535a |
or MODE=5), certain pairs of adjacent elements of D are
|
|
kusano |
7d535a |
interpreted as the real and complex parts of a complex
|
|
kusano |
7d535a |
conjugate pair; A thus becomes block diagonal, with 1x1
|
|
kusano |
7d535a |
and 2x2 blocks.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
3. If UPPER='T', the upper triangle of A is set to random values
|
|
kusano |
7d535a |
out of distribution DIST.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
4. If SIM='T', A is multiplied on the left by a random matrix
|
|
kusano |
7d535a |
X, whose singular values are specified by DS, MODES, and
|
|
kusano |
7d535a |
CONDS, and on the right by X inverse.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
5. If KL < N-1, the lower bandwidth is reduced to KL using
|
|
kusano |
7d535a |
Householder transformations. If KU < N-1, the upper
|
|
kusano |
7d535a |
bandwidth is reduced to KU.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
6. If ANORM is not negative, the matrix is scaled to have
|
|
kusano |
7d535a |
maximum-element-norm ANORM.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
(Note: since the matrix cannot be reduced beyond Hessenberg form,
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
no packing options are available.)
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Arguments
|
|
kusano |
7d535a |
=========
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
N - INTEGER
|
|
kusano |
7d535a |
The number of columns (or rows) of A. Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DIST - CHARACTER*1
|
|
kusano |
7d535a |
On entry, DIST specifies the type of distribution to be used
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
to generate the random eigen-/singular values, and for the
|
|
kusano |
7d535a |
upper triangle (see UPPER).
|
|
kusano |
7d535a |
'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
|
|
kusano |
7d535a |
'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
|
|
kusano |
7d535a |
'N' => NORMAL( 0, 1 ) ( 'N' for normal )
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
ISEED - INTEGER array, dimension ( 4 )
|
|
kusano |
7d535a |
On entry ISEED specifies the seed of the random number
|
|
kusano |
7d535a |
generator. They should lie between 0 and 4095 inclusive,
|
|
kusano |
7d535a |
and ISEED(4) should be odd. The random number generator
|
|
kusano |
7d535a |
uses a linear congruential sequence limited to small
|
|
kusano |
7d535a |
integers, and so should produce machine independent
|
|
kusano |
7d535a |
random numbers. The values of ISEED are changed on
|
|
kusano |
7d535a |
exit, and can be used in the next call to SLATME
|
|
kusano |
7d535a |
to continue the same random number sequence.
|
|
kusano |
7d535a |
Changed on exit.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
D - REAL array, dimension ( N )
|
|
kusano |
7d535a |
This array is used to specify the eigenvalues of A. If
|
|
kusano |
7d535a |
MODE=0, then D is assumed to contain the eigenvalues (but
|
|
kusano |
7d535a |
see the description of EI), otherwise they will be
|
|
kusano |
7d535a |
computed according to MODE, COND, DMAX, and RSIGN and
|
|
kusano |
7d535a |
placed in D.
|
|
kusano |
7d535a |
Modified if MODE is nonzero.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
MODE - INTEGER
|
|
kusano |
7d535a |
On entry this describes how the eigenvalues are to
|
|
kusano |
7d535a |
be specified:
|
|
kusano |
7d535a |
MODE = 0 means use D (with EI) as input
|
|
kusano |
7d535a |
MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
|
|
kusano |
7d535a |
MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
|
|
kusano |
7d535a |
MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
|
|
kusano |
7d535a |
MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
|
|
kusano |
7d535a |
MODE = 5 sets D to random numbers in the range
|
|
kusano |
7d535a |
( 1/COND , 1 ) such that their logarithms
|
|
kusano |
7d535a |
are uniformly distributed. Each odd-even pair
|
|
kusano |
7d535a |
of elements will be either used as two real
|
|
kusano |
7d535a |
eigenvalues or as the real and imaginary part
|
|
kusano |
7d535a |
of a complex conjugate pair of eigenvalues;
|
|
kusano |
7d535a |
the choice of which is done is random, with
|
|
kusano |
7d535a |
50-50 probability, for each pair.
|
|
kusano |
7d535a |
MODE = 6 set D to random numbers from same distribution
|
|
kusano |
7d535a |
as the rest of the matrix.
|
|
kusano |
7d535a |
MODE < 0 has the same meaning as ABS(MODE), except that
|
|
kusano |
7d535a |
the order of the elements of D is reversed.
|
|
kusano |
7d535a |
Thus if MODE is between 1 and 4, D has entries ranging
|
|
kusano |
7d535a |
from 1 to 1/COND, if between -1 and -4, D has entries
|
|
kusano |
7d535a |
ranging from 1/COND to 1,
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
COND - REAL
|
|
kusano |
7d535a |
On entry, this is used as described under MODE above.
|
|
kusano |
7d535a |
If used, it must be >= 1. Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DMAX - REAL
|
|
kusano |
7d535a |
If MODE is neither -6, 0 nor 6, the contents of D, as
|
|
kusano |
7d535a |
computed according to MODE and COND, will be scaled by
|
|
kusano |
7d535a |
DMAX / max(abs(D(i))). Note that DMAX need not be
|
|
kusano |
7d535a |
positive: if DMAX is negative (or zero), D will be
|
|
kusano |
7d535a |
scaled by a negative number (or zero).
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
EI - CHARACTER*1 array, dimension ( N )
|
|
kusano |
7d535a |
If MODE is 0, and EI(1) is not ' ' (space character),
|
|
kusano |
7d535a |
this array specifies which elements of D (on input) are
|
|
kusano |
7d535a |
real eigenvalues and which are the real and imaginary parts
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
of a complex conjugate pair of eigenvalues. The elements
|
|
kusano |
7d535a |
of EI may then only have the values 'R' and 'I'. If
|
|
kusano |
7d535a |
EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
|
|
kusano |
7d535a |
CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
|
|
kusano |
7d535a |
conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th
|
|
kusano |
7d535a |
eigenvalue is D(j) (i.e., real). EI(1) may not be 'I',
|
|
kusano |
7d535a |
nor may two adjacent elements of EI both have the value 'I'.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
If MODE is not 0, then EI is ignored. If MODE is 0 and
|
|
kusano |
7d535a |
EI(1)=' ', then the eigenvalues will all be real.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
RSIGN - CHARACTER*1
|
|
kusano |
7d535a |
If MODE is not 0, 6, or -6, and RSIGN='T', then the
|
|
kusano |
7d535a |
elements of D, as computed according to MODE and COND, will
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
be multiplied by a random sign (+1 or -1). If RSIGN='F',
|
|
kusano |
7d535a |
they will not be. RSIGN may only have the values 'T' or
|
|
kusano |
7d535a |
'F'.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
UPPER - CHARACTER*1
|
|
kusano |
7d535a |
If UPPER='T', then the elements of A above the diagonal
|
|
kusano |
7d535a |
(and above the 2x2 diagonal blocks, if A has complex
|
|
kusano |
7d535a |
eigenvalues) will be set to random numbers out of DIST.
|
|
kusano |
7d535a |
If UPPER='F', they will not. UPPER may only have the
|
|
kusano |
7d535a |
values 'T' or 'F'.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SIM - CHARACTER*1
|
|
kusano |
7d535a |
If SIM='T', then A will be operated on by a "similarity
|
|
kusano |
7d535a |
transform", i.e., multiplied on the left by a matrix X and
|
|
kusano |
7d535a |
on the right by X inverse. X = U S V, where U and V are
|
|
kusano |
7d535a |
random unitary matrices and S is a (diagonal) matrix of
|
|
kusano |
7d535a |
singular values specified by DS, MODES, and CONDS. If
|
|
kusano |
7d535a |
SIM='F', then A will not be transformed.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
DS - REAL array, dimension ( N )
|
|
kusano |
7d535a |
This array is used to specify the singular values of X,
|
|
kusano |
7d535a |
in the same way that D specifies the eigenvalues of A.
|
|
kusano |
7d535a |
If MODE=0, the DS contains the singular values, which
|
|
kusano |
7d535a |
may not be zero.
|
|
kusano |
7d535a |
Modified if MODE is nonzero.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
MODES - INTEGER
|
|
kusano |
7d535a |
CONDS - REAL
|
|
kusano |
7d535a |
Same as MODE and COND, but for specifying the diagonal
|
|
kusano |
7d535a |
of S. MODES=-6 and +6 are not allowed (since they would
|
|
kusano |
7d535a |
result in randomly ill-conditioned eigenvalues.)
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
KL - INTEGER
|
|
kusano |
7d535a |
This specifies the lower bandwidth of the matrix. KL=1
|
|
kusano |
7d535a |
specifies upper Hessenberg form. If KL is at least N-1,
|
|
kusano |
7d535a |
then A will have full lower bandwidth. KL must be at
|
|
kusano |
7d535a |
least 1.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
KU - INTEGER
|
|
kusano |
7d535a |
This specifies the upper bandwidth of the matrix. KU=1
|
|
kusano |
7d535a |
specifies lower Hessenberg form. If KU is at least N-1,
|
|
kusano |
7d535a |
then A will have full upper bandwidth; if KU and KL
|
|
kusano |
7d535a |
are both at least N-1, then A will be dense. Only one of
|
|
kusano |
7d535a |
KU and KL may be less than N-1. KU must be at least 1.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
ANORM - REAL
|
|
kusano |
7d535a |
If ANORM is not negative, then A will be scaled by a non-
|
|
kusano |
7d535a |
negative real number to make the maximum-element-norm of A
|
|
kusano |
7d535a |
to be ANORM.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
A - REAL array, dimension ( LDA, N )
|
|
kusano |
7d535a |
On exit A is the desired test matrix.
|
|
kusano |
7d535a |
Modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
LDA - INTEGER
|
|
kusano |
7d535a |
LDA specifies the first dimension of A as declared in the
|
|
kusano |
7d535a |
calling program. LDA must be at least N.
|
|
kusano |
7d535a |
Not modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
WORK - REAL array, dimension ( 3*N )
|
|
kusano |
7d535a |
Workspace.
|
|
kusano |
7d535a |
Modified.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
INFO - INTEGER
|
|
kusano |
7d535a |
Error code. On exit, INFO will be set to one of the
|
|
kusano |
7d535a |
following values:
|
|
kusano |
7d535a |
0 => normal return
|
|
kusano |
7d535a |
-1 => N negative
|
|
kusano |
7d535a |
-2 => DIST illegal string
|
|
kusano |
7d535a |
-5 => MODE not in range -6 to 6
|
|
kusano |
7d535a |
-6 => COND less than 1.0, and MODE neither -6, 0 nor 6
|
|
kusano |
7d535a |
-8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
two adjacent elements of EI are 'I'.
|
|
kusano |
7d535a |
-9 => RSIGN is not 'T' or 'F'
|
|
kusano |
7d535a |
-10 => UPPER is not 'T' or 'F'
|
|
kusano |
7d535a |
-11 => SIM is not 'T' or 'F'
|
|
kusano |
7d535a |
-12 => MODES=0 and DS has a zero singular value.
|
|
kusano |
7d535a |
-13 => MODES is not in the range -5 to 5.
|
|
kusano |
7d535a |
-14 => MODES is nonzero and CONDS is less than 1.
|
|
kusano |
7d535a |
-15 => KL is less than 1.
|
|
kusano |
7d535a |
-16 => KU is less than 1, or KL and KU are both less than
|
|
kusano |
7d535a |
N-1.
|
|
kusano |
7d535a |
-19 => LDA is less than N.
|
|
kusano |
7d535a |
1 => Error return from SLATM1 (computing D)
|
|
kusano |
7d535a |
2 => Cannot scale to DMAX (max. eigenvalue is 0)
|
|
kusano |
7d535a |
3 => Error return from SLATM1 (computing DS)
|
|
kusano |
7d535a |
4 => Error return from SLARGE
|
|
kusano |
7d535a |
5 => Zero singular value from SLATM1.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
=====================================================================
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
1) Decode and Test the input parameters.
|
|
kusano |
7d535a |
Initialize flags & seed.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Parameter adjustments */
|
|
kusano |
7d535a |
--iseed;
|
|
kusano |
7d535a |
--d;
|
|
kusano |
7d535a |
--ei;
|
|
kusano |
7d535a |
--ds;
|
|
kusano |
7d535a |
a_dim1 = *lda;
|
|
kusano |
7d535a |
a_offset = a_dim1 + 1;
|
|
kusano |
7d535a |
a -= a_offset;
|
|
kusano |
7d535a |
--work;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Function Body */
|
|
kusano |
7d535a |
*info = 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Quick return if possible */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*n == 0) {
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Decode DIST */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (lsame_(dist, "U")) {
|
|
kusano |
7d535a |
idist = 1;
|
|
kusano |
7d535a |
} else if (lsame_(dist, "S")) {
|
|
kusano |
7d535a |
idist = 2;
|
|
kusano |
7d535a |
} else if (lsame_(dist, "N")) {
|
|
kusano |
7d535a |
idist = 3;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
idist = -1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Check EI */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
useei = TRUE_;
|
|
kusano |
7d535a |
badei = FALSE_;
|
|
kusano |
7d535a |
if (lsame_(ei + 1, " ") || *mode != 0) {
|
|
kusano |
7d535a |
useei = FALSE_;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
if (lsame_(ei + 1, "R")) {
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 2; j <= i__1; ++j) {
|
|
kusano |
7d535a |
if (lsame_(ei + j, "I")) {
|
|
kusano |
7d535a |
if (lsame_(ei + (j - 1), "I")) {
|
|
kusano |
7d535a |
badei = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
if (! lsame_(ei + j, "R")) {
|
|
kusano |
7d535a |
badei = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* L10: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
badei = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Decode RSIGN */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (lsame_(rsign, "T")) {
|
|
kusano |
7d535a |
irsign = 1;
|
|
kusano |
7d535a |
} else if (lsame_(rsign, "F")) {
|
|
kusano |
7d535a |
irsign = 0;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
irsign = -1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Decode UPPER */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (lsame_(upper, "T")) {
|
|
kusano |
7d535a |
iupper = 1;
|
|
kusano |
7d535a |
} else if (lsame_(upper, "F")) {
|
|
kusano |
7d535a |
iupper = 0;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
iupper = -1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Decode SIM */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (lsame_(sim, "T")) {
|
|
kusano |
7d535a |
isim = 1;
|
|
kusano |
7d535a |
} else if (lsame_(sim, "F")) {
|
|
kusano |
7d535a |
isim = 0;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
isim = -1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Check DS, if MODES=0 and ISIM=1 */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
bads = FALSE_;
|
|
kusano |
7d535a |
if (*modes == 0 && isim == 1) {
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 1; j <= i__1; ++j) {
|
|
kusano |
7d535a |
if (ds[j] == 0.f) {
|
|
kusano |
7d535a |
bads = TRUE_;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* L20: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Set INFO if an error */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*n < 0) {
|
|
kusano |
7d535a |
*info = -1;
|
|
kusano |
7d535a |
} else if (idist == -1) {
|
|
kusano |
7d535a |
*info = -2;
|
|
kusano |
7d535a |
} else if (abs(*mode) > 6) {
|
|
kusano |
7d535a |
*info = -5;
|
|
kusano |
7d535a |
} else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
|
|
kusano |
7d535a |
*info = -6;
|
|
kusano |
7d535a |
} else if (badei) {
|
|
kusano |
7d535a |
*info = -8;
|
|
kusano |
7d535a |
} else if (irsign == -1) {
|
|
kusano |
7d535a |
*info = -9;
|
|
kusano |
7d535a |
} else if (iupper == -1) {
|
|
kusano |
7d535a |
*info = -10;
|
|
kusano |
7d535a |
} else if (isim == -1) {
|
|
kusano |
7d535a |
*info = -11;
|
|
kusano |
7d535a |
} else if (bads) {
|
|
kusano |
7d535a |
*info = -12;
|
|
kusano |
7d535a |
} else if (isim == 1 && abs(*modes) > 5) {
|
|
kusano |
7d535a |
*info = -13;
|
|
kusano |
7d535a |
} else if (isim == 1 && *modes != 0 && *conds < 1.f) {
|
|
kusano |
7d535a |
*info = -14;
|
|
kusano |
7d535a |
} else if (*kl < 1) {
|
|
kusano |
7d535a |
*info = -15;
|
|
kusano |
7d535a |
} else if (*ku < 1 || *ku < *n - 1 && *kl < *n - 1) {
|
|
kusano |
7d535a |
*info = -16;
|
|
kusano |
7d535a |
} else if (*lda < max(1,*n)) {
|
|
kusano |
7d535a |
*info = -19;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*info != 0) {
|
|
kusano |
7d535a |
i__1 = -(*info);
|
|
kusano |
7d535a |
xerbla_("SLATME", &i__1);
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Initialize random number generator */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
for (i = 1; i <= 4; ++i) {
|
|
kusano |
7d535a |
iseed[i] = (i__1 = iseed[i], abs(i__1)) % 4096;
|
|
kusano |
7d535a |
/* L30: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (iseed[4] % 2 != 1) {
|
|
kusano |
7d535a |
++iseed[4];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* 2) Set up diagonal of A
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Compute D according to COND and MODE */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
slatm1_(mode, cond, &irsign, &idist, &iseed[1], &d[1], n, &iinfo);
|
|
kusano |
7d535a |
if (iinfo != 0) {
|
|
kusano |
7d535a |
*info = 1;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
if (*mode != 0 && abs(*mode) != 6) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Scale by DMAX */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
temp = dabs(d[1]);
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (i = 2; i <= i__1; ++i) {
|
|
kusano |
7d535a |
/* Computing MAX */
|
|
kusano |
7d535a |
r__2 = temp, r__3 = (r__1 = d[i], dabs(r__1));
|
|
kusano |
7d535a |
temp = dmax(r__2,r__3);
|
|
kusano |
7d535a |
/* L40: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (temp > 0.f) {
|
|
kusano |
7d535a |
alpha = *dmax__ / temp;
|
|
kusano |
7d535a |
} else if (*dmax__ != 0.f) {
|
|
kusano |
7d535a |
*info = 2;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
alpha = 0.f;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
sscal_(n, &alpha, &d[1], &c__1);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
slaset_("Full", n, n, &c_b23, &c_b23, &a[a_offset], lda);
|
|
kusano |
7d535a |
i__1 = *lda + 1;
|
|
kusano |
7d535a |
scopy_(n, &d[1], &c__1, &a[a_offset], &i__1);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Set up complex conjugate pairs */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*mode == 0) {
|
|
kusano |
7d535a |
if (useei) {
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 2; j <= i__1; ++j) {
|
|
kusano |
7d535a |
if (lsame_(ei + j, "I")) {
|
|
kusano |
7d535a |
a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
|
|
kusano |
7d535a |
a[j + (j - 1) * a_dim1] = -(doublereal)a[j + j * a_dim1];
|
|
kusano |
7d535a |
a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* L50: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} else if (abs(*mode) == 5) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 2; j <= i__1; j += 2) {
|
|
kusano |
7d535a |
if (slaran_(&iseed[1]) > .5f) {
|
|
kusano |
7d535a |
a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
|
|
kusano |
7d535a |
a[j + (j - 1) * a_dim1] = -(doublereal)a[j + j * a_dim1];
|
|
kusano |
7d535a |
a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* L60: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* 3) If UPPER='T', set upper triangle of A to random numbers.
|
|
kusano |
7d535a |
(but don't modify the corners of 2x2 blocks.) */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (iupper != 0) {
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (jc = 2; jc <= i__1; ++jc) {
|
|
kusano |
7d535a |
if (a[jc - 1 + jc * a_dim1] != 0.f) {
|
|
kusano |
7d535a |
jr = jc - 2;
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
jr = jc - 1;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
slarnv_(&idist, &iseed[1], &jr, &a[jc * a_dim1 + 1]);
|
|
kusano |
7d535a |
/* L70: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* 4) If SIM='T', apply similarity transformation.
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
-1
|
|
kusano |
7d535a |
Transform is X A X , where X = U S V, thus
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
it is U S V A V' (1/S) U' */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (isim != 0) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Compute S (singular values of the eigenvector matrix)
|
|
kusano |
7d535a |
according to CONDS and MODES */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
slatm1_(modes, conds, &c__0, &c__0, &iseed[1], &ds[1], n, &iinfo);
|
|
kusano |
7d535a |
if (iinfo != 0) {
|
|
kusano |
7d535a |
*info = 3;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Multiply by V and V' */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
slarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
|
|
kusano |
7d535a |
if (iinfo != 0) {
|
|
kusano |
7d535a |
*info = 4;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Multiply by S and (1/S) */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 1; j <= i__1; ++j) {
|
|
kusano |
7d535a |
sscal_(n, &ds[j], &a[j + a_dim1], lda);
|
|
kusano |
7d535a |
if (ds[j] != 0.f) {
|
|
kusano |
7d535a |
r__1 = 1.f / ds[j];
|
|
kusano |
7d535a |
sscal_(n, &r__1, &a[j * a_dim1 + 1], &c__1);
|
|
kusano |
7d535a |
} else {
|
|
kusano |
7d535a |
*info = 5;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/* L80: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Multiply by U and U' */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
slarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
|
|
kusano |
7d535a |
if (iinfo != 0) {
|
|
kusano |
7d535a |
*info = 4;
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* 5) Reduce the bandwidth. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*kl < *n - 1) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Reduce bandwidth -- kill column */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
i__1 = *n - 1;
|
|
kusano |
7d535a |
for (jcr = *kl + 1; jcr <= i__1; ++jcr) {
|
|
kusano |
7d535a |
ic = jcr - *kl;
|
|
kusano |
7d535a |
irows = *n + 1 - jcr;
|
|
kusano |
7d535a |
icols = *n + *kl - jcr;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
scopy_(&irows, &a[jcr + ic * a_dim1], &c__1, &work[1], &c__1);
|
|
kusano |
7d535a |
xnorms = work[1];
|
|
kusano |
7d535a |
slarfg_(&irows, &xnorms, &work[2], &c__1, &tau);
|
|
kusano |
7d535a |
work[1] = 1.f;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
sgemv_("T", &irows, &icols, &c_b39, &a[jcr + (ic + 1) * a_dim1],
|
|
kusano |
7d535a |
lda, &work[1], &c__1, &c_b23, &work[irows + 1], &c__1)
|
|
kusano |
7d535a |
;
|
|
kusano |
7d535a |
r__1 = -(doublereal)tau;
|
|
kusano |
7d535a |
sger_(&irows, &icols, &r__1, &work[1], &c__1, &work[irows + 1], &
|
|
kusano |
7d535a |
c__1, &a[jcr + (ic + 1) * a_dim1], lda);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
sgemv_("N", n, &irows, &c_b39, &a[jcr * a_dim1 + 1], lda, &work[1]
|
|
kusano |
7d535a |
, &c__1, &c_b23, &work[irows + 1], &c__1);
|
|
kusano |
7d535a |
r__1 = -(doublereal)tau;
|
|
kusano |
7d535a |
sger_(n, &irows, &r__1, &work[irows + 1], &c__1, &work[1], &c__1,
|
|
kusano |
7d535a |
&a[jcr * a_dim1 + 1], lda);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
a[jcr + ic * a_dim1] = xnorms;
|
|
kusano |
7d535a |
i__2 = irows - 1;
|
|
kusano |
7d535a |
slaset_("Full", &i__2, &c__1, &c_b23, &c_b23, &a[jcr + 1 + ic *
|
|
kusano |
7d535a |
a_dim1], lda);
|
|
kusano |
7d535a |
/* L90: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
} else if (*ku < *n - 1) {
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Reduce upper bandwidth -- kill a row at a time. */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
i__1 = *n - 1;
|
|
kusano |
7d535a |
for (jcr = *ku + 1; jcr <= i__1; ++jcr) {
|
|
kusano |
7d535a |
ir = jcr - *ku;
|
|
kusano |
7d535a |
irows = *n + *ku - jcr;
|
|
kusano |
7d535a |
icols = *n + 1 - jcr;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
scopy_(&icols, &a[ir + jcr * a_dim1], lda, &work[1], &c__1);
|
|
kusano |
7d535a |
xnorms = work[1];
|
|
kusano |
7d535a |
slarfg_(&icols, &xnorms, &work[2], &c__1, &tau);
|
|
kusano |
7d535a |
work[1] = 1.f;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
sgemv_("N", &irows, &icols, &c_b39, &a[ir + 1 + jcr * a_dim1],
|
|
kusano |
7d535a |
lda, &work[1], &c__1, &c_b23, &work[icols + 1], &c__1)
|
|
kusano |
7d535a |
;
|
|
kusano |
7d535a |
r__1 = -(doublereal)tau;
|
|
kusano |
7d535a |
sger_(&irows, &icols, &r__1, &work[icols + 1], &c__1, &work[1], &
|
|
kusano |
7d535a |
c__1, &a[ir + 1 + jcr * a_dim1], lda);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
sgemv_("C", &icols, n, &c_b39, &a[jcr + a_dim1], lda, &work[1], &
|
|
kusano |
7d535a |
c__1, &c_b23, &work[icols + 1], &c__1);
|
|
kusano |
7d535a |
r__1 = -(doublereal)tau;
|
|
kusano |
7d535a |
sger_(&icols, n, &r__1, &work[1], &c__1, &work[icols + 1], &c__1,
|
|
kusano |
7d535a |
&a[jcr + a_dim1], lda);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
a[ir + jcr * a_dim1] = xnorms;
|
|
kusano |
7d535a |
i__2 = icols - 1;
|
|
kusano |
7d535a |
slaset_("Full", &c__1, &i__2, &c_b23, &c_b23, &a[ir + (jcr + 1) *
|
|
kusano |
7d535a |
a_dim1], lda);
|
|
kusano |
7d535a |
/* L100: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Scale the matrix to have norm ANORM */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (*anorm >= 0.f) {
|
|
kusano |
7d535a |
temp = slange_("M", n, n, &a[a_offset], lda, tempa);
|
|
kusano |
7d535a |
if (temp > 0.f) {
|
|
kusano |
7d535a |
alpha = *anorm / temp;
|
|
kusano |
7d535a |
i__1 = *n;
|
|
kusano |
7d535a |
for (j = 1; j <= i__1; ++j) {
|
|
kusano |
7d535a |
sscal_(n, &alpha, &a[j * a_dim1 + 1], &c__1);
|
|
kusano |
7d535a |
/* L110: */
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
return 0;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* End of SLATME */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} /* slatme_ */
|
|
kusano |
7d535a |
|