Parallel Engineering and Scientific Subroutine Library for AIX Version 2 Release 3: Guide and Reference

PSCFT3 and PDCFT3--Complex Fourier Transforms in Three Dimensions

These subroutines compute the mixed-radix three-dimensional discrete Fourier transform of complex data:

for:

k1 = 0, 1, ..., n1-1

k2 = 0, 1, ..., n2-1

k3 = 0, 1, ..., n3-1

where:

and where:

x_j1,j2,j3 are elements of array X.

y_k1,k2,k3 are elements of array Y.

Isign is + or - (determined by argument isign).

scale is a scalar value.

For scale = 1 and isign being positive, you obtain the discrete Fourier transform. For scale = 1/((n1)(n2)(n3)) and isign being negative, you obtain the inverse Fourier transform.

See references [1] and [3].

Table 112. Data Types

Syntax

On Entry

x

is the local array X, containing the three-dimensional data to be transformed that has been block-plane distributed over a 1 × q process grid, where q is the number of processes. (The values of ldx1 and ldx2 are set in the IP array.)

Scope: local

Specified as: an array of (at least) length ldx1 × ldx2 × LOCq(n3), containing numbers of the data type indicated in Table 112. This array must be aligned on a doubleword boundary.

y

See On Return.

n1

is the length of the first dimension of the three-dimensional data of the array to be transformed.

Scope: global

Specified as: a fullword integer; n1 <= 37748736 and must be one of the values listed in Figure 9.

n2

is the length of the second dimension of the three-dimensional data of the array to be transformed.

Scope: global

Specified as: a fullword integer; n2 <= 37748736 and must be one of the values listed in Figure 9.

n3

is the length of the third dimension of the three-dimensional data of the array to be transformed.

Scope: global

Specified as: a fullword integer; n3 <= 37748736 and must be one of the values listed in Figure 9.

isign

controls the direction of the transform, determining the sign, isign, of the exponent of W_n, where:

If isign = positive value, Isign = + (transforming time to frequency).

If isign = negative value, Isign = - (transforming frequency to time).

Scope: global

Specified as: a fullword integer; where isign > 0 or isign < 0.

scale

is the scaling constant scale.

Scope: global

Specified as: a number of the data type indicated in Table 112, where scale > 0.0 or scale < 0.0.

icontxt

is the BLACS context parameter.

Scope: global

Specified as: the fullword integer that was returned by a prior call to BLACS_GRIDINIT or BLACS_GRIDMAP.

ip

is an array of parameters, IP(i), where:

• IP(1) indicates whether the default values for ip are used or you set the values for ip.

  If IP(1) = 0, then the following default values are used:

  • y is returned in transposed form; that is global y has dimensions n3 × n2 × n1.
  • ldx1 and ldx2, the leading dimensions of the array specified for X, equal n1 and n2, respectively.
  • ldy1 and ldy2, the leading dimensions of the array specified for Y, equal n3 and n2, respectively.

    The remaining parameters of array IP are ignored.

  If IP(1) <> 0, then you set the remaining values of ip to indicate whether y is stored in normal or transposed form, and indicate values for the leading dimensions.

• IP(2) indicates whether y is to be stored in normal or transposed form.

  If IP(2)=0, then y is to be stored in transposed form on output.

  If IP(2)=1, then y is to be stored in normal form on output.

• IP(3-19) are reserved.
• IP(20) indicates the values of the leading dimension, ldx1, for the array specified for X, where:

  If IP(20) = 0, then ldx1 = n1.

  If IP(20) <> 0, then ldx1 is this value of IP(20).

• IP(21) indicates the values of the leading dimension, ldx2, for the array specified for X, where:

  If IP(21) = 0, then ldx2 = n2.

  If IP(21) <> 0, then ldx2 is this value of IP(21).

• IP(22) indicates the values of the leading dimension, ldy1, for the array specified for Y, where:

  If IP(22) = 0 and IP(2) = 1, then ldy1 = n1.

  If IP(22) = 0 and IP(2) = 0, then ldy1 = n3.

  If IP(22) <> 0, then ldy1 is this value of IP(22).

• IP(23) indicates the values of the leading dimension, ldy2, for the array specified for Y, where:

  If IP(23) = 0, then ldy2 = n2.

  If IP(23) <> 0, then ldy2 is this value of IP(23).

• IP(24-40) are reserved.

Scope: global

