## C specifications:

**acmlversion**(

*int *major, int *minor, int *patch*);

— function: void

**acmlinfo**(

*void*);

## CFFT1D Routine Documentation

### SUBROUTINE: **CFFT1D** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT1D`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1D`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1D`

.MODE=-2 : (default) initializations and a forward transform are performed.MODE=2 : (default) initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis the length of the complex sequenceX

#### — Input/Output: COMPLEX **X**(*N*)

On input:

Xcontains the complex sequence of lengthNto be transformed.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX **COMM**(*5*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL CFFT1D(0,N,X,COMM,INFO) CALL CFFT1D(-1,N,X,COMM,INFO) CALL CFFT1D(-1,N,Y,COMM,INFO) DO 10 I = 1, N X(I) = X(I)*CONJG(Y(I)) 10 CONTINUE CALL CFFT1D(1,N,X,COMM,INFO)

## CFFT1DX Routine Documentation

### SUBROUTINE: **CFFT1DX** (*MODE,SCALE,INPL,N,X,INCX,Y,INCY,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT1DX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1DX`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1DX`

.MODE=-2 : (default) initializations and a forward transform are performed.MODE=2 : (default) initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequence

#### — Input: LOGICAL **INPL**

On input: if

INPLis .TRUE. thenXis overwritten by the output sequence; otherwise the output sequence is returned inY.

#### — Input: INTEGER **N**

On input:

Nis the number of elements to be transformed

#### — Input/Output: COMPLEX **X**(*1+*(*N-1*)**INCX*)

On input:

Xcontains the complex sequence of lengthNto be transformed, with the ith element stored in X(1+(i-1)*INCX).

On output: ifINPLis .TRUE. thenXcontains the transformed sequence in the same locations as on input; otherwiseXremains unchanged.

#### — Input: INTEGER **INCX**

On input:

INCXis the increment used to store successive elements of a sequence inX.

Constraint:INCX> 0.

#### — Output: COMPLEX **Y**(*1+*(*N-1*)**INCY*)

On output: if

INPLis .FALSE. thenYcontains the transformed sequence, with the ith element stored in Y(1+(i-1)*INCY); otherwiseYis not referenced.

#### — Input: INTEGER **INCY**

On input:

INCYis the increment used to store successive elements of a sequence inY. IfINPLis .TRUE. thenINCYis not referenced.

Constraint:INCY> 0.

#### — Input/Output: COMPLEX **COMM**(*5*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward FFTs are performed unscaled and in-place on contiguous C vectors X and Y following initialization. Manipulations on C resultant Fourier coefficients are stored in X which is then C transformed back. C SCALE = 1.0 INPL = .TRUE. CALL CFFT1DX(0,SCALE,INPL,N,X,1,DUM,1,COMM,INFO) CALL CFFT1DX(-1,SCALE,INPL,N,X,1,DUM,1,COMM,INFO) CALL CFFT1DX(-1,SCALE,INPL,N,Y,1,DUM,1,COMM,INFO) DO 10 I = 1, N X(I) = X(I)*CONJG(Y(I))/REAL(N) 10 CONTINUE CALL CFFT1DX(1,SCALE,INPL,N,X,1,DUM,1,COMM,INFO)

## CFFT1M Routine Documentation

### SUBROUTINE:**CFFT1M** (*MODE,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT1M`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : forward transforms are performed. Initializations are assumed to have been performed by a prior call to`CFFT1M`

.MODE=1 : backward (reverse) transforms are performed. Initializations are assumed to have been performed by a prior call to`CFFT1M`

.MODE=-2 : (default) initializations and forward transforms are performed.MODE=2 : (default) initializations and backward transforms are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the complex sequences inX

#### — Input/Output: COMPLEX **X**(*N*M*)

On input:

Xcontains theMcomplex sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.

On output:Xcontains the transformed sequences.

#### — Input/Output: COMPLEX **COMM**(*5*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL CFFT1M(0,1,N,X,COMM,INFO) CALL CFFT1M(-1,2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*CONJG(X(I,2)) X(I,2) = CMPLX(0.0D0,1.0D0)*X(I,2) 10 CONTINUE CALL CFFT1M(1,2,N,X(1,2),COMM,INFO)

## CFFT1MX Routine Documentation

### SUBROUTINE:**CFFT1MX** (*MODE,SCALE,INPL,NSEQ,N,X,INCX1,INCX2,Y,INCY1,INCY2,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT1MX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1MX`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT1MX`

.MODE=-2 : (default) initializations and a forward transform are performed.MODE=2 : (default) initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **INPL**

On input: if

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **NSEQ**

On input:

NSEQis the number of sequences to be transformed

#### — Input: INTEGER **N**

On input:

Nis the number of elements in each sequence to be transformed

#### — Input/Output: COMPLEX **X**(*1+*(*N-1*)**INCX1+*(*NSEQ-1*)**INCX2*)

On input:

Xcontains theNSEQcomplex sequences of lengthNto be transformed; the ith element of sequence j is stored in X(1+(i-1)*INCX1+(j-1)*INCX2).

On output: ifINPLis .TRUE. thenXcontains the transformed sequences in the same locations as on input; otherwiseXremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the increment used to store successive elements of a given sequence inX(INCX1=1 for contiguous data).

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the increment used to store corresponding elements of successive sequences inX(INCX2=Nfor contiguous data).

Constraint:INCX2> 0.

#### — Output: COMPLEX **Y**(*1+*(*N-1*)**INCY1+*(*NSEQ-1*)**INCY2*)

On output: if

INPLis .FALSE. thenYcontains the transformed sequences with the ith element of sequence j stored in Y(1+(i-1)*INCY1+(j-1)*INCY2); otherwiseYis not referenced.

#### — Input: INTEGER **INCY1**

On input:

INCY1is the increment used to store successive elements of a given sequence inY. IfINPLis .TRUE. thenINCY1is not referenced.

Constraint:INCY1> 0.

#### — Input: INTEGER **INCY2**

On input:

INCY2is the increment used to store corresponding elements of successive sequences inY(INCY2=Nfor contiguous data). IfINPLis .TRUE. thenINCY2is not referenced.

Constraint:INCY2> 0.

#### — Input/Output: COMPLEX **COMM**(*5*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward FFTs are performed unscaled and in-place on two C contiguous vectors stored in the first two columns of X. C Manipulations are stored in 2nd and 3rd columns of X which are C then transformed back. C COMPLEX X(N,3) SCALE = 1.0 INPL = .TRUE. CALL CFFT1MX(0,SCALE,INPL,2,N,X,1,N,DUM,1,N,COMM,INFO) CALL CFFT1MX(-1,SCALE,INPL,2,N,X,1,N,DUM,1,N,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*CONJG(X(I,2))/REAL(N) X(I,2) = CMPLX(0.0D0,1.0D0)*X(I,2)/REAL(N) 10 CONTINUE CALL CFFT1MX(1,SCALE,INPL,2,N,X(1,2),1,N,DUM,1,N,COMM,INFO)

## CFFT2D Routine Documentation

### SUBROUTINE:**CFFT2D** (*MODE,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the direction of transform to be performed by`CFFT2D`

.

On input:

MODE=-1 : a forward 2D transform is performed.MODE=1 : a backward (reverse) 2D transform is performed.

#### — Input: INTEGER **M**

On input:

Mis the number of rows in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the first dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the second dimension ofX.

#### — Input/Output: COMPLEX **X**(*M*N*)

On input:

Xcontains theMbyNcomplex 2D array to be transformed. Element ij is stored in location i+(j-1)*MofX.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX **COMM**(*M*N+5**(*M+N*))

COMMis a communication array used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL CFFT2D(-1,M,N,X,COMM,INFO) DO 20 J = 1, N DO 10 I = 1, MIN(J-1,M) X(I,J) = 0.5D0*(X(I,J) + X(J,I)) X(J,I) = CONJG(X(I,J)) 10 CONTINUE 20 CONTINUE CALL CFFT2D(1,M,N,X,COMM,INFO)

## CFFT2DX Routine Documentation

### SUBROUTINE:**CFFT2DX** (*MODE,SCALE,LTRANS,INPL,M,N,X,INCX1,INCX2,Y,INCY1,INCY2,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT2DX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 2D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT2DX`

.MODE=1 : a backward (reverse) 2D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT2DX`

.MODE=-2 : (default) initializations and a forward 2D transform are performed.MODE=2 : (default) initializations and a backward 2D transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the values ofNandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofNandM.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **LTRANS**

On input: if

LTRANSis .TRUE. then a normal final transposition is performed internally to return transformed data consistent with the values for argumentsINPL,INCX1,INCX2,INCY1andINCY2. IfLTRANSis .FALSE. then the final transposition is not performed explicitly; the storage format on output is determined by whether the output data is stored contiguously or not – please see the output specifications forXandYfor details.

#### — Input: LOGICAL **INPL**

On input: if

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **M**

On input:

Mis the first dimension of the 2D transform.

#### — Input: INTEGER **N**

On input:

Nis the second dimension of the 2D transform.

#### — Input/Output: COMPLEX **X**(*1+*(*M-1*)**INCX1+*(*N-1*)**INCX2*)

On input:

Xcontains theMbyNcomplex 2D data array to be transformed; the (ij)th element is stored in X(1+(i-1)*INCX1+(j-1)*INCX2).

On output: ifINPLis .TRUE. thenXcontains the transformed data, either in the same locations as on input whenLTRANS=.TRUE.; in locations X((i-1)*N+j) whenLTRANS=.FALSE.,INCX1=1 andINCX2=M; and otherwise in the same locations as on input. IfINPLis .FALSE.Xremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the increment used to store, inX, successive elements in the first dimension (INCX1=1 for contiguous data).

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the increment used to store, inX, successive elements in the second dimension (INCX2=M for contiguous data).

Constraint:INCX2> 0;

INCX2> (M-1)*INCX1 ifN> 1.

#### — Output: COMPLEX **Y**(*1+*(*M-1*)**INCY1+*(*N-1*)**INCY2*)

On output: if

INPLis .FALSE. thenYcontains the transformed data. IfLTRANS=.TRUE. then the (ij)th data element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2); ifLTRANS=.FALSE.,INCY1=1 andINCY2=N then the (ij)th data element is stored in Y((i-1)*N+j); and otherwise the (ij)th element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2). IfINPLis .TRUE. thenYis not referenced.

#### — Input: INTEGER **INCY1**

On input:

INCY1is the increment used to store successive elements in the first dimension inY(INCY1=1 for contiguous data). IfINPLis .TRUE. thenINCY1is not referenced.

Constraint:INCY1> 0.

#### — Input: INTEGER **INCY2**

On input:

INCY2is the increment used to store successive elements in the second dimension inY(for contiguous data,INCY2=MwhenLTRANSis .TRUE. orINCY2=NwhenLTRANSis .FALSE.). IfINPLis .TRUE. thenINCY2is not referenced.

Constraints:INCY2> 0;

INCY2> (M-1)*INCY1 ifN> 1 andLTRANSis .TRUE.;

INCY2= N ifM> 1 andLTRANSis .FALSE..

#### — Input/Output: COMPLEX **COMM**(*M*N+5*M+5*N+200*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same dimensionsMandN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 2D FFT is performed unscaled, without final transpose C and out-of-place on data stored in array X and output to Y. C Manipulations are stored in vector Y which is then transformed C back, with scaling, into the first M rows of X. C COMPLEX X(M,N), Y(N,M) SCALE = 1.0 INPL = .FALSE. LTRANS = .FALSE. CALL CFFT2DX(0,SCALE,LTRANS,INPL,M,N,X,1,M,Y,1,N,COMM,INFO) CALL CFFT2DX(-1,SCALE,LTRANS,INPL,M,N,X,1,M,Y,1,N,COMM,INFO) DO 20 I = M DO 10 J = 1, N Y(J,I) = 0.5*Y(J,I)*EXP(-0.001*REAL(I+J-2)) IY = IY + 1 10 CONTINUE 20 CONTINUE SCALE = 1.0/REAL(M*N) CALL CFFT2DX(1,SCALE,LTRANS,INPL,N,M,Y,1,N,X,1,M,COMM,INFO)

## CFFT3D Routine Documentation

### SUBROUTINE:**CFFT3D** (*MODE,L,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the direction of transform to be performed by`CFFT3D`

.

On input:

MODE=-1 : forward 3D transform is performed.MODE=1 : backward (reverse) 3D transform is performed.

#### — Input: INTEGER **L**

On input: the length of the first dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenLis the first dimension ofX.

#### — Input: INTEGER **M**

On input: the length of the second dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenMis the second dimension ofX.

#### — Input: INTEGER **N**

On input: the length of the third dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenNis the third dimension ofX.

#### — Input/Output: COMPLEX **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNcomplex 3D array to be transformed. Element ijk is stored in location i+(j-1)*L+(k-1)*L*MofX.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX **COMM**(*5**(*L+M+N*)*+4*)

COMMis a communication array used as temporary store. Note that the amount of store explicitly required here is less than in some versions prior to this release (version 4.1 and older). Some further workspace will be allocated internally; the amount of allocated memory requested will be MAX(L*N+N, L*M + L + M, N).

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL CFFT3D(-1,L,M,N,X,COMM,INFO) DO 30 K = 1, N DO 20 J = 1, M DO 10 I = 1, L X(I,J) = X(I,J)*EXP(-0.001D0*REAL(I+J+K)) 10 CONTINUE 20 CONTINUE 30 CONTINUE CALL CFFT3D(1,L,M,N,X,COMM,INFO)

## CFFT3DX Routine Documentation

### SUBROUTINE:**CFFT3DX** (*MODE,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT3DX`

.

On input:

MODE=0 : only initializations (specific to the values ofL,MandN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 3D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT3DX`

.MODE=1 : a backward (reverse) 3D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT3DX`

.MODE=100 : similar toMODE=0; only initializations (specific to the values ofL,MandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofL,MandN.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **LTRANS**

On input: if

LTRANSis .TRUE. then a normal final transposition is performed internally to return transformed data using the same storage format as the input data. IfLTRANSis .FALSE. then the final transposition is not performed and transformed data is stored, inXorY, in transposed form.

#### — Input: LOGICAL **INPL**

On input: if

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **L**

On input:

Lis the first dimension of the 3D transform.

#### — Input: INTEGER **M**

On input:

Mis the second dimension of the 3D transform.

#### — Input: INTEGER **N**

On input:

Nis the third dimension of the 3D transform.

#### — Input/Output: COMPLEX **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNcomplex 3D data array to be transformed; the (ijk)th element is stored in X(i+(j-1)*L+(k-1)*L*M).

On output: ifINPLis .TRUE. thenXcontains the transformed data, either in the same locations as on input whenLTRANS=.TRUE.; or in locations X(k+(j-1)*N+(i-1)*N*M) whenLTRANS=.FALSE. IfINPLis .FALSE.Xremains unchanged.

#### — Output: COMPLEX **Y**(*L*M*N*)

On output: if

INPLis .FALSE. thenYcontains the three-dimensional transformed data. IfLTRANS=.TRUE. then the (ijk)th data element is stored in Y(i+(j-1)*L+(k-1)*L*M); otherwise, the (ijk)th data element is stored in Y(k+(j-1)*N+(k-1)*N*M). IfINPLis .TRUE. thenYis not referenced.

#### — Input/Output: COMPLEX **COMM**(***)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence dimensions. The remainder is used as temporary store. The amount of store required depends on the values of the argumentsL,M,NandLTRANS. IfLTRANS=.TRUE. then for a genuine 3D transform (all ofL,M,Ngreater than 1) the dimension ofCOMMneed only be 5*(L+M+N) + 150; in this case some further workspace will be allocated internally; the amount of allocated memory requested will be MAX(L*N+N, L*M + L + M, N). IfLTRANS=.FALSE. then for a genuine 3D transform the workspace requirement is considerably more to allow for the storage of an intermediate 3D transposed array; in this case the dimension ofCOMMmust be at least L*M*N+5*(L+M+N)+150. It is recommended that the appropriate 1D or 2D FFT routine be called when at least one ofL,MorNis 1.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 3D FFT is performed unscaled, without final transpose C and out-of-place on data stored in array X and output to Y. C Manipulations are stored in vector Y which is then transformed C back, with scaling, into the first M rows of X. C SCALE = 1.0 INPL = .FALSE. LTRANS = .FALSE. CALL CFFT3DX(0,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO) CALL CFFT3DX(-1,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO) IY = 1 DO 20 I = 1, L DO 40 J = 1, M DO 10 K = 1, N Y(IY) = Y(IY)*EXP(-0.001*REAL(I+J+K-3)) IY = IY + 1 10 CONTINUE 20 CONTINUE SCALE = 1.0/REAL(L*M*N) CALL CFFT3DX(1,SCALE,LTRANS,INPL,N,M,L,Y,X,COMM,INFO)

## CFFT3DY Routine Documentation

### SUBROUTINE:**CFFT3DY** (*MODE,SCALE,INPL,L,M,N,X,INCX1,INCX2,INCX3,Y,INCY1,INCY2,INCY3,COMM,LCOMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CFFT3DY`

.

On input:

MODE=0 : only initializations (specific to the values ofL,MandN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 3D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT3DY`

.MODE=1 : a backward (reverse) 3D transform is performed. Initializations are assumed to have been performed by a prior call to`CFFT3DY`

.MODE=-2 : (default) initializations and a forward 3D transform are performed.MODE=2 : (default) initializations and a backward 3D transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the values ofL,MandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofL,MandN.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **INPL**

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **L**

On input:

Lis the first dimension of the 3D transform.

#### — Input: INTEGER **M**

On input:

Mis the second dimension of the 3D transform.

#### — Input: INTEGER **N**

On input:

Nis the third dimension of the 3D transform.

#### — Input/Output: COMPLEX **X**(***)

On input:

Xcontains theLbyMbyNcomplex 3D data array to be transformed; the (ijk)th element is stored in X(1+(i-1)*INCX1+(j-1)*INCX2+(k-1)*INCX3).

On output: ifINPLis .TRUE. thenXcontains the transformed data in the same locations as on input. IfINPLis .FALSE.Xremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the step in index ofXbetween successive data elements in the first dimension of the 3D data. UsuallyINCX1=1 so that succesive elements in the first dimension are stored contiguously.

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the step in index ofXbetween successive data elements in the second dimension of the 3D data. For completely contiguous data (no gaps inX)INCX2should be set toL.

Constraint:INCX2> 0;

INCX2> (L-1)*INCX1 if max(M,N) > 1.

#### — Input: INTEGER **INCX3**

On input:

INCX3is the step in index ofXbetween successive data elements in the third dimension of the 3D data. For completely contiguous data (no gaps inX)INCX3should be set toL*M.

Constraint:INCX3> 0;

INCX3> (L-1)*INCX1+(M-1)*INCX2if N > 1.

#### — Output: COMPLEX **Y**(***)

On output: if

INPLis .FALSE. thenYcontains the three-dimensional transformed data. IfLTRANS=.TRUE. then the the (ijk)th element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2+(k-1)*INCY3).

IfINPLis .TRUE. thenYis not referenced.

#### — Input: INTEGER **INCY1**

On input: if

INPLis .FALSE. thenINCY1is the step in index ofYbetween successive data elements in the first dimension of the 3D transformed data. UsuallyINCY1=1 so that succesive elements in the first dimension are stored contiguously.

IfINPLis .TRUE. thenINCY1is not referenced. Constraint: IfINPLis .FALSE. thenINCY1> 0.

#### — Input: INTEGER **INCY2**

On input: if

INPLis .FALSE. thenINCY2is the step in index ofYbetween successive data elements in the second dimension of the 3D transformed data. For completely contiguous data (no gaps inY)INCY2should be set toL.

Constraint:INCY2> 0 ifINPLis .FALSE.;

INCY2> (L-1)*INCY1, ifINPLis .FALSE. and max(M,N) > 1.

#### — Input: INTEGER **INCY3**

On input: if

INPLis .FALSE. thenINCY3is the step in index ofYbetween successive data elements in the third dimension of the 3D transformed data. For completely contiguous data (no gaps inY)INCY3should be set toL*M.

Constraint:INCY3> 0 ifINPLis .FALSE.;

INCY3> (L-1)*INCY1+(M-1)*INCY2, ifINPLis .FALSE. and N > 1.

#### — Input/Output: COMPLEX **COMM**(*LCOMM*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence dimensions. The remainder is used as temporary store; if this is not sufficient for the requirements of the routine then temporary storage space will be dynamically allocated internally.

#### — Input: INTEGER **LCOMM**

On input:

LCOMMis the length of the communication arrayCOMM. The amount of internal dynamic allocation of temporary storage is dependent on the values of the increment arguments for arraysXandY. The amount is minimized when the increments for the output array are 1,LandL*Mrespectively, since this represents a contiguous array in which calculations can be performedd in-place.

Constraint: If min(L,M,N) > 1,LCOMM>= 5*(L+M+N) + 150 (otherwise see see CFFT2DX).

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 3D FFT is performed unscaled and in-place, on the leading C 10x10x10 submatrix of a larger 100x100x100 array of data. C The result is transformed back with scaling. C SCALE = 1.0 INPL = .TRUE. L = 10 M = 10 N = 10 LCOMM = 2000000 CALL CFFT3DY(0,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO) CALL CFFT3DY(-1,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO) IY = 1 DO 20 I = 1, L DO 40 J = 1, M DO 10 K = 1, N X(I,J,K) = X(I,J,K)*EXP(-0.001*REAL(I+J+K-3)) 10 CONTINUE 20 CONTINUE SCALE = 1.0/REAL(L*M*N) CALL CFFT3DY(1,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO)

## CSFFT Routine Documentation

### SUBROUTINE:**CSFFT** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CSFFT`

.

On input:

MODE=0 : only initializations (specific to the values ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real transform is performed. Initializations are assumed to have been performed by a prior call to`CSFFT`

.MODE=2 : (default) initializations and a real transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **N**

On input:

Nis the length of the sequence inX

#### — Input/Output: REAL **X**(*N*)

On input:

Xcontains the Hermitian sequence of lengthNto be transformed.

On output:Xcontains the transformed real sequence.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT(0,N,X,COMM,INFO) CALL SCFFT(1,N,X,COMM,INFO) DO 10 I = N/2+2, N X(I) = -X(I) 10 CONTINUE CALL CSFFT(2,N,X,COMM,INFO)

## CSFFT1D Routine Documentation

### SUBROUTINE:**CSFFT1D** (*MODE,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CSFFT1D`

.

On input:

MODE=0 : only initializations (specific to the values ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`CSFFT1D`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inY

#### — Input: COMPLEX **X**(*N/2+1*)

On input:

Xcontains the complex sequence to be transformed.

#### — Output: REAL **Y**(*N*)

On output:

Ycontains the transformed real sequence of lengthN.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT1D(0,N,X,Y,COMM,INFO) CALL SCFFT1D(1,N,X,Y,COMM,INFO) DO 10 I = 1, N/2+1 Y(I) = -Y(I)*EXP(-REAL(I-1)/REAL(N)) 10 CONTINUE CALL CSFFT1D(2,N,Y,X,COMM,INFO)

## CSFFT1M Routine Documentation

### SUBROUTINE:**CSFFT1M** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CSFFT1M`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`CSFFT1M`

.MODE=2 : (default) initializations and a complex-to-real transform a re performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inY

#### — Input: COMPLEX **X**((*N/2+1*)**M*)

On input:

Xcontains the Hermitian sequences to be transformed in complex-Hermitian storage.

#### — Output: REAL **Y**(*N*M*)

On output:

Ycontains theMreal sequences of lengthNtransformed. Elementiof sequencejis stored in locationi + (j – 1) * NofY.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT1M(0,M,N,X,Y,COMM,INFO) CALL SCFFT1M(1,M,N,X,Y,COMM,INFO) DO J = 1, M DO I = 1, N/2+1 Y(I,J) = -Y(I,J)*EXP(-REAL(I-1)/REAL(N)) END DO EMD DO CALL CSFFT1M(2,M,N,Y,X,COMM,INFO)

## CSFFT2D Routine Documentation

### SUBROUTINE:**CSFFT2D** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CSFFT2D`

.

On input:

MODE=0 : only default initializations (specific toMandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`CSFFT2D`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of rows in the 2D array of the real data obtained from the transform. IfYis declared as a 2D array thenMis the first dimension ofY.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of the real data obtained from the transform. IfYis declared as a 2D array thenNis the second dimension ofY.

#### — Input: COMPLEX **X**((*M/2+1*)**N*)

On input:

Xcontains the Hermitian sequences in complex-Hermitian storage to be transformed.

#### — Output: REAL **Y**(*M*N*)

On output:

Ycontains theMbyNreal 2D array obtained from the transform. Elementijis stored in locationi + (j – 1) * MofY.

#### — Input/Output: REAL **COMM**(*4*M+10*N+300*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthM*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT2D(0,M,N,X,Y,COMM,INFO) CALL SCFFT2D(1,M,N,X,Y,COMM,INFO) DO J = 1, N DO I = 1, M/2+1 Y(I,J) = -Y(I,J)/SQRT(REAL(M*N)) END DO EMD DO CALL CSFFT2D(2,M,N,Y,X,COMM,INFO)

## CSFFT3D Routine Documentation

### SUBROUTINE:**CSFFT3D** (*MODE,L,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`CSFFT3D`

.

On input:

MODE=0 : only default initializations (specific toL,MandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`CSFFT3D`

.MODE=2 : (default) initializations and a complex-to-real transform a re performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **L**

On input:

Lis the length of the first dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the first dimension ofY.

#### — Input: INTEGER **M**

On input:

Mis the length of the second dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the second dimension ofY.

#### — Input: INTEGER **N**

On input:

Nis the length of the third dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the third dimension ofY.

#### — Input: COMPLEX **X**((*L/2+1*)**M*N*)

On input:

Xcontains the Hermitian sequences in complex-Hermitian storage to be transformed.

#### — Output: REAL **Y**(*L*M*N*)

On output:

Ycontains theLbyMbyNreal 3D array obtained from the transform. Elementijkis stored in locationi + (j – 1) * L + (k – 1) * L * MofY.

#### — Input/Output: REAL **COMM**(*4*L+10*M+10*N+500*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthL*M*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT3D(0,L,M,N,X,Y,COMM,INFO) CALL SCFFT3D(1,L,M,N,X,Y,COMM,INFO) DO K = 1, N DO J = 1, M DO I = 1, L/2+1 Y(I,J,K) = -Y(I,J,K)/SQRT(REAL(L*M*N)) END DO END DO EMD DO CALL CSFFT3D(2,L,M,N,Y,X,COMM,INFO)

## CSFFTM Routine Documentation

### SUBROUTINE:**CSFFTM** (*M,N,X,COMM,INFO*)

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the sequences inX

#### — Input/Output: REAL **X**(*N*M*)

On input:

Xcontains theMHermitian sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.On output:

Xcontains the transformed real sequences.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFTM(1,N,X,COMM,INFO) CALL SCFFTM(2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*X(N-I+1,2) 10 CONTINUE CALL CSFFTM(1,N,X(1,3),COMM,INFO)

## DRANDBETA / SRANDBETA

Generates a vector of random variates from a beta distribution with probability density function, f(X), where: f(X) = [Gamma(A + B) * X^(A – 1) * (1 – X)^(B – 1)] / [Gamma(A) * Gamma(B)], if 0 <= X <= 1; A,B > 0.0, otherwise f(X) = 0.

*(Note that SRANDBETA is the single precision version of DRANDBETA. The argument lists of both routines are identical except that any double precision arguments of DRANDBETA are replaced in SRANDBETA by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDBETA** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: first parameter for the distribution.

Constraint:A>0.

#### — Input: DOUBLE PRECISION **B**

On input: second parameter for the distribution.

Constraint:B>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDBETA`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Beta distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Beta distribution CALL DRANDBETA(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDBINOMIAL / SRANDBINOMIAL

Generates a vector of random variates from a Binomial distribution with probability, f(X), defined by: f(X) = [M! * P^X * (1 – P)^(M – X)] / [X! * (M – 1)!], I=0,1,…,M

*(Note that SRANDBINOMIAL is the single precision version of DRANDBINOMIAL. The argument lists of both routines are identical except that any double precision arguments of DRANDBINOMIAL are replaced in SRANDBINOMIAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDBINOMIAL** (*N,M,P,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of trials.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **P**

On input: probability of success.

Constraint: 0<=P<1.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDBINOMIAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Binomial distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER M DOUBLE PRECISION P INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Binomial distribution CALL DRANDBINOMIAL(N,M,P,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDBINOMIALREFERENCE / SRANDBINOMIALREFERENCE

Initializes a reference vector for use with `DRANDGENERALDISCRETE`

. Reference vector is for a Binomial distribution with probability, f(X), defined by: f(X) = [M! * P^X * (1 – P)^(M – X)] / [X! * (M – 1)!], I=0,1,…,M

*(Note that SRANDBINOMIALREFERENCE is the single precision version of DRANDBINOMIALREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDBINOMIALREFERENCE are replaced in SRANDBINOMIALREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDBINOMIALREFERENCE** (*M,P,REF,LREF,INFO*)

#### — Input: INTEGER **M**

On input: number of trials.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **P**

On input: probability of success.

Constraint: 0<=P<1.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Binomial distribution using DRANDGENERALDISCRETE.

#### — Input/Output: INTEGER **LREF**

On input: either the length of the reference vector

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the Binomial distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER M DOUBLE PRECISION P INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDBINOMIALREFERENCE(M,P,REF,LREF,INFO) C Generate N variates from the Binomial distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDBLUMBLUMSHUB / SRANDBLUMBLUMSHUB

Allows direct access to the bit stream generated by the Blum-Blum-Shub generator.

*(Note that SRANDBLUMBLUMSHUB is the single precision version of DRANDBLUMBLUMSHUB. The argument lists of both routines are identical except that any double precision arguments of DRANDBLUMBLUMSHUB are replaced in SRANDBLUMBLUMSHUB by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDBLUMBLUMSHUB** (*N,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required. The total number of bits generated is 24

N.

Constraint:N>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDBLUMBLUMSHUB`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector holding the bit stream. The least significant 24 bits of each of the

X(i) contain the bit stream as generated by the Blum-Blum-Shub generator. The least significant bit ofX(1) is the first bit generated, the second least significant bit ofX(1) is the second bit generated etc.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

## DRANDCAUCHY / SRANDCAUCHY

Generates a vector of random variates from a Cauchy distribution with probability density function, f(X), where: f(X) = 1 / [Pi * B * (1 + (X – A / B)^2)].

*(Note that SRANDCAUCHY is the single precision version of DRANDCAUCHY. The argument lists of both routines are identical except that any double precision arguments of DRANDCAUCHY are replaced in SRANDCAUCHY by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDCAUCHY** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: median of the distribution.

#### — Input: DOUBLE PRECISION **B**

On input: semi-quartile range of the distribution.

Constraint:B>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDCAUCHY`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Cauchy distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Cauchy distribution CALL DRANDCAUCHY(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDCHISQUARED / SRANDCHISQUARED

Generates a vector of random variates from a chi-squared distribution with probability density function, f(X), where: f(X) = [X^(v / 2 – 1) * exp(-X / 2)] / [2^(v / 2) * (v / 2 – 1)!], if X > 0, otherwise f(X) = 0. Here v is the degrees of freedom, *DF*.

*(Note that SRANDCHISQUARED is the single precision version of DRANDCHISQUARED. The argument lists of both routines are identical except that any double precision arguments of DRANDCHISQUARED are replaced in SRANDCHISQUARED by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDCHISQUARED** (*N,DF,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **DF**

On input: degrees of freedom of the distribution.

Constraint:DF>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDCHISQUARED`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Chi-squared distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER DF DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) DF C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Chi-squared distribution CALL DRANDCHISQUARED(N,DF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDDISCRETEUNIFORM / SRANDDISCRETEUNIFORM

Generates a vector of random variates from a Uniform distribution with probability f(X) defined by: f(X) = 1 / [B – A], X=A,A+1,…,B

*(Note that SRANDDISCRETEUNIFORM is the single precision version of DRANDDISCRETEUNIFORM. The argument lists of both routines are identical except that any double precision arguments of DRANDDISCRETEUNIFORM are replaced in SRANDDISCRETEUNIFORM by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDDISCRETEUNIFORM** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **A**

On input: minimum for the distribution.

#### — Input: INTEGER **B**

On input: maximum for the distribution.

Constraint:B>=A.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDDISCRETEUNIFORM`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Uniform distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER A,B INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Uniform distribution CALL DRANDDISCRETEUNIFORM(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDEXPONENTIAL / SRANDEXPONENTIAL

Generates a vector of random variates from an exponential distribution with probability density function, f(X), where: f(X) = exp(-X / A) / A, if X>0, otherwise f(X)=0.

*(Note that SRANDEXPONENTIAL is the single precision version of DRANDEXPONENTIAL. The argument lists of both routines are identical except that any double precision arguments of DRANDEXPONENTIAL are replaced in SRANDEXPONENTIAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDEXPONENTIAL** (*N,A,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: exponential parameter.

Constraint:A>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDEXPONENTIAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Exponential distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Exponential distribution CALL DRANDEXPONENTIAL(N,A,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDF / SRANDF

Generates a vector of random variates from an F distribution, also called the Fisher’s variance ratio distribution, with probability density function, f(X), where: f(X) = [((m + v – 2) / 2)! * X^(m / 2 – 1) * m^(m / 2)] / [(m / 2 – 1)! * (v / 2 – 1)! * (1 + m * X / v)^((m + v) / 2) * v^(m / 2)], if X > 0, otherwise f(X) = 0. Here m is the first degrees of freedom, (*DF1*) and v is the second degrees of freedom, (*DF2*).

*(Note that SRANDF is the single precision version of DRANDF. The argument lists of both routines are identical except that any double precision arguments of DRANDF are replaced in SRANDF by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDF** (*N,DF1,DF2,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **DF1**

On input: first degrees of freedom.

Constraint:DF1>=0.

#### — Input: INTEGER **DF2**

On input: second degrees of freedom.

Constraint:DF2>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDF`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the F distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER DF1,DF2 DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) DF1,DF2 C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the F distribution CALL DRANDF(N,DF1,DF2,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDGAMMA / SRANDGAMMA

Generates a vector of random variates from a Gamma distribution with probability density function, f(X), where: f(X) = [X^(A – 1) * exp(-X / B)] / [B^A * Gamma(A)], if X >= 0 and A,B > 0.0, otherwise f(X) = 0.

*(Note that SRANDGAMMA is the single precision version of DRANDGAMMA. The argument lists of both routines are identical except that any double precision arguments of DRANDGAMMA are replaced in SRANDGAMMA by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDGAMMA** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: first parameter of the distribution.

Constraint:A>0.

#### — Input: DOUBLE PRECISION **B**

On input: second parameter of the distribution.

Constraint:B>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDGAMMA`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Gamma distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Gamma distribution CALL DRANDGAMMA(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDGAUSSIAN / DRANDGAUSSIAN

Generates a vector of random variates from a Gaussian distribution with probability density function, f(X), where: f(X) = [exp([-(X – m)^2] / [2s^2])] / [s * sqrt(2 * Pi)]. Here m is the mean, (*XMU*), and s^2 is the variance, (*VAR*) of the distribution.

*(Note that SRANDGAUSSIAN is the single precision version of DRANDGAUSSIAN. The argument lists of both routines are identical except that any double precision arguments of DRANDGAUSSIAN are replaced in SRANDGAUSSIAN by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDGAUSSIAN** (*N,XMU,VAR,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **XMU**

On input: mean of the distribution.

#### — Input: DOUBLE PRECISION **VAR**

On input: variance of the distribution.

Constraint:VAR>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDGAUSSIAN`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Gaussian distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION XMU,VAR DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) XMU,VAR C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Gaussian distribution CALL DRANDGAUSSIAN(N,XMU,VAR,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDGENERALDISCRETE / SRANDGENERALDISCRETE

Takes a reference vector initialized via one of `DRANDBINOMIALREFERENCE`

, `DRANDGEOMETRICREFERENCE`

, `DRANDHYPERGEOMETRICREFERENCE`

, `DRANDNEGATIVEBINOMIALREFERENCE`

, `DRANDPOISSONREFERENCE`

and generates a vector of random variates from it.

*(Note that SRANDGENERALDISCRETE is the single precision version of DRANDGENERALDISCRETE. The argument lists of both routines are identical except that any double precision arguments of DRANDGENERALDISCRETE are replaced in SRANDGENERALDISCRETE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDGENERALDISCRETE** (*N,REF,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **REF**(***)

On input: reference vector generated by one of the following: DRANDBINOMIALREFERENCE, DRANDGEOMETRICREFERENCE, DRANDHYPERGEOMETRICREFERENCE, DRANDNEGATIVEBINOMIALREFERENCE, DRANDPOISSONREFERENCE.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDGENERALDISCRETE`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Binomial distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER M DOUBLE PRECISION P INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDBINOMIALREFERENCE(M,P,REF,LREF,INFO) C Generate N variates from the Binomial distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDGEOMETRIC / SRANDGEOMETRIC

Generates a vector of random variates from a Geometric distribution with probability, f(X), defined by: f(X) = P * (1 – P)^X, X=0,1,….

*(Note that SRANDGEOMETRIC is the single precision version of DRANDGEOMETRIC. The argument lists of both routines are identical except that any double precision arguments of DRANDGEOMETRIC are replaced in SRANDGEOMETRIC by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDGEOMETRIC** (*N,P,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **P**

On input: distribution parameter.

Constraint: 0<=P<1.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDGEOMETRIC`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Geometric distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION P INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Geometric distribution CALL DRANDGEOMETRIC(N,P,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDGEOMETRICREFERENCE / SRANDGEOMETRICREFERENCE

Initializes a reference vector for use with `DRANDGENERALDISCRETE`

. Reference vector is for a Geometric distribution with probability, f(X), defined by: f(X) = P * (1 – P)^X, X=0,1,….

*(Note that SRANDGEOMETRICREFERENCE is the single precision version of DRANDGEOMETRICREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDGEOMETRICREFERENCE are replaced in SRANDGEOMETRICREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDGEOMETRICREFERENCE** (*P,REF,LREF,INFO*)

#### — Input: DOUBLE PRECISION **P**

On input: distribution parameter.

Constraint: 0<=P<1.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Geometric distribution using DRANDGENERALDISCRETE.

#### — Input/Output: INTEGER **LREF**

On input: either the length of the reference vector

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the Geometric distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION P INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDGEOMETRICREFERENCE(P,REF,LREF,INFO) C Generate N variates from the Geometric distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDHYPERGEOMETRIC / SRANDHYPERGEOMETRIC

Generates a vector of random variates from a Hypergeometric distribution with probability, f(X), defined by: f(X) = [s! * m! * (p – s)!] / [X! * (s – X)! * (m – X)! * (p – m – s + X)! * p!], if X = max(0,m+s-p),…,min(l,m), otherwise f(X) = 0. Here p is the size of the population, (*NP*), s is the size of the sample taken from the population, (*NS*) and m is the number of labeled, or specified, items in the population, (*M*).

*(Note that SRANDHYPERGEOMETRIC is the single precision version of DRANDHYPERGEOMETRIC. The argument lists of both routines are identical except that any double precision arguments of DRANDHYPERGEOMETRIC are replaced in SRANDHYPERGEOMETRIC by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDHYPERGEOMETRIC** (*N,NP,NS,M,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **NP**

On input: size of population.

Constraint:NP>=0.

#### — Input: INTEGER **NS**

On input: size of sample being taken from population.

Constraint: 0<=NS<=NP.

#### — Input: INTEGER **M**

On input: number of specified items in the population.

Constraint: 0<=M<=NP.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDHYPERGEOMETRIC`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Hypergeometric distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER NP,NS,M INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) NP,NS,M C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Hypergeometric distribution CALL DRANDHYPERGEOMETRIC(N,NP,NS,M,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDHYPERGEOMETRICREFERENCE / SRANDHYPERGEOMETRICREFERENCE

Initializes a reference vector for use with `DRANDGENERALDISCRETE`

. Reference vector is for a Hypergeometric distribution with probability, f(X), defined by: f(X) = [s! * m! * (p – s)!] / [X! * (s – X)! * (m – X)! * (p – m – s + X)! * p!], if X = max(0,m+s-p),…,min(l,m), otherwise f(X) = 0. Here p is the size of the population, (*NP*), s is the size of the sample taken from the population, (*NS*) and m is the number of labeled, or specified, items in the population, (*M*).

*(Note that SRANDHYPERGEOMETRICREFERENCE is the single precision version of DRANDHYPERGEOMETRICREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDHYPERGEOMETRICREFERENCE are replaced in SRANDHYPERGEOMETRICREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDHYPERGEOMETRICREFERENCE** (*NP,NS,M,REF,LREF,INFO*)

#### — Input: INTEGER **NP**

On input: size of population.

Constraint:NP>=0.

#### — Input: INTEGER **NS**

On input: size of sample being taken from population.

Constraint: 0<=NS<=NP.

#### — Input: INTEGER **M**

On input: number of specified items in the population.

Constraint: 0<=M<=NP.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Hypergeometric distribution using DRANDGENERALDISCRETE.

#### — Input/Output: INTEGER **LREF**

On input: either the length of the reference vector

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the Hypergeometric distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER NP, NS,M INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) NP, NS,M C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDHYPERGEOMETRICREFERENCE(NP, NS,M,REF,LREF,INFO) C Generate N variates from the Hypergeometric distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDINITIALIZE / SRANDINITIALIZE

Initialize one of the five supplied base generators; NAG basic generator, Wichmann-Hill generator, Mersenne Twister, L’Ecuyer’s combined recursive generator (MRG32k3a) or the Blum-Blum-Shub generator.

*(Note that SRANDINITIALIZE is the single precision version of DRANDINITIALIZE. The argument lists of both routines are identical except that any double precision arguments of DRANDINITIALIZE are replaced in SRANDINITIALIZE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDINITIALIZE** (*GENID,SUBID,SEED,LSEED,STATE, LSTATE,INFO*)

#### — Input: INTEGER **GENID**

On input: a numerical code indicating which of the five base generators to initialize.

- 1 = NAG basic generator (Basic NAG Generator).
- 2 = Wichmann-Hill generator (Wichmann-Hill Generator).
- 3 = Mersenne Twister (Mersenne Twister).
- 4 = L’Ecuyer’s Combined Recursive generator (L’Ecuyer’s Combined Recursive Generator).
- 5 = Blum-Blum-Shub generator (Blum-Blum-Shub Generator).
Constraint: 1<=

GENID<=5.

#### — Input: INTEGER **SUBID**

On input: if

GENID=2, thenSUBIDindicates which of the 273 Wichmann-Hill generators to use. IfGENID=5 thenSUBIDindicates the number of bits to use (v) from each of iteration of the Blum-Blum-Shub generator. In all other casesSUBIDis not referenced.

Constraint: IfGENID=2 then 1<=SUBID<=273 .

#### — Input: INTEGER **SEED**(*LSEED*)

On input: if

GENIDis not equal to 5 , thenSEEDis a vector of initial values for the base generator. These values must be positive integers. The number of values required depends on the base generator being used. The NAG basic generator requires one initial value, the Wichmann-Hill generator requires four initial values, the L’Ecuyer combined recursive generator requires six initial values and the Mersenne Twister requires 624 initial values. If the number of seeds required by the chosen generator is >LSEEDthenSEED(1) is used to initialize the NAG basic generator. This is then used to generate all of the remaining seed values required. In general it is best not to set all the elements ofSEEDto anything too obvious, such as a single repeated value or a simple sequence. Using such a seed array may lead to several similar values being created in a row when the generator is subsequently called. This is particularly true for the Mersenne Twister generator.In order to initialize the Blum-Blum-Shub generator two large prime values, p and q are required as well as an initial value s. As p, q and s can be of an arbitrary size, these values are factored into a expressed as a polynomial in B, where B=2^(24). For example, p can be factored into a polynomial of order l_p, with p = p_1 + p_2B + p_3B^2 + … + p_(l_p) * B^(l_p-1). The elements of

SEEDshould then be set to the following:

SEED(1) = l_pSEED(2) toSEED(l_p+1) = p_1 to p_(l_p)SEED(l_p+2) = l_qSEED(l_p+3) toSEED(l_p+l_q+2) = q_1 to q_(l_q)SEED(l_p+l_q+3) = l_sSEED(l_p+l_q+4) toSEED(l_p+l_q+l_s+3) = s_1 to s_(l_s)Constraint: If

GENIDis not equal to 5 thenSEED(i)>0 for i=1,2,…. IfGENID=5 thenSEEDmust take the values described above.

#### — Input/Output: INTEGER **LSEED**

On input: either the length of the seed vector,

SEED, or a value <=0 .

On output: ifLSEED<=0 on input, thenLSEEDis set to the number of initial values required by the selected generator, and the routine returns. OtherwiseLSEEDis left unchanged.

#### — Output: INTEGER **STATE**(*LSTATE*)

On output: the state vector required by all of the supplied distributional and base generators.

#### — Input/Output: INTEGER **LSTATE**

On input: either the length of the state vector,

STATE, or a value <=0 .

On output: ifLSTATE<=0 on input, thenLSTATEis set to the minimum length of the state vectorSTATEfor the base generator chosen, and the routine returns. OtherwiseLSTATEis left unchanged.

Constraint:LSTATE<=0 or the minimum length for the chosen base generator, given by:

GENID=1:LSTATE>=16,GENID=2:LSTATE>=20,GENID=3:LSTATE>=633,GENID=4:LSTATE>=61,GENID=5:LSTATE>=l_p+l_q+l_s+6 , where l_p, l_q and l_s are the order of the polynomials used to express the parameters p,q and s respectively.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, then either, or both ofLSEEDand / orLSTATEhave been set to the required length for vectorsSEEDandSTATErespectively. Of the two variablesLSEEDandLSTATE, only those which had an input value <=0 will have been set. TheSTATEvector will not have been initialized. IfINFO= 0 then the state vector,STATE, has been successfully initialized.

**Example:**

C Generate 100 values from the Beta distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Beta distribution CALL DRANDBETA(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDINITIALIZEBBS / SRANDINITIALIZEBBS

Alternative initialization routine for the Blum-Blum-Shub generator. Unlike the other base generators supplied with the ACML, the Blum-Blum-Shub generator requires two additional parameters, p and q as well as an initial state, s. As p, q and s can be of an arbitrary size. In order to avoid overflow these values are expressed as a polynomial in B, where B=2^(24). For example, p can be factored into a polynomial of order l_p, with p = p_1 + p_2B + p_3B^2 + … + p_(l_p) * B^(l_p-1), similarly q = q_1 + q_2B + q_3B^2 + … + q_(l_q) * B^(l_q-1) and s = s_1 + s_2B + s_3B^2 + \cdots + s_(l_s) * B^(l_s-1).

*(Note that SRANDINITIALIZEBBS is the single precision version of DRANDINITIALIZEBBS. The argument lists of both routines are identical except that any double precision arguments of DRANDINITIALIZEBBS are replaced in SRANDINITIALIZEBBS by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDINITIALIZEBBS** (*NBITS,LP,P,LQ,Q,LS,S,STATE,LSTATE, INFO*)

#### — Input: INTEGER **NBITS**

On input: the number of bits, v, to use from each iteration of the Blum-Blum-Shub generator. If

NBITS<1 thenNBITS=1. IfNBITS>15 thenNBITS=15.

#### — Input: INTEGER **LP**

On input: the order of the polynomial used to express p (l_p).

Constraint: 1 <=LP<=25.

#### — Input: INTEGER **P**(*LP*)

On input: the coefficients of the polynomial used to express p.

P(i)=p_i, i=1 to l_p.

Constraint: 0 <=P(i) < 2^(24)

#### — Input: INTEGER **LQ**

On input: the order of the polynomial used to express q (l_q).

Constraint: 1 <=LQ<=25.

#### — Input: INTEGER **Q**(*LQ*)

On input: the coefficients of the polynomial used to express q.

Q(i)=q_i, i=1 to l_q.

Constraint: 0 <=Q(i) < 2^(24)

#### — Input: INTEGER **LS**

On input: the order of the polynomial used to express s (l_s).

Constraint: 1 <=LS<=25.

#### — Input: INTEGER **S**(*LS*)

On input: the coefficients of the polynomial used to express s.

S(i)=s_i, i=1 to l_s.

Constraint: 0 <=S(i) < 2^(24)

#### — Output: INTEGER **STATE**(***)

On output: the initial state for the Blum-Blum-Shub generator with parameters

P,Q,SandNBITS.

#### — Input/Output: INTEGER **LSTATE**

On input: either the length of the state vector,

STATE, or a value <=0 .

On output: ifLSTATE<=0 on input, thenLSTATEis set to the minimum length of the state vectorSTATEfor the parameters chosen, and the routine returns. OtherwiseLSTATEis left unchanged.

Constraint:LSTATE<=0 orLSTATE>=l_p+l_q+l_s+6

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLSTATEhas been set to the required length for theSTATEvector. IfINFO= 0 then the state vector,STATE, has been successfully initialized.

## DRANDINITIALIZEUSER / SRANDINITIALIZEUSER

Registers a user supplied base generator so that it can be used with the ACML distributional generators.

*(Note that SRANDINITIALIZEUSER is the single precision version of DRANDINITIALIZEUSER. The argument lists of both routines are identical except that any double precision arguments of DRANDINITIALIZEUSER are replaced in SRANDINITIALIZEUSER by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDINITIALIZEUSER** (*UINI,UGEN,GENID,SUBID,SEED,LSEED, STATE,LSTATE,INFO*)

#### — Input: SUBROUTINE **UINI**

On input: routine that will be used to initialize the user supplied generator,

UGEN.

#### — Input: SUBROUTINE **UGEN**

On input: user supplied base generator.

#### — Input: INTEGER **GENID**

On input: parameter is passed directly to

UINI. Its function therefore depends on that routine.

#### — Input: INTEGER **SUBID**

On input: parameter is passed directly to

UINI. Its function therefore depends on that routine.

#### — Input: INTEGER **SEED**(*LSEED*)

On input: parameter is passed directly to

UINI. Its function therefore depends on that routine.

#### — Input/Output: INTEGER **LSEED**

On input: length of the vector

SEED. This parameter is passed directly toUINIand therefore its required value depends on that routine.

On output: whetherLSEEDchanges will depend onUINI.

#### — Output: INTEGER **STATE**(*LSTATE*)

On output: the state vector required by all of the supplied distributional generators. The value of

STATEreturned byUINIhas some housekeeping elements appended to the end before being returned by`DRANDINITIALIZEUSER`

. See User Supplied Generators for details about the form ofSTATE.

#### — Input/Output: INTEGER **LSTATE**

On input: length of the vector

STATE. This parameter is passed directly toUINIand therefore its required value depends on that routine.

On output: whetherLSTATEchanges will depend onUINI. IfLSTATE<=0 then it is assumed that a request for the required length ofSTATEhas been made. The value ofLSTATEreturned fromUINIis therefore adjusted to allow for housekeeping elements to be added to the end of theSTATEvector. This results in the value ofLSTATEreturned by`DRANDINITIALIZEUSER`

being 3 larger than that returned byUINI.

#### — Output: INTEGER **INFO**

On output:

INFOis an error indicator.`DRANDINITIALIZEUSER`

will return a value of -6 if the value ofLSTATEis between 1 and 3. OtherwiseINFOis passed directly back fromUINI. It is recommended that the value ofINFOreturned byUINIis kept consistent with the rest of the ACML, that is ifINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, then either, or both ofLSEEDand / orLSTATEhave been set to the required length for vectorsSEEDandSTATErespectively and theSTATEvector has not have been initialized. IfINFO= 0 then the state vector,STATE, has been successfully initialized.

**Example:**

C Generate 100 values from the Uniform distribution using C a user supplied base generator INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,NSKIP,SEED(1),STATE(LSTATE) INTEGER X(N) DOUBLE PRECISION A,B C Set the seed SEED(1) = 1234 C Set the distributional parameters A = 0.0D0 B = 1.0D0 C Initialize the base generator. Here ACMLRNGNB0GND is a user C supplied generator and ACMLRNGNB0INI its initializer CALL DRANDINITIALIZEUSER(ACMLRNGNB0INI,ACMLRNGNB0GND,1,0,SEED, * LSEED,STATE,LSTATE,INFO) C Generate N variates from the Univariate distribution CALL DRANDUNIFORM(N,A,B,STATE,X,LDX,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDLEAPFROG / SRANDLEAPFROG

Amend a generator so that it will generate every *K*th value.

*(Note that SRANDLEAPFROG is the single precision version of DRANDLEAPFROG. The argument lists of both routines are identical except that any double precision arguments of DRANDLEAPFROG are replaced in SRANDLEAPFROG by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDLEAPFROG** (*N,K,STATE,INFO*)

#### — Input: INTEGER **N**

On input: total number of streams being used.

Constraint:N>0.

#### — Input: INTEGER **K**

On input: number of the current stream

Constraint: 0<K<=N.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDLEAPFROG`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: TheSTATEvector for a generator that has been advancedK-1 places and will return everyNth value.

Constraint: TheSTATEarray must be for either the NAG basic, Wichmann-Hill or L’Ecuyer Combined Recursive base generators.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 3 * 100 values from the Uniform distribution C Multiple streams generated using the Leap Frog method INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO INTEGER SEED(1),STATE1(LSTATE),STATE2(LSTATE),STATE3(LSTATE) INTEGER X1(N),X2(N),X3(N) DOUBLE PRECISION A,B C Set the seed SEED(1) = 1234 C Set the distributional parameters A = 0.0D0 B = 1.0D0 C Initialize the STATE1 vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE1,LSTATE,INFO) C Copy the STATE1 vector into other state vectors DO 20 I = 1,LSTATE STATE2(I) = STATE1(I) STATE3(I) = STATE1(I) 20 CONTINUE C Update each stream so they generate every 3rd value CALL DRANDLEAPFROG(3,1,STATE1,INFO) CALL DRANDLEAPFROG(3,2,STATE2,INFO) CALL DRANDLEAPFROG(3,3,STATE3,INFO) C Generate 3 sets of N variates from the Univariate distribution CALL DRANDUNIFORM(N,A,B,STATE1,X1,LDX,INFO) CALL DRANDUNIFORM(N,A,B,STATE2,X2,LDX,INFO) CALL DRANDUNIFORM(N,A,B,STATE3,X3,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) X1(I),X2(I),X3(I) 40 CONTINUE

## DRANDLOGISTIC / SRANDLOGISTIC

Generates a vector of random variates from a logistic distribution with probability density function, f(X), where: f(X) = [exp((X – A) / B)] / [B * (1 + exp((X – A) / B))^2].

*(Note that SRANDLOGISTIC is the single precision version of DRANDLOGISTIC. The argument lists of both routines are identical except that any double precision arguments of DRANDLOGISTIC are replaced in SRANDLOGISTIC by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDLOGISTIC** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: mean of the distribution.

#### — Input: DOUBLE PRECISION **B**

On input: spread of the distribution.

B= SQRT(3) * SIGMA / PI, where SIGMA is the standard deviation of the distribution.

Constraint:B>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDLOGISTIC`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Logistic distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Logistic distribution CALL DRANDLOGISTIC(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDLOGNORMAL / SRANDLOGNORMAL

Generates a vector of random variates from a lognormal distribution with probability density function, f(X), where: f(X) = [exp([-(log(X) – m)^2] / [2s^2])] / [X * s * sqrt(2 * Pi)], if X > 0, otherwise f(X) = 0. Here m is the mean, (*XMU*), and s^2 is the variance, (*VAR*) of the underlying Gaussian distribution.

*(Note that SRANDLOGNORMAL is the single precision version of DRANDLOGNORMAL. The argument lists of both routines are identical except that any double precision arguments of DRANDLOGNORMAL are replaced in SRANDLOGNORMAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDLOGNORMAL** (*N,XMU,VAR,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **XMU**

On input: mean of the underlying Gaussian distribution.

#### — Input: DOUBLE PRECISION **VAR**

On input: variance of the underlying Gaussian distribution.

Constraint:VAR>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDLOGNORMAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Lognormal distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION XMU,VAR DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) XMU,VAR C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Lognormal distribution CALL DRANDLOGNORMAL(N,XMU,VAR,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDMULTINOMIAL / SRANDMULTINOMIAL

Generates a matrix of random variates from a Multinomial distribution with probability, f(X), defined by: f(X) = [M! * P1^(X1) * P2^(X2) * … * PK^(XK)] / [X1! * X2! * … * XK!], where X = (X1,X2,…,XK), P = (P1,P2,…,PK), X1 + X2 + … + XK = 1 and P1 + P2 + … + PK = 1.

*(Note that SRANDMULTINOMIAL is the single precision version of DRANDMULTINOMIAL. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTINOMIAL are replaced in SRANDMULTINOMIAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTINOMIAL** (*N,M,P,K,STATE,X,LDX,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of trials.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **P**(*K*)

On input: vector of probabilities for each of the

Kpossible outcomes.

Constraint: 0 <= P <= 1, for all elements of P and P1 + P2 + … + PK = 1.

#### — Input: INTEGER **K**

On input: number of possible outcomes.

Constraint:K>=2.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDBINOMIAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*LDX,K*)

On output: matrix of variates from the specified distribution.

#### — Input: INTEGER **LDX**

On input: leading dimension of

Xin the calling routine.

Constraint:LDX>=N.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Multinomial distribution INTEGER LSTATE,N, MK PARAMETER (LSTATE=16,N=100,MK=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,K,M INTEGER X(N,MK) DOUBLE PRECISION P(MK) C Set array sizes LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) K READ(5,*) (P(I),I=1,K) C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Multinomial distribution CALL DRANDMULTINOMIAL(N,M,P,K,STATE,X,LDX,INFO) C Print the results DO 20 I = 1,N WRITE(6,*) (X(I,J),J=1,K) 20 CONTINUE

## DRANDMULTINORMAL / SRANDMULTINORMAL

Generates an array of random variates from a Multivariate Normal distribution with probability density function, f(X), where: f(X) = sqrt([det(C^(-1))] / [(2 * Pi)^M]) * exp(-(X – u)^T * C^(-1) * (X – u)), where u is the vector of means, *XMU*.

*(Note that SRANDMULTINORMAL is the single precision version of DRANDMULTINORMAL. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTINORMAL are replaced in SRANDMULTINORMAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTINORMAL** (*N,M,XMU,C,LDC,STATE,X,LDX,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of dimensions for the distribution.

Constraint:M>=1.

#### — Input: DOUBLE PRECISION **XMU**(*M*)

On input: vector of means for the distribution.

#### — Input: DOUBLE PRECISION **C**(*LDC,M*)

On input: variance / covariance matrix for the distribution.

#### — Input: INTEGER **LDC**

On input: leading dimension of

Cin the calling routine.

Constraint:LDC>=M.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDMULTINORMAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*LDX,M*)

On output: matrix of variates from the specified distribution.

#### — Input: INTEGER **LDX**

On input: leading dimension of

Xin the calling routine.

Constraint:LDX>=N.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the C Multivariate Normal distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the C Multivariate Normal distribution CALL DRANDMULTINORMAL(N,M,XMU,C,LDC,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDMULTINORMALR / SRANDMULTINORMALR

Generates an array of random variates from a Multivariate Normal distribution using a reference vector initialized by `DRANDMULTINORMALREFERENCE`

.

*(Note that SRANDMULTINORMALR is the single precision version of DRANDMULTINORMALR. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTINORMALR are replaced in SRANDMULTINORMALR by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTINORMALR** (*N,REF,STATE,X,LDX,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **REF**(***)

On input: a reference vector generated by DRANDMULTINORMALREFERENCE.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDMULTINORMALR`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*LDX,M*)

On output: matrix of variates from the specified distribution.

#### — Input: INTEGER **LDX**

On input: leading dimension of

Xin the calling routine.

Constraint:LDX>=N.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the C Multivariate Normal distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) INTEGER LREF DOUBLE PRECISION REF(1000) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDMULTINORMALREFERENCE(M,XMU,C,LDC,REF,LREF,INFO) C Generate N variates from the C Multivariate Normal distribution CALL DRANDMULTINORMALR(N,REF,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDMULTINORMALREFERENCE / SRANDMULTINORMALREFERENCE

Initializes a reference vector for use with `DRANDMULTINORMALR`

. Reference vector is for a Multivariate Normal distribution with probability density function, f(X), where: f(X) = sqrt([det(C^(-1))] / [(2 * Pi)^M]) * exp(-(X – u)^T * C^(-1) * (X – u)), where u is the vector of means, *XMU*.

*(Note that SRANDMULTINORMALREFERENCE is the single precision version of DRANDMULTINORMALREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTINORMALREFERENCE are replaced in SRANDMULTINORMALREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTINORMALREFERENCE** (*M,XMU,C,LDC,REF,LREF,INFO*)

#### — Input: INTEGER **M**

On input: number of dimensions for the distribution.

Constraint:M>=1.

#### — Input: DOUBLE PRECISION **XMU**(*M*)

On input: vector of means for the distribution.

#### — Input: DOUBLE PRECISION **C**(*LDC,M*)

On input: variance / covariance matrix for the distribution.

#### — Input: INTEGER **LDC**

On input: leading dimension of

Cin the calling routine.

Constraint:LDC>=M.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Multivariate Normal distribution using DRANDMULTINORMALR.

#### — Input/Output: INTEGER **LREF**

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the C Multivariate Normal distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) INTEGER LREF DOUBLE PRECISION REF(1000) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDMULTINORMALREFERENCE(M,XMU,C,LDC,REF,LREF,INFO) C Generate N variates from the C Multivariate Normal distribution CALL DRANDMULTINORMALR(N,REF,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDMULTISTUDENTST / SRANDMULTISTUDENTST

Generates an array of random variates from a Multivariate Students T distribution with probability density function, f(X), where: f(X) = [Gamma([v + M] / 2)] / [(Pi * v)^(M / 2) * Gamma(v / 2) * det(C)^(1 / 2)] * (1 + [(X – u)^T * C^(-1) * (X – u)] / v)^(-[v + M] / 2), where u us the vector of means, *XMU* and v is the degrees of freedom, *DF*.

*(Note that SRANDMULTISTUDENTST is the single precision version of DRANDMULTISTUDENTST. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTISTUDENTST are replaced in SRANDMULTISTUDENTST by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTISTUDENTST** (*N,M,DF,XMU,C,LDC,STATE,X,LDX,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of dimensions for the distribution.

Constraint:M>=1.

#### — Input: INTEGER **DF**

On input: degrees of freedom.

Constraint:DF>2.

#### — Input: DOUBLE PRECISION **XMU**(*M*)

On input: vector of means for the distribution.

#### — Input: DOUBLE PRECISION **C**(*LDC,M*)

On input: matrix defining the variance / covariance for the distribution. The variance / covariance matrix is given by [DF * C]/ [DF – 2].

#### — Input: INTEGER **LDC**

On input: leading dimension of

Cin the calling routine.

Constraint:LDC>=M.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDMULTISTUDENTST`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*LDX,M*)

On output: matrix of variates from the specified distribution.

#### — Input: INTEGER **LDX**

On input: leading dimension of

Xin the calling routine.

Constraint:LDX>=N.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the C Multivariate Students T distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M,DF DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,DF READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the C Multivariate Students T distribution CALL DRANDMULTISTUDENTST(N,M,DF,XMU,C,LDC,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDMULTISTUDENTSTR / SRANDMULTISTUDENTSTR

Generates an array of random variates from a Multivariate Students T distribution using a reference vector initialized by `DRANDMULTISTUDENTSTREFERENCE`

.

*(Note that SRANDMULTISTUDENTSTR is the single precision version of DRANDMULTISTUDENTSTR. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTISTUDENTSTR are replaced in SRANDMULTISTUDENTSTR by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTISTUDENTSTR** (*N,REF,STATE,X,LDX,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **REF**(***)

On input: a reference vector generated by DRANDMULTISTUDENTSTREFERENCE.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDMULTISTUDENTSTR`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*LDX,M*)

On output: matrix of variates from the specified distribution.

#### — Input: INTEGER **LDX**

On input: leading dimension of

Xin the calling routine.

Constraint:LDX>=N.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the C Multivariate Students T distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M,DF DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) INTEGER LREF DOUBLE PRECISION REF(1000) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,DF READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDMULTISTUDENTSTREFERENCE(M,DF,XMU,C,LDC,REF,LREF,INFO) C Generate N variates from the C Multivariate Students T distribution CALL DRANDMULTISTUDENTSTR(N,REF,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDMULTISTUDENTSTREFERENCE / SRANDMULTISTUDENTSTREFERENCE

Initializes a reference vector for use with `DRANDMULTISTUDENTSTR`

. Reference vector is for a Multivariate Students T distribution with probability density function, f(X), where: f(X) = [Gamma([v + M] / 2)] / [(Pi * v)^(M / 2) * Gamma(v / 2) * det(C)^(1 / 2)] * (1 + [(X – u)^T * C^(-1) * (X – u)] / v)^(-[v + M] / 2), where u us the vector of means, *XMU* and v is the degrees of freedom, *DF*.

*(Note that SRANDMULTISTUDENTSTREFERENCE is the single precision version of DRANDMULTISTUDENTSTREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDMULTISTUDENTSTREFERENCE are replaced in SRANDMULTISTUDENTSTREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDMULTISTUDENTSREFERENCE** (*M,DF,XMU,C,LDC,REF,LREF,INFO*)

#### — Input: INTEGER **M**

On input: number of dimensions for the distribution.

Constraint:M>=1.

#### — Input: INTEGER **DF**

On input: degrees of freedom.

Constraint:DF>2.

#### — Input: DOUBLE PRECISION **XMU**(*M*)

On input: vector of means for the distribution.

#### — Input: DOUBLE PRECISION **C**(*LDC,M*)

On input: matrix defining the variance / covariance for the distribution. The variance / covariance matrix is given by [DF * C]/ [DF – 2].

#### — Input: INTEGER **LDC**

On input: leading dimension of

Cin the calling routine.

Constraint:LDC>=M.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Multivariate Students T distribution using DRANDMULTISTUDENTSTR.

#### — Input/Output: INTEGER **LREF**

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the C Multivariate Students T distribution INTEGER LSTATE,N, MM PARAMETER (LSTATE=16,N=100,MM=10) INTEGER I,J,INFO,SEED(1),STATE(LSTATE) INTEGER LDC,LDX,M,DF DOUBLE PRECISION X(N,MM),XMU(MM),C(MM,MM) INTEGER LREF DOUBLE PRECISION REF(1000) C Set array sizes LDC = MM LDX = N C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,DF READ(5,*) (XMU(I),I=1,M) DO 20 I = 1,M READ(5,*) (C(I,J),J=1,M) 20 CONTINUE C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDMULTISTUDENTSTREFERENCE(M,DF,XMU,C,LDC,REF,LREF,INFO) C Generate N variates from the C Multivariate Students T distribution CALL DRANDMULTISTUDENTSTR(N,REF,STATE,X,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) (X(I,J),J=1,M) 40 CONTINUE

## DRANDNEGATIVEBINOMIAL / SRANDNEGATIVEBINOMIAL

Generates a vector of random variates from a Negative Binomial distribution with probability f(X) defined by: f(X) = [(M + X – 1)! * P^X * (1 – P)^M] / [X! * (M – 1)!], X=0,1,…

*(Note that SRANDNEGATIVEBINOMIAL is the single precision version of DRANDNEGATIVEBINOMIAL. The argument lists of both routines are identical except that any double precision arguments of DRANDNEGATIVEBINOMIAL are replaced in SRANDNEGATIVEBINOMIAL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDNEGATIVEBINOMIAL** (*N,M,P,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of failures.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **P**

On input: probability of success.

Constraint: 0<=P<1.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDNEGATIVEBINOMIAL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Negative Binomial distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER M DOUBLE PRECISION P INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Negative Binomial distribution CALL DRANDNEGATIVEBINOMIAL(N,M,P,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDNEGATIVEBINOMIALREFERENCE / SRANDNEGATIVEBINOMIALREFERENCE

Initializes a reference vector for use with `DRANDGENERALDISCRETE`

. Reference vector is for a Negative Binomial distribution with probability f(X) defined by: f(X) = [(M + X – 1)! * P^X * (1 – P)^M] / [X! * (M – 1)!], X=0,1,…

*(Note that SRANDNEGATIVEBINOMIALREFERENCE is the single precision version of DRANDNEGATIVEBINOMIALREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDNEGATIVEBINOMIALREFERENCE are replaced in SRANDNEGATIVEBINOMIALREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDNEGATIVEBINOMIALREFERENCE** (*M,P,REF,LREF,INFO*)

#### — Input: INTEGER **M**

On input: number of failures.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **P**

>

On input: probability of success.

Constraint: 0<=P<1.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Negative Binomial distribution using DRANDGENERALDISCRETE.

#### — Input/Output: INTEGER **LREF**

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the Negative Binomial distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER M DOUBLE PRECISION P INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) M,P C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDNEGATIVEBINOMIALREFERENCE(M,P,REF,LREF,INFO) C Generate N variates from the Negative Binomial distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDPOISSON / SRANDPOISSON

Generates a vector of random variates from a Poisson distribution with probability f(X) defined by: f(X) = [l^X * exp(-l)] / X!, X=0,1,…, where l is the mean of the distribution, *LAMBDA*.

*(Note that SRANDPOISSON is the single precision version of DRANDPOISSON. The argument lists of both routines are identical except that any double precision arguments of DRANDPOISSON are replaced in SRANDPOISSON by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDPOISSON** (*N,LAMBDA,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **M**

On input: number of failures.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **LAMBDA**

On input: mean of the distribution.

Constraint:LAMBDA>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDPOISSON`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: INTEGER **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Poisson distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION LAMBDA INTEGER X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) LAMBDA C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Poisson distribution CALL DRANDPOISSON(N,LAMBDA,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDPOISSONREFERENCE / SRANDPOISSONREFERENCE

Initializes a reference vector for use with `DRANDGENERALDISCRETE`

. Reference vector is for a Poisson distribution with probability f(X) defined by: f(X) = [l^X * exp(-l)] / X!, X=0,1,…, where l is the mean of the distribution, *LAMBDA*.

*(Note that SRANDPOISSONREFERENCE is the single precision version of DRANDPOISSONREFERENCE. The argument lists of both routines are identical except that any double precision arguments of DRANDPOISSONREFERENCE are replaced in SRANDPOISSONREFERENCE by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDPOISSONREFERENCE** (*LAMBDA,REF,LREF,INFO*)

#### — Input: INTEGER **M**

On input: number of failures.

Constraint:M>=0.

#### — Input: DOUBLE PRECISION **LAMBDA**

On input: mean of the distribution.

Constraint:LAMBDA>=0.

#### — Output: DOUBLE PRECISION **REF**(*LREF*)

On output: if

INFOreturns with a value of 0 thenREFcontains reference information required to generate values from a Poisson distribution using DRANDGENERALDISCRETE.

#### — Input/Output: INTEGER **LREF**

REF, or -1.

On output: ifLREF=-1 on input, thenLREFis set to the recommended length of the reference vector and the routine returns. OtherwiseLREFis left unchanged.

#### — Output: INTEGER **INFO**

INFOis an error indicator. IfINFO= -i on exit, the i-th argument had an illegal value. IfINFO=1 on exit, thenLREFhas been set to the recommended length for the reference vectorREF. IfINFO= 0 then the reference vector,REF, has been successfully initialized.

**Example:**

C Generate 100 values from the Poisson distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION LAMBDA INTEGER X(N) INTEGER LREF DOUBLE PRECISION REF(1000) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) LAMBDA C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Initialize the reference vector LREF = 1000 CALL DRANDPOISSONREFERENCE(LAMBDA,REF,LREF,INFO) C Generate N variates from the Poisson distribution CALL DRANDGENERALDISCRETE(N,REF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDSKIPAHEAD / SRANDSKIPAHEAD

Advance a generator *N* places.

*(Note that SRANDSKIPAHEAD is the single precision version of DRANDSKIPAHEAD. The argument lists of both routines are identical except that any double precision arguments of DRANDSKIPAHEAD are replaced in SRANDSKIPAHEAD by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDSKIPAHEAD** (*N,STATE,INFO*)

#### — Input: INTEGER **N**

On input: number of places to skip ahead.

Constraint:N>=0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDSKIPAHEAD`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: TheSTATEvector for a generator that has been advancedNplaces.

Constraint: TheSTATEvector must be for either the NAG basic, Wichmann-Hill or L’Ecuyer Combined Recursive base generators.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 3 * 100 values from the Uniform distribution C Multiple streams generated using the Skip Ahead method INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,NSKIP INTEGER SEED(1),STATE1(LSTATE),STATE2(LSTATE),STATE3(LSTATE) INTEGER X1(N),X2(N),X3(N) DOUBLE PRECISION A,B C Set the seed SEED(1) = 1234 C Set the distributional parameters A = 0.0D0 B = 1.0D0 C Initialize the STATE1 vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE1,LSTATE,INFO) C Copy the STATE1 vector into other state vectors DO 20 I = 1,LSTATE STATE2(I) = STATE1(I) STATE3(I) = STATE1(I) 20 CONTINUE C Calculate how many places we want to skip, this C should be >> than the number of variates we C wish to generate from each stream NSKIP = N * N C Advance each stream, first does not need changing CALL DRANDSKIPAHEAD(NSKIP,STATE2,INFO) CALL DRANDSKIPAHEAD(2*NSKIP,STATE3,INFO) C Generate 3 sets of N variates from the Univariate distribution CALL DRANDUNIFORM(N,A,B,STATE1,X1,LDX,INFO) CALL DRANDUNIFORM(N,A,B,STATE2,X2,LDX,INFO) CALL DRANDUNIFORM(N,A,B,STATE3,X3,LDX,INFO) C Print the results DO 40 I = 1,N WRITE(6,*) X1(I),X2(I),X3(I) 40 CONTINUE

## DRANDSTUDENTST / SRANDSTUDENTST

Generates a vector of random variates from a Students T distribution with probability density function, f(X), where: f(X) = [((v – 1) / 2)!] / [(v / 2 – 1)! * sqrt(Pi * v) * (1 + X^2 / v)^((v + 1) / 2). Here v is the degrees of freedom, *DF*.

*(Note that SRANDSTUDENTST is the single precision version of DRANDSTUDENTST. The argument lists of both routines are identical except that any double precision arguments of DRANDSTUDENTST are replaced in SRANDSTUDENTST by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDSTUDENTST** (*N,DF,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: INTEGER **DF**

On input: degrees of freedom.

Constraint:DF>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDSTUDENTST`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Students T distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) INTEGER DF DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) DF C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Students T distribution CALL DRANDSTUDENTST(N,DF,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDTRIANGULAR / SRANDTRIANGULAR

Generates a vector of random variates from a Triangular distribution with probability density function, f(X), where: f(X) = [2 * (X – XMIN)] / [(XMAX – XMIN) * (XMED – XMIN)], if XMIN < X <= XMED, else f(X) = [2 * (XMAX – X)] / [(XMAX – XMIN) * (XMAX – XMED)], if XMED < X <= XMAX, otherwise f(X) = 0.

*(Note that SRANDTRIANGULAR is the single precision version of DRANDTRIANGULAR. The argument lists of both routines are identical except that any double precision arguments of DRANDTRIANGULAR are replaced in SRANDTRIANGULAR by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDTRIANGULAR** (*N,XMIN,XMED,XMAX,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **XMIN**

On input: minimum value for the distribution.

#### — Input: DOUBLE PRECISION **XMED**

On input: median value for the distribution.

Constraint:XMIN<=XMED<=XMAX.

#### — Input: DOUBLE PRECISION **XMAX**

On input: maximum value for the distribution.

Constraint:XMAX>=XMIN.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDTRIANGULAR`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Triangular distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION XMIN,XMAX,XMED DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) XMIN,XMAX,XMED C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Triangular distribution CALL DRANDTRIANGULAR(N,XMIN,XMAX,XMED,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDUNIFORM / SRANDUNIFORM

Generates a vector of random variates from a Uniform distribution with probability density function, f(X), where: f(X) = 1 / (B – A).

*(Note that SRANDUNIFORM is the single precision version of DRANDUNIFORM. The argument lists of both routines are identical except that any double precision arguments of DRANDUNIFORM are replaced in SRANDUNIFORM by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDUNIFORM** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: minimum value for the distribution.

#### — Input: DOUBLE PRECISION **B**

On input: maximum value for the distribution.

Constraint:B>=A.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDUNIFORM`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Uniform distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Uniform distribution CALL DRANDUNIFORM(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDVONMISES / SRANDVONMISES

Generates a vector of random variates from a Von Mises distribution with probability density function, f(X), where: f(X) = [exp(k * cos(X))] / [2 * Pi * I0(k)], where X is reduced modulo 2 * Pi, so that it lies between +/- Pi, and k is the concentration parameter *VK*.

*(Note that SRANDVONMISES is the single precision version of DRANDVONMISES. The argument lists of both routines are identical except that any double precision arguments of DRANDVONMISES are replaced in SRANDVONMISES by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDVONMISES** (*N,VK,,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **VK**

On input: concentration parameter.

Constraint:VK>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDVONMISES`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Von Mises distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION VK DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) VK C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Von Mises distribution CALL DRANDVONMISES(N,VK,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DRANDWEIBULL / SRANDWEIBULL

Generates a vector of random variates from a Weibull distribution with probability density function, f(X), where: f(X) = [A * X^(A – 1) * exp(-(X^A) / B)] / B, if X > 0, otherwise f(X)=0.

*(Note that SRANDWEIBULL is the single precision version of DRANDWEIBULL. The argument lists of both routines are identical except that any double precision arguments of DRANDWEIBULL are replaced in SRANDWEIBULL by single precision arguments – type REAL in FORTRAN or type float in C).*

### SUBROUTINE:**DRANDWEIBULL** (*N,A,B,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: number of variates required.

Constraint:N>=0.

#### — Input: DOUBLE PRECISION **A**

On input: shape parameter for the distribution.

Constraint:A>0.

#### — Input: DOUBLE PRECISION **B**

On input: scale parameter for the distribution.

Constraint:B>0.

#### — Input/Output: INTEGER **STATE**(***)

The

STATEvector holds information on the state of the base generator being used and as such its minimum length varies. Prior to calling`DRANDWEIBULL`

STATEmust have been initialized. See Initialization of the Base Generators for information on initialization of theSTATEvariable.

On input: the current state of the base generator.

On output: the updated state of the base generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: vector of variates from the specified distribution.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Generate 100 values from the Weibull distribution INTEGER LSTATE,N PARAMETER (LSTATE=16,N=100) INTEGER I,INFO,SEED(1),STATE(LSTATE) DOUBLE PRECISION A,B DOUBLE PRECISION X(N) C Set the seed SEED(1) = 1234 C Read in the distributional parameters READ(5,*) A,B C Initialize the STATE vector CALL DRANDINITIALIZE(1,1,SEED,1,STATE,LSTATE,INFO) C Generate N variates from the Weibull distribution CALL DRANDWEIBULL(N,A,B,STATE,X,INFO) C Print the results WRITE(6,*) (X(I),I=1,N)

## DZFFT Routine Documentation

### SUBROUTINE:**DZFFT** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`DZFFT`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real transform is performed. Initializations are assumed to have been performed by a prior call to`DZFFT`

.MODE=2 : (default) initializations and a real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis the length of the real sequenceX

#### — Input/Output: DOUBLE PRECISION **X**(*N*)

On input:

Xcontains the real sequence of lengthNto be transformed.

On output:Xcontains the transformed Hermitian sequence.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT(0,N,X,COMM,INFO) CALL DZFFT(1,N,X,COMM,INFO) DO 10 I = N/2+2, N X(I) = -X(I) 10 CONTINUE CALL ZDFFT(2,N,X,COMM,INFO)

## DZFFT1D Routine Documentation

### SUBROUTINE:**DZFFT1D** (*MODE,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`DZFFT1D`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`DZFFT1D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inX

#### — Input: DOUBLE PRECISION **X**(*N*)

On input:

Xcontains the real sequence of lengthNto be transformed.

#### — Output: DOUBLE COMPLEX **Y**(*N/2+1*)

On output:

Ycontains the tansformed complex sequence with roughly half redundant information due to complex conjugate removed.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT1D(0,N,X,Y,COMM,INFO) CALL DZFFT1D(1,N,X,Y,COMM,INFO) DO 10 I = 1, N/2+1 Y(I) = -Y(I)*EXP(-DBLE(I-1)/DBLE(N)) 10 CONTINUE CALL ZDFFT1D(2,N,Y,X,COMM,INFO)

## DZFFT1M Routine Documentation

### SUBROUTINE:**DZFFT1M** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`DZFFT1M`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`DZFFT1M`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inX

#### — Input: DOUBLE PRECISION **X**(*N*M*)

On input:

Xcontains theMreal sequences of lengthNto be transformed. Elementiof sequencejis stored in locationi + (j – 1) * NofX.

#### — Output: DOUBLE COMPLEX **Y**((*N/2+1*)**M*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT1M(0,M,N,X,Y,COMM,INFO) CALL DZFFT1M(1,M,N,X,Y,COMM,INFO) DO J = 1, M DO I = 1, N/2+1 Y(I,J) = -Y(I,J)*EXP(-DBLE(I-1)/DBLE(N)) END DO EMD DO CALL ZDFFT1M(2,M,N,Y,X,COMM,INFO)

## DZFFT2D Routine Documentation

### SUBROUTINE:**DZFFT2D** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`DZFFT2D`

.

On input:

MODE=0 : only default initializations (specific toMandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`DZFFT2D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of rows in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the first dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of data to be transformed. IfXis declared as a 2D array thenNis the second dimension ofX.

#### — Input: DOUBLE PRECISION **X**(*M*N*)

On input:

Xcontains theMbyNreal 2D array to be transformed. Elementijis stored in locationi + (j – 1) * MofX.

#### — Output: DOUBLE COMPLEX **Y**((*M/2+1*)**N*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*M+6*N+300*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthM*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT2D(0,M,N,X,Y,COMM,INFO) CALL DZFFT2D(1,M,N,X,Y,COMM,INFO) DO J = 1, N DO I = 1, M/2+1 Y(I,J) = -Y(I,J)/SQRT(DBLE(M*N)) END DO EMD DO CALL ZDFFT2D(2,M,N,Y,X,COMM,INFO)

## DZFFT3D Routine Documentation

### SUBROUTINE:**DZFFT3D** (*MODE,L,M,N,X,Y,COMM,INFO*)

**MODE**

The value of

MODEon input determines the operation performed by`DZFFT3D`

.

On input:

MODE=0 : only default initializations (specific toL,MandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`DZFFT3D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **L**

On input:

Lis the length of the first dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the first dimension ofX.

#### — Input: INTEGER **M**

On input:

Mis the length of the second dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the second dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the length of the third dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the third dimension ofX.

#### — Input: DOUBLE PRECISION **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNreal 3D array to be transformed. Elementijkis stored in locationi + (j – 1)* L + (k – 1) * L * MofX.

#### — Output: DOUBLE COMPLEX **Y**((*L/2+1*)**M*N*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*L+6*M+6*N+500*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthL*M*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT3D(0,L,M,N,X,Y,COMM,INFO) CALL DZFFT3D(1,L,M,N,X,Y,COMM,INFO) DO K = 1, N DO J = 1, M DO I = 1, L/2+1 Y(I,J,K) = -Y(I,J,K)/SQRT(DBLE(L*M*N)) END DO END DO EMD DO CALL ZDFFT3D(2,L,M,N,Y,X,COMM,INFO)

## DZFFTM Routine Documentation

### SUBROUTINE:**DZFFTM** (*M,N,X,COMM,INFO*)

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the real sequences inX

#### — Input/Output: DOUBLE PRECISION **X**(*N*M*)

On input:

Xcontains theMreal sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.

On output:Xcontains the transformed Hermitian sequences.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFTM(1,N,X,COMM,INFO) CALL DZFFTM(2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*X(N-I+1,2) 10 CONTINUE CALL ZDFFTM(2,N,X(1,3),COMM,INFO)

## FORTRAN specifications:

### SUBROUTINE:**ACMLVERSION** (*MAJOR, MINOR, PATCH*)

#### — INTEGER: *MAJOR, MINOR, PATCH*

### SUBROUTINE:**ACMLINFO** ()

## SCFFT Routine Documentation

### SUBROUTINE:**SCFFT** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`SCFFT`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real transform is performed. Initializations are assumed to have been performed by a prior call to`SCFFT`

.MODE=2 : (default) initializations and a real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis the length of the real sequenceX

#### — Input/Output: REAL **X**(*N*)

On input:

Xcontains the real sequence of lengthNto be transformed.

On output:Xcontains the transformed Hermitian sequence.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT(0,N,X,COMM,INFO) CALL SCFFT(1,N,X,COMM,INFO) DO 10 I = N/2+2, N X(I) = -X(I) 10 CONTINUE CALL CSFFT(2,N,X,COMM,INFO)

## SCFFT1D Routine Documentation

### SUBROUTINE:**SCFFT1D** (*MODE,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`SCFFT1D`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`SCFFT1D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inX

#### — Input: REAL **X**(*N*)

On input:

Xcontains the real sequence of lengthNto be transformed.

#### — Output: COMPLEX **Y**(*N/2+1*)

On output:

Ycontains the tansformed complex sequence with roughly half redundant information due to complex conjugate removed.

#### — Input/Output: REAL **COMM**(*4*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT1D(0,N,X,Y,COMM,INFO) CALL SCFFT1D(1,N,X,Y,COMM,INFO) DO 10 I = 1, N/2+1 Y(I) = -Y(I)*EXP(-REAL(I-1)/REAL(N)) 10 CONTINUE CALL CSFFT1D(2,N,Y,X,COMM,INFO)

## SCFFT1M Routine Documentation

### SUBROUTINE:**SCFFT1M** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`SCFFT1M`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`SCFFT1M`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inX

#### — Input: REAL **X**(*N*M*)

On input:

Xcontains theMreal sequences of lengthNto be transformed. Elementiof sequencejis stored in locationi + (j – 1) * NofX.

#### — Output: COMPLEX **Y**((*N/2+1*)**M*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: REAL **COMM**(*4*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT1M(0,M,N,X,Y,COMM,INFO) CALL SCFFT1M(1,M,N,X,Y,COMM,INFO) DO J = 1, M DO I = 1, N/2+1 Y(I,J) = -Y(I,J)*EXP(-REAL(I-1)/REAL(N)) END DO EMD DO CALL CSFFT1M(2,M,N,Y,X,COMM,INFO)

## SCFFT2D Routine Documentation

### SUBROUTINE:**SCFFT2D** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`SCFFT2D`

.

On input:

MODE=0 : only default initializations (specific toMandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`SCFFT2D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of rows in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the first dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of data to be transformed. IfXis declared as a 2D array thenNis the second dimension ofX.

#### — Input: REAL **X**(*M*N*)

On input:

Xcontains theMbyNreal 2D array to be transformed. Elementijis stored in locationi + (j – 1) * MofX.

#### — Output: COMPLEX **Y**((*M/2+1*)**N*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: REAL **COMM**(*4*M+10*N+300*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthM*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT2D(0,M,N,X,Y,COMM,INFO) CALL SCFFT2D(1,M,N,X,Y,COMM,INFO) DO J = 1, N DO I = 1, M/2+1 Y(I,J) = -Y(I,J)/SQRT(REAL(M*N)) END DO EMD DO CALL CSFFT2D(2,M,N,Y,X,COMM,INFO)

## SCFFT3D Routine Documentation

### SUBROUTINE:**SCFFT3D** (*MODE,L,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`SCFFT3D`

.

On input:

MODE=0 : only default initializations (specific toL,MandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real-to-complex transform is performed. Initializations are assumed to have been performed by a prior call to`SCFFT3D`

.MODE=2 : (default) initializations and a real-to-complex transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **L**

On input:

Lis the length of the first dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the first dimension ofX.

#### — Input: INTEGER **M**

On input:

Mis the length of the second dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the second dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the length of the third dimension of the 3D array of data to be transformed. IfXis declared as a 3D array thenLis the third dimension ofX.

#### — Input: REAL **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNreal 3D array to be transformed. Elementijkis stored in locationi + (j – 1)* L + (k – 1) * L * MofX.

#### — Output: COMPLEX **Y**((*L/2+1*)**M*N*)

On output:

Ycontains the transformed Hermitian sequences in complex-Hermitian storage.

#### — Input/Output: REAL **COMM**(*4*L+10*M+10*N+500*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthL*M*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFT3D(0,L,M,N,X,Y,COMM,INFO) CALL SCFFT3D(1,L,M,N,X,Y,COMM,INFO) DO K = 1, N DO J = 1, M DO I = 1, L/2+1 Y(I,J,K) = -Y(I,J,K)/SQRT(REAL(L*M*N)) END DO END DO EMD DO CALL CSFFT3D(2,L,M,N,Y,X,COMM,INFO)

## SCFFTM Routine Documentation

### SUBROUTINE:**SCFFTM** (*M,N,X,COMM,INFO*)

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the real sequences inX

#### — Input/Output: REAL **X**(*N*M*)

On input:

Xcontains theMreal sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.

On output:Xcontains the transformed Hermitian sequences.

#### — Input/Output: REAL **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL SCFFTM(1,N,X,COMM,INFO) CALL SCFFTM(2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*X(N-I+1,2) 10 CONTINUE CALL CSFFTM(1,N,X(1,3),COMM,INFO)

## UGEN

Specification for a user supplied base generator.

### SUBROUTINE:**UGEN** (*N,STATE,X,INFO*)

#### — Input: INTEGER **N**

On input: the number of random numbers to be generated.

#### — Input/Output: INTEGER **STATE**(***)

On input: the internal state of your generator.

#### — Output: DOUBLE PRECISION **X**(*N*)

On output: the array of

Nuniform distributed random numbers, each in the semi-open interval (0.0, 1.0] – i.e. 1.0 is a legitimate return value, but 0.0 is not.

#### — Output: INTEGER **INFO**

On output: a flag which you can use to signal an error in the call of

`UGEN`

– for example, if`UGEN`

is called without being initialized by`UINI`

.

## UINI

Specification for a user supplied initialization routine.

### SUBROUTINE:**UINI** (*GENID,SUBID,SEED,LSEED,STATE,LSTATE,INFO*)

#### — Input: INTEGER **GENID**

On input: the ID associated with the generator. It may be used for anything you like.

#### — Input: INTEGER **SUBID**

On input: the sub-ID associated with the generator. It may be used for anything you like.

#### — Input: INTEGER **SEED**(*LSEED*)

On input: an array containing the initial seed for your generator.

#### — Input/Output: INTEGER **LSEED**

On input: either the size of the

SEEDarray, or a value < 1.

On output: ifLSEED<1 on entry,LSEEDmust be set to the required size of theSEEDarray. This allows a caller of`UINI`

to query the required size.

#### — Output: INTEGER **STATE**(*LSTATE*)

On output: if

LSTATE<1 on entry,STATEshould be unchanged.

Otherwise,STATEis a state vector holding internal details required by your generator. On exit from`UINI`

, the arraySTATEmust hold the following information:

`STATE(1)`

=`ESTATE`

, where`ESTATE`

is your minimum allowed size of arraySTATE.

`STATE(2)`

=`MAGIC`

, where`MAGIC`

is a magic number of your own choice. This can be used by your routine`UGEN`

as a check that`UINI`

has previously been called.

`STATE(3)`

=`GENID`

`STATE(4)`

=`SUBID`

`STATE(5)`

…`STATE(ESTATE-1)`

= internal state values required by your generator routine`UGEN`

; for example, the current value of your seed.

`STATE(ESTATE)`

=`MAGIC`

, i.e. the same value as`STATE(2)`

.

#### — Input/Output: INTEGER **LSTATE**

On input: either the size of the

STATEarray, or a value < 1.

On output: ifLSTATE<1 on entry,LSTATEshould be set to the required size of theSTATEarray, i.e. the value`ESTATE`

as described above. This allows the caller of`UINI`

to query the required size.

Constraint: eitherLSTATE<1 orLSTATE>=ESTATE .

#### — Output: INTEGER **INFO**

On output: an error code, to be used in whatever way you wish; for example to flag an incorrect argument to

`UINI`

. If no error is encountered,`UINI`

must setINFOto 0.

## ZDFFT Routine Documentation

### SUBROUTINE:**ZDFFT** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZDFFT`

.

On input:

MODE=0 : only initializations (specific to the values ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a real transform is performed. Initializations are assumed to have been performed by a prior call to`ZDFFT`

.MODE=2 : (default) initializations and a real transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inX

#### — Input/Output: DOUBLE PRECISION **X**(*N*)

On input:

Xcontains the Hermitian sequence of lengthNto be transformed.

On output:Xcontains the transformed real sequence.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT(0,N,X,COMM,INFO) CALL DZFFT(1,N,X,COMM,INFO) DO 10 I = N/2+2, N X(I) = -X(I) 10 CONTINUE CALL ZDFFT(2,N,X,COMM,INFO)

## ZDFFT1D Routine Documentation

### SUBROUTINE:**ZDFFT1D** (*MODE,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZDFFT1D`

.

On input:

MODE=0 : only initializations (specific to the values ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`ZDFFT1D`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inY

#### — Input: DOUBLE COMPLEX **X**(*N/2+1*)

On input:

Xcontains the complex sequence to be transformed.

#### — Output: DOUBLE PRECISION **Y**(*N*)

On output:

Ycontains the transformed real sequence of lengthN.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT1D(0,N,X,Y,COMM,INFO) CALL DZFFT1D(1,N,X,Y,COMM,INFO) DO 10 I = 1, N/2+1 Y(I) = -Y(I)*EXP(-DBLE(I-1)/DBLE(N)) 10 CONTINUE CALL ZDFFT1D(2,N,Y,X,COMM,INFO)

## ZDFFT1M Routine Documentation

### SUBROUTINE:**ZDFFT1M** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZDFFT1M`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`ZDFFT1M`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis length of the sequence inY

#### — Input: DOUBLE COMPLEX **X**((*N/2+1*)**M*)

On input:

Xcontains the Hermitian sequences to be transformed in complex-Hermitian storage.

#### — Output: DOUBLE PRECISION **Y**(*N*M*)

On output:

Ycontains theMreal sequences of lengthNtransformed. Elementiof sequencejis stored in locationi + (j – 1) * NofY.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT1M(0,M,N,X,Y,COMM,INFO) CALL DZFFT1M(1,M,N,X,Y,COMM,INFO) DO J = 1, M DO I = 1, N/2+1 Y(I,J) = -Y(I,J)*EXP(-DBLE(I-1)/DBLE(N)) END DO EMD DO CALL ZDFFT1M(2,M,N,Y,X,COMM,INFO)

## ZDFFT2D Routine Documentation

### SUBROUTINE:**ZDFFT2D** (*MODE,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZDFFT2D`

.

On input:

MODE=0 : only default initializations (specific toMandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`ZDFFT2D`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **M**

On input:

Mis the number of rows in the 2D array of the real data obtained from the transform. IfYis declared as a 2D array thenMis the first dimension ofY.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of the real data obtained from the transform. IfYis declared as a 2D array thenNis the second dimension ofY.

#### — Input: DOUBLE COMPLEX **X**((*M/2+1*)**N*)

On input:

Xcontains the Hermitian sequences in complex-Hermitian storage to be transformed.

#### — Output: DOUBLE PRECISION **Y**(*M*N*)

On output:

Ycontains theMbyNreal 2D array obtained from the transform. Elementijis stored in locationi + (j – 1) * MofY.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*M+6*N+300*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthM*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT2D(0,M,N,X,Y,COMM,INFO) CALL DZFFT2D(1,M,N,X,Y,COMM,INFO) DO J = 1, N DO I = 1, M/2+1 Y(I,J) = -Y(I,J)/SQRT(DBLE(M*N)) END DO EMD DO CALL ZDFFT2D(2,M,N,Y,X,COMM,INFO)

## ZDFFT3D Routine Documentation

### SUBROUTINE:**ZDFFT3D** (*MODE,L,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZDFFT3D`

.

On input:

MODE=0 : only default initializations (specific toL,MandN) are performed; this is usually followed by calls to the same routine withMODE=1.MODE=1 : a complex-to-real transform is performed. Initializations are assumed to have been performed by a prior call to`ZDFFT3D`

.MODE=2 : (default) initializations and a complex-to-real transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **L**

On input:

Lis the length of the first dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the first dimension ofY.

#### — Input: INTEGER **M**

On input:

Mis the length of the second dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the second dimension ofY.

#### — Input: INTEGER **N**

On input:

Nis the length of the third dimension of the 3D array of data obtained from the transform. IfYis declared as a 3D array thenLis the third dimension ofY.

#### — Input: DOUBLE COMPLEX **X**((*L/2+1*)**M*N*)

On input:

Xcontains the Hermitian sequences in complex-Hermitian storage to be transformed.

#### — Output: DOUBLE PRECISION **Y**(*L*M*N*)

On output:

Ycontains theLbyMbyNreal 3D array obtained from the transform. Elementijkis stored in locationi + (j – 1) * L + (k – 1) * L * MofY.

#### — Input/Output: DOUBLE PRECISION **COMM**(*4*L+6*M+6*N+500*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthL*M*N. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFT3D(0,L,M,N,X,Y,COMM,INFO) CALL DZFFT3D(1,L,M,N,X,Y,COMM,INFO) DO K = 1, N DO J = 1, M DO I = 1, L/2+1 Y(I,J,K) = -Y(I,J,K)/SQRT(DBLE(L*M*N)) END DO END DO EMD DO CALL ZDFFT3D(2,L,M,N,Y,X,COMM,INFO)

## ZDFFTM Routine Documentation

### SUBROUTINE:**ZDFFTM** (*M,N,X,COMM,INFO*)

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the sequences inX

#### — Input/Output: DOUBLE PRECISION **X**(*N*M*)

On input:

Xcontains theMHermitian sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.

On output:Xcontains the transformed real sequences.

#### — Input/Output: DOUBLE PRECISION **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL DZFFTM(1,N,X,COMM,INFO) CALL DZFFTM(2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*X(N-I+1,2) 10 CONTINUE CALL ZDFFTM(1,N,X(1,3),COMM,INFO)

## ZFFT1D Routine Documentation

### SUBROUTINE:**ZFFT1D** (*MODE,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT1D`

.

On input:

MODE=0 : only default initializations (specific toN) are performed; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1D`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1D`

.MODE=-2 : initializations and a forward transform are performed.MODE=2 : initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations are performed, but first a plan is generated. This plan is chosen based on the fastest FFT computation for a subset of all possible plans.

#### — Input: INTEGER **N**

On input:

Nis the length of the complex sequenceX

#### — Input/Output: COMPLEX*16 **X**(*N*)

On input:

Xcontains the complex sequence of lengthNto be transformed.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX*16 **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL ZFFT1D(0,N,X,COMM,INFO) CALL ZFFT1D(-1,N,X,COMM,INFO) CALL ZFFT1D(-1,N,Y,COMM,INFO) DO 10 I = 1, N X(I) = X(I)*DCONJG(Y(I)) 10 CONTINUE CALL ZFFT1D(1,N,X,COMM,INFO)

## ZFFT1DX Routine Documentation

### SUBROUTINE:**ZFFT1DX** (*MODE,SCALE,INPL,N,X,INCX,Y,INCY,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT1DX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1DX`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1DX`

.MODE=-2 : (default) initializations and a forward transform are performed.MODE=2 : (default) initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: DOUBLE PRECISION **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequence

#### — Input: LOGICAL **INPL**

On input: if

INPLis .TRUE. thenXis overwritten by the output sequence; otherwise the output sequence is returned inY.

#### — Input: INTEGER **N**

On input:

Nis the number of elements to be transformed

#### — Input/Output: COMPLEX*16 **X**(*1+*(*N-1*)**INCX*)

On input:

Xcontains the complex sequence of lengthNto be transformed, with the ith element stored in X(1+(i-1)*INCX).

On output: ifINPLis .TRUE. thenXcontains the transformed sequence in the same locations as on input; otherwiseXremains unchanged.

#### — Input: INTEGER **INCX**

On input:

INCXis the increment used to store successive elements of a sequence inX.

Constraint:INCX> 0.

#### — Output: COMPLEX*16 **Y**(*1+*(*N-1*)**INCY*)

On output: if

INPLis .FALSE. thenYcontains the transformed sequence, with the ith element stored in Y(1+(i-1)*INCY); otherwiseYis not referenced.

#### — Input: INTEGER **INCY**

On input:

INCYis the increment used to store successive elements of a sequence inY. IfINPLis .TRUE. thenINCYis not referenced.

Constraint:INCY> 0.

#### — Input/Output: COMPLEX*16 **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

**INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward FFTs are performed unscaled and in-place on contiguous C vectors X and Y following initialization. Manipulations on C resultant Fourier coefficients are stored in X which is then C transformed back. C SCALE = 1.0D0 INPL = .TRUE. CALL ZFFT1DX(0,SCALE,INPL,N,X,1,DUM,1,COMM,INFO) CALL ZFFT1DX(-1,SCALE,INPL,N,X,1,DUM,1,COMM,INFO) CALL ZFFT1DX(-1,SCALE,INPL,N,Y,1,DUM,1,COMM,INFO) DO 10 I = 1, N X(I) = X(I)*DCONJG(Y(I))/DBLE(N) 10 CONTINUE CALL ZFFT1DX(1,SCALE,INPL,N,X,1,DUM,1,COMM,INFO)

## ZFFT1M Routine Documentation

### SUBROUTINE:**ZFFT1M** (*MODE,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT1M`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : forward transforms are performed. Initializations are assumed to have been performed by a prior call to`ZFFT1M`

.MODE=1 : backward (reverse) transforms are performed. Initializations are assumed to have been performed by a prior call to`ZFFT1M`

.MODE=-2 : (default) initializations and forward transforms are performed.MODE=2 : (default) initializations and backward transforms are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: INTEGER **M**

On input:

Mis the number of sequences to be transformed.

#### — Input: INTEGER **N**

On input:

Nis the length of the complex sequences inX

#### — Input/Output: COMPLEX*16 **X**(*N*M*)

On input:

Xcontains theMcomplex sequences of lengthNto be transformed. Element i of sequence j is stored in location i+(j-1)*NofX.

On output:Xcontains the transformed sequences.

#### — Input/Output: COMPLEX*16 **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL ZFFT1M(0,1,N,X,COMM,INFO) CALL ZFFT1M(-1,2,N,X,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*DCONJG(X(I,2)) X(I,2) = DCMPLX(0.0D0,1.0D0)*X(I,2) 10 CONTINUE CALL ZFFT1M(1,2,N,X(1,2),COMM,INFO)

## ZFFT1MX Routine Documentation

### SUBROUTINE:**ZFFT1MX** (*MODE,SCALE,INPL,NSEQ,N,X,INCX1,INCX2,Y,INCY1,INCY2,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT1MX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1MX`

.MODE=1 : a backward (reverse) transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT1MX`

.MODE=-2 : (default) initializations and a forward transform are performed.MODE=2 : (default) initializations and a backward transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the value ofN) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the value ofN.

#### — Input: DOUBLE PRECISION **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **INPL**

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **NSEQ**

On input:

NSEQis the number of sequences to be transformed

#### — Input: INTEGER **N**

On input:

Nis the number of elements in each sequence to be transformed

#### — Input/Output: COMPLEX*16 **X**(*1+*(*N-1*)**INCX1+*(*NSEQ-1*)**INCX2*)

On input:

Xcontains theNSEQcomplex sequences of lengthNto be transformed; the ith element of sequence j is stored in X(1+(i-1)*INCX1+(j-1)*INCX2).

On output: ifINPLis .TRUE. thenXcontains the transformed sequences in the same locations as on input; otherwiseXremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the increment used to store successive elements of a given sequence inX(INCX1=1 for contiguous data).

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the increment used to store corresponding elements of successive sequences inX(INCX2=Nfor contiguous data).

Constraint:INCX2> 0.

#### — Output: COMPLEX*16 **Y**(*1+*(*N-1*)**INCY1+*(*NSEQ-1*)**INCY2*)

On output: if

INPLis .FALSE. thenYcontains the transformed sequences with the ith element of sequence j stored in Y(1+(i-1)*INCY1+(j-1)*INCY2); otherwiseYis not referenced.

#### — Input: INTEGER **INCY1**

On input:

INCY1is the increment used to store successive elements of a given sequence inY. IfINPLis .TRUE. thenINCY1is not referenced.

Constraint:INCY1> 0.

#### — Input: INTEGER **INCY2**

On input:

INCY2is the increment used to store corresponding elements of successive sequences inY(INCY2=Nfor contiguous data). IfINPLis .TRUE. thenINCY2is not referenced.

Constraint:INCY2> 0.

#### — Input/Output: COMPLEX*16 **COMM**(*3*N+100*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence lengthN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward FFTs are performed unscaled and in-place on two C contiguous vectors stored in the first two columns of X. C Manipulations are stored in 2nd and 3rd columns of X which are C then transformed back. C COMPLEX *16 X(N,3) SCALE = 1.0D0 INPL = .TRUE. CALL ZFFT1MX(0,SCALE,INPL,2,N,X,1,N,DUM,1,N,COMM,INFO) CALL ZFFT1MX(-1,SCALE,INPL,2,N,X,1,N,DUM,1,N,COMM,INFO) DO 10 I = 1, N X(I,3) = X(I,1)*DCONJG(X(I,2))/DBLE(N) X(I,2) = DCMPLX(0.0D0,1.0D0)*X(I,2)/DBLE(N) 10 CONTINUE CALL ZFFT1MX(1,SCALE,INPL,2,N,X(1,2),1,N,DUM,1,N,COMM,INFO)

## ZFFT2D Routine Documentation

### SUBROUTINE:**ZFFT2D** (*MODE,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the direction of transform to be performed by`ZFFT2D`

.

On input:

MODE=-1 : forward 2D transform is performed.MODE=1 : backward (reverse) 2D transform is performed.

#### — Input: INTEGER **M**

Mis the number of rows in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the first dimension ofX.

#### — Input: INTEGER **N**

On input:

Nis the number of columns in the 2D array of data to be transformed. IfXis declared as a 2D array thenMis the second dimension ofX.

#### — Input/Output: COMPLEX*16 **X**(*M*N*)

On input:

Xcontains theMbyNcomplex 2D array to be transformed. Element ij is stored in location i+(j-1)*MofX.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX*16 **COMM**(*M*N+3**(*M+N*)*+100*)

COMMis a communication array used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL ZFFT2D(-1,M,N,X,COMM,INFO) DO 20 J = 1, N DO 10 I = 1, MIN(J-1,M) X(I,J) = 0.5D0*(X(I,J) + X(J,I)) X(J,I) = DCONJG(X(I,J)) 10 CONTINUE 20 CONTINUE CALL ZFFT2D(1,M,N,X,COMM,INFO)

## ZFFT2DX Routine Documentation

### SUBROUTINE:**ZFFT2DX** (*MODE,SCALE,LTRANS,INPL,M,N,X,INCX1,INCX2,Y,INCY1,INCY2,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT2DX`

.

On input:

MODE=0 : only initializations (specific to the value ofN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 2D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT2DX`

.MODE=1 : a backward (reverse) 2D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT2DX`

.MODE=-2 : (default) initializations and a forward 2D transform are performed.MODE=2 : (default) initializations and a backward 2D transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the values ofNandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofNandM.

#### — Input: DOUBLE PRECISION **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **LTRANS**

On input: if

LTRANSis .TRUE. then a normal final transposition is performed internally to return transformed data consistent with the values for argumentsINPL,INCX1,INCX2,INCY1andINCY2. IfLTRANSis .FALSE. then the final transposition is not performed explicitly; the storage format on output is determined by whether the output data is stored contiguously or not – please see the output specifications forXandYfor details.

#### — Input: LOGICAL **INPL**

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **M**

On input:

Mis the first dimension of the 2D transform.

#### — Input: INTEGER **N**

On input:

Nis the second dimension of the 2D transform.

#### — Input/Output: COMPLEX*16 **X**(*1+*(*M-1*)**INCX1+*(*N-1*)**INCX2*)

On input:

Xcontains theMbyNcomplex 2D data array to be transformed; the (ij)th element is stored in X(1+(i-1)*INCX1+(j-1)*INCX2).

On output: ifINPLis .TRUE. thenXcontains the transformed data, either in the same locations as on input whenLTRANS=.TRUE.; in locations X((i-1)*N+j) whenLTRANS=.FALSE.,INCX1=1 andINCX2=M; and otherwise in the same locations as on input. IfINPLis .FALSE.Xremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the increment used to store, inX, successive elements in the first dimension (INCX1=1 for contiguous data).

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the increment used to store, inX, successive elements in the second dimension (INCX2=M for contiguous data).

Constraint:INCX2> 0;

INCX2> (M-1)*INCX1 ifN> 1.

#### — Output: COMPLEX*16 **Y**(*1+*(*M-1*)**INCY1+*(*N-1*)**INCY2*)

On output: if

INPLis .FALSE. thenYcontains the transformed data. IfLTRANS=.TRUE. then the (ij)th data element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2); ifLTRANS=.FALSE.,INCY1=1 andINCY2=N then the (ij)th data element is stored in Y((i-1)*N+j); and otherwise the (ij)th element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2). IfINPLis .TRUE. thenYis not referenced.

#### — Input: INTEGER **INCY1**

On input:

INCY1is the increment used to store successive elements in the first dimension inY(INCY1=1 for contiguous data). IfINPLis .TRUE. thenINCY1is not referenced.

Constraint:INCY1> 0.

#### — Input: INTEGER **INCY2**

On input:

INCY2is the increment used to store successive elements in the second dimension inY(for contiguous data,INCY2=MwhenLTRANSis .TRUE. orINCY2=NwhenLTRANSis .FALSE.). IfINPLis .TRUE. thenINCY2is not referenced.

Constraints:INCY2> 0;

INCY2> (M-1)*INCY1 ifN> 1 andLTRANSis .TRUE.;

INCY2= N ifM> 1 andLTRANSis .FALSE..

#### — Input/Output: COMPLEX*16 **COMM**(*M*N+3*M+3*N+200*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same dimensionsMandN. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 2D FFT is performed unscaled, without final transpose C and out-of-place on data stored in array X and output to Y. C Manipulations are stored in vector Y which is then transformed C back, with scaling, into the first M rows of X. C COMPLEX *16 X(M,N), Y(N,M) SCALE = 1.0D0 INPL = .FALSE. LTRANS = .FALSE. CALL ZFFT2DX(0,SCALE,LTRANS,INPL,M,N,X,1,M,Y,1,N,COMM,INFO) CALL ZFFT2DX(-1,SCALE,LTRANS,INPL,M,N,X,1,M,Y,1,N,COMM,INFO) DO 20 I = M DO 10 J = 1, N Y(J,I) = 0.5D0*Y(J,I)*EXP(0.001D0*(I+J-2)) 10 CONTINUE 20 CONTINUE SCALE = 1.0D0/DBLE(M*N) CALL ZFFT2DX(1,SCALE,LTRANS,INPL,N,M,Y,1,N,X,1,M,COMM,INFO)

## ZFFT3D Routine Documentation

### SUBROUTINE:**ZFFT3D** (*MODE,L,M,N,X,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the direction of transform to be performed by`ZFFT3D`

.

On input:

MODE=-1 : forward 3D transform is performed.MODE=1 : backward (reverse) 3D transform is performed.

#### — Input: INTEGER **L**

On input: the length of the first dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenLis the first dimension ofX.

#### — Input: INTEGER **M**

On input: the length of the second dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenMis the second dimension ofX.

#### — Input: INTEGER **N**

On input: the length of the third dimension of the 3D array of data to be transformed. If

Xis declared as a 3D array thenNis the third dimension ofX.

#### — Input/Output: COMPLEX*16 **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNcomplex 3D array to be transformed. Element ijk is stored in location i+(j-1)*L+(k-1)*L*MofX.

On output:Xcontains the transformed sequence.

#### — Input/Output: COMPLEX*16 **COMM**(*L*M*N+3**(*L+M+N*))

COMMis a communication array used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

CALL ZFFT3D(-1,L,M,N,X,COMM,INFO) DO 30 K = 1, N DO 20 J = 1, M DO 10 I = 1, L X(I,J) = X(I,J)*EXP(-0.001D0*DBLE(I+J+K)) 10 CONTINUE 20 CONTINUE 30 CONTINUE CALL ZFFT3D(1,L,M,N,X,COMM,INFO)

## ZFFT3DX Routine Documentation

### SUBROUTINE:**ZFFT3DX** (*MODE,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT3DX`

.

On input:

MODE=0 : only initializations (specific to the values ofL,MandN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 3D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT3DX`

.MODE=1 : a backward (reverse) 3D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT3DX`

.MODE=100 : similar toMODE=0; only initializations (specific to the values ofL,MandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofL,MandN.

#### — Input: DOUBLE PRECISION **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **LTRANS**

On input: if

LTRANSis .TRUE. then a normal final transposition is performed internally to return transformed data using the same storage format as the input data. IfLTRANSis .FALSE. then the final transposition is not performed and transformed data is stored, inXorY, in transposed form.

#### — Input: LOGICAL **INPL**

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **L**

On input:

Lis the first dimension of the 3D transform.

#### — Input: INTEGER **M**

On input:

Mis the second dimension of the 3D transform.

#### — Input: INTEGER **N**

On input:

Nis the third dimension of the 3D transform.

#### — Input/Output: COMPLEX*16 **X**(*L*M*N*)

On input:

Xcontains theLbyMbyNcomplex 3D data array to be transformed; the (ijk)th element is stored in X(i+(j-1)*L+(k-1)*L*M).

On output: ifINPLis .TRUE. thenXcontains the transformed data, either in the same locations as on input whenLTRANS=.TRUE.; or in locations X(k+(j-1)*N+(i-1)*N*M) whenLTRANS=.FALSE. IfINPLis .FALSE.Xremains unchanged.

#### — Output: COMPLEX*16 **Y**(*L*M*N*)

On output: if

INPLis .FALSE. thenYcontains the three-dimensional transformed data. IfLTRANS=.TRUE. then the (ijk)th data element is stored in Y(i+(j-1)*L+(k-1)*L*M); otherwise, the (ijk)th data element is stored in Y(k+(j-1)*N+(i-1)*N*M). IfINPLis .TRUE. thenYis not referenced.

#### — Input/Output: COMPLEX*16 **COMM**(*L*M*N+3**(*L+M+N*)*+300*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence dimensions. The remainder is used as temporary store.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 3D FFT is performed unscaled, without final transpose C and out-of-place on data stored in array X and output to Y. C Manipulations are stored in vector Y which is then transformed C back, with scaling, into the first M rows of X. C COMPLEX *16 X(L*M*N), Y(L*M*N) SCALE = 1.0D0 INPL = .FALSE. LTRANS = .FALSE. CALL ZFFT3DX(0,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO) CALL ZFFT3DX(-1,SCALE,LTRANS,INPL,L,M,N,X,Y,COMM,INFO) IY = 1 DO 20 I = 1, L DO 40 J = 1, M DO 10 K = 1, N Y(IY) = Y(IY)*EXP(-0.001D0*DBLE(I+J+K-3)) IY = IY + 1 10 CONTINUE 20 CONTINUE SCALE = 1.0D0/DBLE(L*M*N) CALL ZFFT3DX(1,SCALE,LTRANS,INPL,N,M,L,Y,X,COMM,INFO)

## ZFFT3DY Routine Documentation

### SUBROUTINE:**ZFFT3DY** (*MODE,SCALE,INPL,L,M,N,X,INCX1,INCX2,INCX3,Y,INCY1,INCY2,INCY3,COMM,LCOMM,INFO*)

#### — Input: INTEGER **MODE**

The value of

MODEon input determines the operation performed by`ZFFT3DY`

.

On input:

MODE=0 : only initializations (specific to the values ofL,MandN) are performed using a default plan; this is usually followed by calls to the same routine withMODE=-1 or 1.MODE=-1 : a forward 3D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT3DY`

.MODE=1 : a backward (reverse) 3D transform is performed. Initializations are assumed to have been performed by a prior call to`ZFFT3DY`

.MODE=-2 : (default) initializations and a forward 3D transform are performed.MODE=2 : (default) initializations and a backward 3D transform are performed.MODE=100 : similar toMODE=0; only initializations (specific to the values ofL,MandM) are performed, but these are based on a plan that is first generated by timing a subset of all possible plans and choosing the quickest (i.e. the FFT computation was timed as fastest based on the chosen plan). The plan generation phase may take a significant amount of time depending on the values ofL,MandN.

#### — Input: REAL **SCALE**

On input:

SCALEis the scaling factor to apply to the output sequences

#### — Input: LOGICAL **INPL**

INPLis .TRUE. thenXis overwritten by the output sequences; otherwise the output sequences are returned inY.

#### — Input: INTEGER **L**

On input:

Lis the first dimension of the 3D transform.

#### — Input: INTEGER **M**

On input:

Mis the second dimension of the 3D transform.

#### — Input: INTEGER **N**

On input:

Nis the third dimension of the 3D transform.

#### — Input/Output: COMPLEX*16 **X**(***)

On input:

Xcontains theLbyMbyNcomplex 3D data array to be transformed; the (ijk)th element is stored in X(1+(i-1)*INCX1+(j-1)*INCX2+(k-1)*INCX3).

On output: ifINPLis .TRUE. thenXcontains the transformed data in the same locations as on input. IfINPLis .FALSE.Xremains unchanged.

#### — Input: INTEGER **INCX1**

On input:

INCX1is the step in index ofXbetween successive data elements in the first dimension of the 3D data. UsuallyINCX1=1 so that succesive elements in the first dimension are stored contiguously.

Constraint:INCX1> 0.

#### — Input: INTEGER **INCX2**

On input:

INCX2is the step in index ofXbetween successive data elements in the second dimension of the 3D data. For completely contiguous data (no gaps inX)INCX2should be set toL.

Constraint:INCX2> 0;

INCX2> (L-1)*INCX1 if max(M,N) > 1.

#### — Input: INTEGER **INCX3**

On input:

INCX3is the step in index ofXbetween successive data elements in the third dimension of the 3D data. For completely contiguous data (no gaps inX)INCX3should be set toL*M.

Constraint:INCX3> 0;

INCX3> (L-1)*INCX1+(M-1)*INCX2if N > 1.

#### — Output: COMPLEX*16 **Y**(***)

On output: if

INPLis .FALSE. thenYcontains the three-dimensional transformed data. IfLTRANS=.TRUE. then the the (ijk)th element is stored in Y(1+(i-1)*INCY1+(j-1)*INCY2+(k-1)*INCY3).

IfINPLis .TRUE. thenYis not referenced.

#### — Input: INTEGER **INCY1**

On input: if

INPLis .FALSE. thenINCY1is the step in index ofYbetween successive data elements in the first dimension of the 3D transformed data. UsuallyINCY1=1 so that succesive elements in the first dimension are stored contiguously.

IfINPLis .TRUE. thenINCY1is not referenced. Constraint: IfINPLis .FALSE. thenINCY1> 0.

#### — Input: INTEGER **INCY2**

On input: if

INPLis .FALSE. thenINCY2is the step in index ofYbetween successive data elements in the second dimension of the 3D transformed data. For completely contiguous data (no gaps inY)INCY2should be set toL.

Constraint:INCY2> 0 ifINPLis .FALSE.;

INCY2> (L-1)*INCY1, ifINPLis .FALSE. and max(M,N) > 1.

#### — Input: INTEGER **INCY3**

On input: if

INPLis .FALSE. thenINCY3is the step in index ofYbetween successive data elements in the third dimension of the 3D transformed data. For completely contiguous data (no gaps inY)INCY3should be set toL*M.

Constraint:INCY3> 0 ifINPLis .FALSE.;

INCY3> (L-1)*INCY1+(M-1)*INCY2, ifINPLis .FALSE. and N > 1.

#### — Input/Output: COMPLEX*16 **COMM**(*LCOMM*)

COMMis a communication array. Some portions of the array are used to store initializations for subsequent calls with the same sequence dimensions. The remainder is used as temporary store; if this is not sufficient for the requirements of the routine then temporary storage space will be dynamically allocated internally.

#### — Input: INTEGER **LCOMM**

On input:

LCOMMis the length of the communication arrayCOMM. The amount of internal dynamic allocation of temporary storage can be reduced significantly by declaringCOMMto be of length at least L*M*N + 4*(L+M+N) + 300.

Constraint:LCOMM> 3*(L+M+N) + 150.

#### — Output: INTEGER **INFO**

INFOis an error indicator. On successful exit,INFOcontains 0. IfINFO= -i on exit, the i-th argument had an illegal value.

**Example:**

C Forward 3D FFT is performed unscaled and in-place, on the leading C 10x10x10 submatrix of a larger 100x100x100 array of data. C The result is transformed back with scaling. C SCALE = 1.0D0 INPL = .TRUE. L = 10 M = 10 N = 10 LCOMM = 2000000 CALL ZFFT3DY(0,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO) CALL ZFFT3DY(-1,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO) IY = 1 DO 20 I = 1, L DO 40 J = 1, M DO 10 K = 1, N X(I,J,K) = X(I,J,K)*EXP(-1.0D-3*DBLE(I+J+K-3)) 10 CONTINUE 20 CONTINUE SCALE = 1.0/DBLE(L*M*N) CALL ZFFT3DY(1,SCALE,INPL,L,M,N,X,1,100,10000,Y,1,1,1, * COMM,LCOMM,INFO)