clatms.c

SuperLU is a general purpose library for the direct solution of large, sparse, nonsymmetric systems

/*  -- translated by f2c (version 19940927).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Table of constant values */

static complex c_b1 = {0.f,0.f};
static integer c__1 = 1;
static integer c__5 = 5;
static logical c_true = TRUE_;
static logical c_false = FALSE_;

/* Subroutine */ int clatms_(integer *m, integer *n, char *dist, integer *
	iseed, char *sym, real *d, integer *mode, real *cond, real *dmax__, 
	integer *kl, integer *ku, char *pack, complex *a, integer *lda, 
	complex *work, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
    real r__1, r__2, r__3;
    doublereal d__1;
    complex q__1, q__2, q__3;
    logical L__1;

    /* Builtin functions */
    double cos(doublereal), sin(doublereal);
    void r_cnjg(complex *, complex *);

    /* Local variables */
    static integer ilda, icol;
    static real temp;
    static logical csym;
    static integer irow, isym;
    static complex c;
    static integer i, j, k;
    static complex s;
    static real alpha, angle;
    static integer ipack;
    static real realc;
    static integer ioffg;
    extern logical lsame_(char *, char *);
    static integer iinfo;
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static complex ctemp;
    static integer idist, mnmin, iskew;
    static complex extra, dummy;
    extern /* Subroutine */ int slatm1_(integer *, real *, integer *, integer 
	    *, integer *, real *, integer *, integer *);
    static integer ic, jc, nc;
    extern /* Subroutine */ int clagge_(integer *, integer *, integer *, 
	    integer *, real *, complex *, integer *, integer *, complex *, 
	    integer *), claghe_(integer *, integer *, real *, complex *, 
	    integer *, integer *, complex *, integer *);
    static integer il;
    static complex ct;
    static integer iendch, ir, jr, ipackg, mr;
    extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
    static integer minlda;
    static complex st;
    extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
	    *, complex *, complex *, integer *), clartg_(complex *, 
	    complex *, real *, complex *, complex *), xerbla_(char *, integer 
	    *), clagsy_(integer *, integer *, real *, complex *, 
	    integer *, integer *, complex *, integer *);
    extern doublereal slarnd_(integer *, integer *);
    extern /* Subroutine */ int clarot_(logical *, logical *, logical *, 
	    integer *, complex *, complex *, complex *, integer *, complex *, 
	    complex *);
    static logical iltemp, givens;
    static integer ioffst, irsign;
    static logical ilextr, topdwn;
    static integer ir1, ir2, isympk, jch, llb, jkl, jku, uub;


/*  -- LAPACK test routine (version 2.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       September 30, 1994   


    Purpose   
    =======   

       CLATMS generates random matrices with specified singular values   
       (or hermitian with specified eigenvalues)   
       for testing LAPACK programs.   

       CLATMS 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 hermitian or 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 hermitian or 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 hermitian)   
            zero out lower half (if hermitian)   
            store the upper half columnwise (if hermitian or upper   
                  triangular)   
            store the lower half columnwise (if hermitian or lower   
                  triangular)   
            store the lower triangle in banded format (if hermitian or   
                  lower triangular)   
            store the upper triangle in banded format (if hermitian 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.   

    Arguments   
    =========   

    M      - INTEGER   
             The number of rows of A. Not modified.   

    N      - INTEGER   
             The number of columns of A. N must equal M if the matrix   
             is symmetric or hermitian (i.e., if SYM is not 'N')   
             Not modified.   

    DIST   - 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.   

    ISEED  - 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 CLATMS   
             to continue the same random number sequence.   
             Changed on exit.   

    SYM    - CHARACTER*1   
             If SYM='H', the generated matrix is hermitian, with   
               eigenvalues specified by D, COND, MODE, and DMAX; they   
               may be positive, negative, or zero.   
             If SYM='P', the generated matrix is hermitian, 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.   
             If SYM='S', the generated matrix is (complex) symmetric,   
               with singular values specified by D, COND, MODE, and   
               DMAX; they will not be negative.   
             Not modified.   

    D      - REAL 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.   

    MODE   - 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='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.   

    COND   - REAL   
             On entry, this is used as described under MODE above.   
             If used, it must be >= 1. Not modified.   

    DMAX   - REAL   
             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.   

    KL     - 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 or hermitian.   
             Not modified.   

    KU     - 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 or hermitian.   
             Not modified.   

    PACK   - CHARACTER*1   
             This specifies packing of matrix as follows:   
             'N' => no packing   
             'U' => zero out all subdiagonal entries (if symmetric   
                    or hermitian)   
             'L' => zero out all superdiagonal entries (if symmetric   
                    or hermitian)   
             'C' => store the upper triangle columnwise (only if the   
                    matrix is symmetric, hermitian, or upper triangular) 
  
             'R' => store the lower triangle columnwise (only if the   
                    matrix is symmetric, hermitian, or lower triangular) 
  
             'B' => store the lower triangle in band storage scheme   
                    (only if the matrix is symmetric, hermitian, or   
                    lower triangular)   
             'Q' => store the upper triangle in band storage scheme   
                    (only if the matrix is symmetric, hermitian, 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, HB, or TB     - use 'B' or 'Q'   
             PP, SP, HB, or TP     - use 'C' or 'R'   

             If two calls to CLATMS differ only in the PACK parameter,   
             they will generate mathematically equivalent matrices.   
             Not modified.   

    A      - COMPLEX 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.   

    LDA    - 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.   

    WORK   - COMPLEX array, dimension ( 3*MAX( N, M ) )   
             Workspace.   
             Modified.   

    INFO   - 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 is not 'N' and KU is 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 SLATM1   
              2  => Cannot scale to DMAX (max. sing. value is 0)   
              3  => Error return from CLAGGE, CLAGHE or CLAGSY   

    ===================================================================== 
  


       1)      Decode and Test the input parameters.   
               Initialize flags & seed.   

       Parameter adjustments */
    --iseed;
    --d;
    a_dim1 = *lda;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --work;

    /* Function Body */
    *info = 0;

/*     Quick return if possible */

    if (*m == 0 || *n == 0) {
	return 0;
    }

/*     Decode DIST */

    if (lsame_(dist, "U")) {
	idist = 1;
    } else if (lsame_(dist, "S")) {
	idist = 2;
    } else if (lsame_(dist, "N")) {
	idist = 3;
    } else {
	idist = -1;
    }

/*     Decode SYM */

    if (lsame_(sym, "N")) {
	isym = 1;
	irsign = 0;
	csym = FALSE_;
    } else if (lsame_(sym, "P")) {
	isym = 2;
	irsign = 0;
	csym = FALSE_;
    } else if (lsame_(sym, "S")) {
	isym = 2;
	irsign = 0;
	csym = TRUE_;
    } else if (lsame_(sym, "H")) {
	isym = 2;
	irsign = 1;
	csym = FALSE_;
    } else {
	isym = -1;
    }

/*     Decode PACK */

    isympk = 0;
    if (lsame_(pack, "N")) {
	ipack = 0;
    } else if (lsame_(pack, "U")) {
	ipack = 1;
	isympk = 1;
    } else if (lsame_(pack, "L")) {
	ipack = 2;
	isympk = 1;