Specified as: a one-dimensional array of (at least) length 40, containing fullword integers, where:

IP(1) is any integer

IP(2) = 0 or 1

IP(20) >= n1 or IP(20)=0

IP(21) >= n2 or IP(21)=0

IP(22) >= n1 (for normal form) or IP(22)=0

IP(22) >= n3 (for transposed form) or IP(22) = 0

IP(23) >= n2 or IP(23)=0

On Return

y

is the local array Y that is block-plane distributed and contains the results of the computation, where:

If IP(1) = 0, the local array Y is stored in transposed form and has dimensions n3 × n2 × LOCq(n1).

If IP(1) <> 0 and IP(2)=0, then the local array Y is stored in transposed form and has dimensions ldy1 × ldy2 × LOCq(n1).

If IP(1) <> 0 and IP(2)=1, then the local array Y is stored in normal form and has dimensions ldy1 × ldy2 × LOCq(n3).

Scope: local

Returned as: an ldy1 × ldy2 × LOCq(n3) array (for normal form) or an ldy1 × ldy2 × LOCq(n1) array (for transposed form), containing the numbers of the data type indicated in Table 112. This array must be aligned on a doubleword boundary.

Notes and Coding Rules

1. For the output array Y, these subroutines may use any extra space available when ldy1 and ldy2 are greater than their minimum value.
2. You may specify the same array for X and Y. In this case, output overwrites input. If you specify different arrays X and Y, they must have no common elements; otherwise, results are unpredictable.
3. For more information on LOCq(_) and how sequences are block-plane distributed, see Three-Dimensional Sequences.

   In general, distributing your data evenly provides the best work load balance among the processes and allows the use of the most efficient collective communication. However, for your specific problem size and number of processes available, experimentation is necessary to achieve optimal performance.

Error Conditions

Computational Errors

None

Resource Errors

1. Unable to allocate work space.

Input-Argument and Miscellaneous Errors

Stage 1:  

1. icontxt is invalid

Stage 2:  

1. Process grid is not 1 × q
2. The subroutine was called from outside the process grid.

Stage 3:  

1. n1 > 37748736
2. n2 > 37748736
3. n3 > 37748736
4. The length of n1, n2, or n3 is not an allowable transform length.
5. isign = 0
6. scale = 0.0
7. IP(1) <> 0 and IP(2) <> 0 or 1
8. IP(1) <> 0 and IP(20) <> 0 and IP(20) < n1 (that is, ldx1 < n1)
9. IP(1) <> 0 and IP(21) <> 0 and IP(21) < n2 (that is, ldx2 < n2)
10. IP(1) <> 0 and IP(2)=0 (for transpose mode) and IP(22) <> 0 and IP(22) < n3 (that is, ldy1 < n3)
11. IP(1) <> 0 and IP(2)=1 (for normal mode) and IP(22) <> 0 and IP(22) < n1 (that is, ldy1 < n1)
12. IP(1) <> 0 and IP(23) <> 0 and IP(23) < n2 (that is, ldy2 < n2)

Example 1

This example shows how to compute a three-dimensional transform. The data is block-plane distributed over a 1 × 2 process grid. The arrays are declared as follows:

  COMPLEX*16 X(0:3,0:3,0)
  COMPLEX*16 Y(0:3,0:3,0)
  INTEGER*4 IP(40)
  REAL*8 SCALE

Call Statements and Input

ORDER = 'R'
NPROW = 1
NPCOL = 2
CALL BLACS_GET(0, 0, ICONTXT)
CALL BLACS_GRIDINIT(ICONTXT, ORDER, NPROW, NPCOL)
CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYROW, MYCOL)
IP(1) = 1
IP(2) = 1
IP(20) = 4
IP(21) = 4
IP(22) = 4
IP(23) = 4
 
             X   Y   N1  N2  N3  ISIGN  SCALE   ICONTXT   IP
             |   |   |   |   |     |      |        |      |
CALL PDCFT3( X , Y , 4 , 4 , 2  ,  1  , 1.0D0 , ICONTXT , IP)

Global matrix X:

Plane 0:

B,D                          0
     *                                                *
     |  (1.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
 0   |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     *                                                *

Plane 1:

B,D                          1
     *                                                *
     |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
 0   |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     |  (0.0,0.0)   (0.0,0.0)   (0.0,0.0)   (0.0,0.0) |
     *                                                *

The following is the 1 × 2 process grid:

B,D  |    0    |  1  
-----| ------- |-----
0    |   P00   |  P01

Local arrays for X:

p,q  |                      0                       |                       1
-----|----------------------------------------------|----------------------------------------------
     |  (1.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  |   (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)
     |  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  |   (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)
 0   |  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  |   (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)
     |  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  |   (0.0,0.0)  (0.0,0.0)  (0.0,0.0)  (0.0,0.0)

Output

Global matrix Y:

Plane 0:

B,D                          0
     *                                                *
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
 0   |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     *                                                *

Plane 1:

B,D                          1
     *                                                *
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
 0   |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     |  (1.0,0.0)   (1.0,0.0)   (1.0,0.0)   (1.0,0.0) |
     *                                                *

The following is the 1 × 2 process grid:

B,D  |    0    |  1  
-----| ------- |-----
0    |   P00   |  P01

Local arrays for Y:

p,q  |                      0                       |                       1
-----|----------------------------------------------|----------------------------------------------
     |  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  |   (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)
     |  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  |   (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)
 0   |  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  |   (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)
     |  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  |   (1.0,0.0)  (1.0,0.0)  (1.0,0.0)  (1.0,0.0)

Example 2

This example shows how to compute a three-dimensional transform. In this example, the IP array is set to 0, which means array Y is returned in transposed form, ldx=n1, and ldy=n3. This is an example of uneven block-plane distribution over a 1 × 3 process grid. The arrays are declared as follows:

  COMPLEX*16 X(0:3,0:1,0:1), Y(0:5,0:1,0:1)
  INTEGER*4  IP(40)
  REAL*8     SCALE

Call Statements and Input

ORDER = 'R'
NPROW = 1
NPCOL = 3
CALL BLACS_GET(0, 0, ICONTXT)
CALL BLACS_GRIDINIT(ICONTXT, ORDER, NPROW, NPCOL)
CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYROW, MYCOL)
IP(1) = 0
 
             X   Y   N1  N2  N3  ISIGN     SCALE     ICONTXT   IP
             |   |   |   |   |     |         |          |      |
CALL PDCFT3( X , Y , 4 , 2 , 6  ,  1  , 1.0D0/8.0D0, ICONTXT , IP)

Global matrix X:

          Plane 0:                 Plane 1:
----------------------------------------------------
B,D                          0
----------------------------------------------------
     *                                             *
     | (4.9,3.9)  (5.9,3.1) | (2.2,5.1)  (5.9,1.5) |
0    | (6.8,4.6)  (6.7,9.8) | (2.1,0.5)  (3.5,2.9) |
     | (4.9,1.9)  (1.6,4.9) | (7.0,6.8)  (6.4,1.1) |
     | (2.9,7.6)  (7.5,5.5) | (0.6,4.6)  (9.0,7.6) |
     *                                             *
 
          Plane 2:                 Plane 3:
----------------------------------------------------
B,D                         1
----------------------------------------------------
     *                                             *
     | (8.3,2.3)  (7.3,7.3) | (4.6,5.9)  (3.3,3.0) |
0    | (4.5,8.9)  (0.2,0.9) | (3.6,5.0)  (3.4,0.4) |
     | (1.8,6.5)  (2.5,1.7) | (8.5,9.3)  (0.2,1.7) |
     | (6.4,5.2)  (7.8,5.3) | (8.7,9.1)  (9.1,5.5) |
     *                                             *
 
          Plane 4:                 Plane 5:
----------------------------------------------------
B,D                         2
----------------------------------------------------
     *                                             *
     | (1.0,1.0)  (2.2,8.2) | (5.4,3.6)  (4.8,4.0) |
0    | (8.8,1.0)  (8.9,5.8) | (2.0,2.6)  (5.2,9.5) |
     | (3.3,0.0)  (8.9,7.5) | (3.1,9.7)  (5.4,7.5) |
     | (8.7,8.3)  (6.3,5.4) | (1.5,9.6)  (4.4,4.4) |
     *                                             *

The following is the 1 × 3 process grid:

B,D  |    0    |    1    |    2    
-----| ------- | ------- |------- 
0    |   P00   |   P01   |   P02

Local arrays for X:

p,q|                   0                    |                   1                    |                   2
---|----------------------------------------|----------------------------------------|---------------------------------------
   |(4.9,3.9) (5.9,3.1) (2.2,5.1) (5.9,1.5) |(8.3,2.3) (7.3,7.3) (4.6,5.9) (3.3,3.0) |(1.0,1.0) (2.2,8.2) (5.4,3.6) (4.8,4.0)
   |(6.8,4.6) (6.7,9.8) (2.1,0.5) (3.5,2.9) |(4.5,8.9) (0.2,0.9) (3.6,5.0) (3.4,0.4) |(8.8,1.0) (8.9,5.8) (2.0,2.6) (5.2,9.5)
 0 |(4.9,1.9) (1.6,4.9) (7.0,6.8) (6.4,1.1) |(1.8,6.5) (2.5,1.7) (8.5,9.3) (0.2,1.7) |(3.3,0.0) (8.9,7.5) (3.1,9.7) (5.4,7.5)
   |(2.9,7.6) (7.5,5.5) (0.6,4.6) (9.0,7.6) |(6.4,5.2) (7.8,5.3) (8.7,9.1) (9.1,5.5) |(8.7,8.3) (6.3,5.4) (1.5,9.6) (4.4,4.4)

Output

Global matrix Y:

          Plane 0:                    Plane 1:
--------------------------------------------------------
B,D                           0
--------------------------------------------------------
     *                                                   *
     | (29.8,29.7)  (-1.9,1.1) |  (-3.0,0.9) (-3.0,-3.8) |
     |  (-3.3,1.0)  (0.4,-2.5) |  (4.1,-3.9) (-4.6,-1.4) |
0    | (-1.7,-1.2)   (0.1,3.4) |   (0.9,5.0)  (-0.6,1.6) |
     |  (2.3,-0.5)  (1.0,-4.6) |   (3.1,0.8)   (1.7,0.4) |
     |   (3.0,1.9)   (4.5,0.6) |  (0.5,-3.8) (-3.0,-0.1) |
     |   (1.0,0.1) (-5.7,-1.9) | (-1.4,-1.2)   (0.8,2.6) |
     *                                                   *
 
          Plane 2:                      Plane 3:
--------------------------------------------------------
B,D                            1
--------------------------------------------------------
     *                                                   *
     | (-2.4,-2.8)   (2.0,0.1) |  (3.6,-3.4)   (1.4,0.0) |
     |  (2.1,-3.5)   (2.8,1.3) | (-1.9,-0.1)   (2.3,6.3) |
0    |  (-1.7,0.7)   (2.8,1.0) |   (2.0,1.4)  (-0.6,1.0) |
     | (-3.3,-2.2) (-3.2,-5.1) |  (-0.3,3.3)   (0.8,0.5) |
     | (-1.5,-3.1)   (1.4,0.1) | (-1.3,-1.7)  (-2.7,0.7) |
     |   (1.8,0.6)  (-0.7,3.2) |   (0.2,2.9)  (1.0,-2.1) |
     *                                                   *

The following is the 1 × 3 process grid:

B,D  |    0    |    1    |    2    
-----| ------- | ------- |------- 
0    |   P00   |   P01   |   P02

Local arrays for Y:

p,q  |                      0                          |                       1
-----|-------------------------------------------------|-----------------------------------------------
     | (29.8,29.7)  (-1.9,1.1)  (-3.0,0.9) (-3.0,-3.8) | (-2.4,-2.8)   (2.0,0.1)  (3.6,-3.4)  (1.4,0.0)
     |  (-3.3,1.0)  (0.4,-2.5)  (4.1,-3.9) (-4.6,-1.4) |  (2.1,-3.5)   (2.8,1.3) (-1.9,-0.1)  (2.3,6.3)
 0   | (-1.7,-1.2)   (0.1,3.4)   (0.9,5.0)  (-0.6,1.6) |  (-1.7,0.7)   (2.8,1.0)   (2.0,1.4) (-0.6,1.0)
     |  (2.3,-0.5)  (1.0,-4.6)   (3.1,0.8)   (1.7,0.4) | (-3.3,-2.2) (-3.2,-5.1)  (-0.3,3.3)  (0.8,0.5)
     |   (3.0,1.9)   (4.5,0.6)  (0.5,-3.8) (-3.0,-0.1) | (-1.5,-3.1)   (1.4,0.1) (-1.3,-1.7) (-2.7,0.7)
     |   (1.0,0.1) (-5.7,-1.9) (-1.4,-1.2)   (0.8,2.6) |   (1.8,0.6)  (-0.7,3.2)   (0.2,2.9) (1.0,-2.1)

There is not any data located on P₀₂.