NAG FL Interface
d03faf (dim3_​ellip_​helmholtz)

1 Purpose

d03faf solves the Helmholtz equation in Cartesian coordinates in three dimensions using the standard seven-point finite difference approximation. This routine is designed to be particularly efficient on vector processors.

2 Specification

Fortran Interface
Subroutine d03faf ( xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, ldf, ldf2, f, pertrb, w, lwrk, ifail)
Integer, Intent (In) :: l, lbdcnd, m, mbdcnd, n, nbdcnd, ldf, ldf2, lwrk
Integer, Intent (Inout) :: ifail
Real (Kind=nag_wp), Intent (In) :: xs, xf, bdxs(ldf2,n+1), bdxf(ldf2,n+1), ys, yf, bdys(ldf,n+1), bdyf(ldf,n+1), zs, zf, bdzs(ldf,m+1), bdzf(ldf,m+1), lambda
Real (Kind=nag_wp), Intent (Inout) :: f(ldf,ldf2,n+1)
Real (Kind=nag_wp), Intent (Out) :: pertrb, w(lwrk)
C Header Interface
#include <nag.h>
void  d03faf_ (const double *xs, const double *xf, const Integer *l, const Integer *lbdcnd, const double bdxs[], const double bdxf[], const double *ys, const double *yf, const Integer *m, const Integer *mbdcnd, const double bdys[], const double bdyf[], const double *zs, const double *zf, const Integer *n, const Integer *nbdcnd, const double bdzs[], const double bdzf[], const double *lambda, const Integer *ldf, const Integer *ldf2, double f[], double *pertrb, double w[], const Integer *lwrk, Integer *ifail)
The routine may be called by the names d03faf or nagf_pde_dim3_ellip_helmholtz.

3 Description

d03faf solves the three-dimensional Helmholtz equation in Cartesian coordinates:
2u x2 + 2u y2 + 2u z2 +λu=fx,y,z.  
This subroutine forms the system of linear equations resulting from the standard seven-point finite difference equations, and then solves the system using a method based on the fast Fourier transform (FFT) described by Swarztrauber (1984). This subroutine is based on the routine HW3CRT from FISHPACK (see Swarztrauber and Sweet (1979)).
More precisely, the routine replaces all the second derivatives by second-order central difference approximations, resulting in a block tridiagonal system of linear equations. The equations are modified to allow for the prescribed boundary conditions. Either the solution or the derivative of the solution may be specified on any of the boundaries, or the solution may be specified to be periodic in any of the three dimensions. By taking the discrete Fourier transform in the x- and y-directions, the equations are reduced to sets of tridiagonal systems of equations. The Fourier transforms required are computed using the multiple FFT routines found in Chapter C06.

4 References

Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America
Swarztrauber P N and Sweet R A (1979) Efficient Fortran subprograms for the solution of separable elliptic partial differential equations ACM Trans. Math. Software 5 352–364

5 Arguments

1: xs Real (Kind=nag_wp) Input
On entry: the lower bound of the range of x, i.e., xsxxf.
Constraint: xs<xf.
2: xf Real (Kind=nag_wp) Input
On entry: the upper bound of the range of x, i.e., xsxxf.
Constraint: xs<xf.
3: l Integer Input
On entry: the number of panels into which the interval xs,xf is subdivided. Hence, there will be l+1 grid points in the x-direction given by xi=xs+i-1×δx, for i=1,2,,l+1, where δx=xf-xs/l is the panel width.
Constraint: l5.
4: lbdcnd Integer Input
On entry: indicates the type of boundary conditions at x=xs and x=xf.
lbdcnd=0
If the solution is periodic in x, i.e., uxs,y,z=uxf,y,z.
lbdcnd=1
If the solution is specified at x=xs and x=xf.
lbdcnd=2
If the solution is specified at x=xs and the derivative of the solution with respect to x is specified at x=xf.
lbdcnd=3
If the derivative of the solution with respect to x is specified at x=xs and x=xf.
lbdcnd=4
If the derivative of the solution with respect to x is specified at x=xs and the solution is specified at x=xf.
Constraint: 0lbdcnd4.
5: bdxsldf2n+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to x at x=xs. When lbdcnd=3 or 4, bdxsjk=uxxs,yj,zk, for j=1,2,,m+1 and k=1,2,,n+1.
When lbdcnd has any other value, bdxs is not referenced.
6: bdxfldf2n+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to x at x=xf. When lbdcnd=2 or 3, bdxfjk=uxxf,yj,zk, for j=1,2,,m+1 and k=1,2,,n+1.
When lbdcnd has any other value, bdxf is not referenced.
7: ys Real (Kind=nag_wp) Input
On entry: the lower bound of the range of y, i.e., ysyyf.
Constraint: ys<yf.
8: yf Real (Kind=nag_wp) Input
On entry: the upper bound of the range of y, i.e., ysyyf.
Constraint: ys<yf.
9: m Integer Input
On entry: the number of panels into which the interval ys,yf is subdivided. Hence, there will be m+1 grid points in the y-direction given by yj=ys+j-1×δy, for j=1,2,,m+1, where δy=yf-ys/m is the panel width.
Constraint: m5.
10: mbdcnd Integer Input
On entry: indicates the type of boundary conditions at y=ys and y=yf.
mbdcnd=0
If the solution is periodic in y, i.e., ux,yf,z=ux,ys,z.
mbdcnd=1
If the solution is specified at y=ys and y=yf.
mbdcnd=2
If the solution is specified at y=ys and the derivative of the solution with respect to y is specified at y=yf.
mbdcnd=3
If the derivative of the solution with respect to y is specified at y=ys and y=yf.
mbdcnd=4
If the derivative of the solution with respect to y is specified at y=ys and the solution is specified at y=yf.
Constraint: 0mbdcnd4.
11: bdysldfn+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to y at y=ys. When mbdcnd=3 or 4, bdysik=uyxi,ys,zk, for i=1,2,,l+1 and k=1,2,,n+1.
When mbdcnd has any other value, bdys is not referenced.
12: bdyfldfn+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to y at y=yf. When mbdcnd=2 or 3, bdyfik=uyxi,yf,zk, for i=1,2,,l+1 and k=1,2,,n+1.
When mbdcnd has any other value, bdyf is not referenced.
13: zs Real (Kind=nag_wp) Input
On entry: the lower bound of the range of z, i.e., zszzf.
Constraint: zs<zf.
14: zf Real (Kind=nag_wp) Input
On entry: the upper bound of the range of z, i.e., zszzf.
Constraint: zs<zf.
15: n Integer Input
On entry: the number of panels into which the interval zs,zf is subdivided. Hence, there will be n+1 grid points in the z-direction given by zk=zs+k-1×δz, for k=1,2,,n+1, where δz=zf-zs/n is the panel width.
Constraint: n5.
16: nbdcnd Integer Input
On entry: specifies the type of boundary conditions at z=zs and z=zf.
nbdcnd=0
if the solution is periodic in z, i.e., ux,y,zf=ux,y,zs.
nbdcnd=1
if the solution is specified at z=zs and z=zf.
nbdcnd=2
if the solution is specified at z=zs and the derivative of the solution with respect to z is specified at z=zf.
nbdcnd=3
if the derivative of the solution with respect to z is specified at z=zs and z=zf.
nbdcnd=4
if the derivative of the solution with respect to z is specified at z=zs and the solution is specified at z=zf.
Constraint: 0nbdcnd4.
17: bdzsldfm+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to z at z=zs. When nbdcnd=3 or 4, bdzsij=uzxi,yj,zs, for i=1,2,,l+1 and j=1,2,,m+1.
When nbdcnd has any other value, bdzs is not referenced.
18: bdzfldfm+1 Real (Kind=nag_wp) array Input
On entry: the values of the derivative of the solution with respect to z at z=zf. When nbdcnd=2 or 3, bdzfij=uzxi,yj,zf, for i=1,2,,l+1 and j=1,2,,m+1.
When nbdcnd has any other value, bdzf is not referenced.
19: lambda Real (Kind=nag_wp) Input
On entry: the constant λ in the Helmholtz equation. For certain positive values of λ a solution to the differential equation may not exist, and close to these values the solution of the discretized problem will be extremely ill-conditioned. If λ>0, d03faf will set ifail=3, but will still attempt to find a solution. However, since in general the values of λ for which no solution exists cannot be predicted a priori, you are advised to treat any results computed with λ>0 with great caution.
20: ldf Integer Input
On entry: the first dimension of the arrays f, bdys, bdyf, bdzs and bdzf as declared in the (sub)program from which d03faf is called.
Constraint: ldfl+1.
21: ldf2 Integer Input
On entry: the second dimension of the array f and the first dimension of the arrays bdxs and bdxf as declared in the (sub)program from which d03faf is called.
Constraint: ldf2m+1.
22: fldfldf2n+1 Real (Kind=nag_wp) array Input/Output
On entry: the values of the right-side of the Helmholtz equation and boundary values (if any).
fijk = fxi,yj,zki=2,3,,l, j=2,3,,m ​ and ​ k=2,3,,n.  
On the boundaries f is defined by
lbdcnd f1jk fl+1jk  
0 fxs,yj,zk fxs,yj,zk  
1 uxs,yj,zk uxf,yj,zk  
2 uxs,yj,zk fxf,yj,zk j=1,2,,m+1
3 fxs,yj,zk fxf,yj,zk k=1,2,,n+1
4 fxs,yj,zk uxf,yj,zk  
       
mbdcnd fi1k fim+1k  
0 fxi,ys,zk fxi,ys,zk  
1 uxi,ys,zk uxi,yf,zk  
2 uxi,ys,zk fxi,yf,zk i=1,2,,l+1
3 fxi,ys,zk fxi,yf,zk k=1,2,,n+1
4 fxi,ys,zk uxi,yf,zk  
       
nbdcnd fij1 fijn+1  
0 fxi,yj,zs fxi,yj,zs  
1 uxi,yj,zs uxi,yj,zf  
2 uxi,yj,zs fxi,yj,zf i=1,2,,l+1
3 fxi,yj,zs fxi,yj,zf j=1,2,,m+1
4 fxi,yj,zs uxi,yj,zf  
Note: if the table calls for both the solution u and the right-hand side f on a boundary, the solution must be specified.
On exit: contains the solution ui,j,k of the finite difference approximation for the grid point xi,yj,zk, for i=1,2,,l+1, j=1,2,,m+1 and k=1,2,,n+1.
23: pertrb Real (Kind=nag_wp) Output
On exit: pertrb=0, unless a solution to Poisson's equation λ=0 is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and nbdcnd=0 or 3). In this case a solution may not exist. pertrb is a constant, calculated and subtracted from the array f, which ensures that a solution exists. d03faf then computes this solution, which is a least squares solution to the original approximation. This solution is not unique and is unnormalized. The value of pertrb should be small compared to the right-hand side f, otherwise a solution has been obtained to an essentially different problem. This comparison should always be made to ensure that a meaningful solution has been obtained.
24: wlwrk Real (Kind=nag_wp) array Output
25: lwrk Integer Input
On entry: the dimension of the array w as declared in the (sub)program from which d03faf is called. w is no longer used as workspace, lwrk can therefore be safely set to zero.
26: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, l=value.
Constraint: l5.
On entry, lbdcnd=value.
Constraint: lbdcnd0 and lbdcnd4.
On entry, ldf2=value and m=value.
Constraint: ldf2m+1.
On entry, ldf=value and l=value.
Constraint: ldfl+1.
On entry, m=value.
Constraint: m5.
On entry, mbdcnd=value.
Constraint: mbdcnd0 and mbdcnd4.
On entry, n=value.
Constraint: n5.
On entry, nbdcnd=value.
Constraint: nbdcnd0 and nbdcnd4.
On entry, xs=value and xf=value.
Constraint: xs<xf.
On entry, ys=value and yf=value.
Constraint: ys<yf
On entry, zs=value and zf=value.
Constraint: zs<zf.
ifail=3
λ>0.0 in Helmholtz's equation – a solution may not exist.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

Not applicable.

8 Parallelism and Performance

d03faf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d03faf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The execution time is roughly proportional to l×m×n× log2l + log2m +5 , but also depends on input arguments lbdcnd and mbdcnd.

10 Example

This example solves the Helmholz equation
2u x2 + 2u y2 + 2u z2 +λ u=fx,y,z  
for x,y,z0,1×0,2π×0,π2, where λ=-2, and fx,y,z is derived from the exact solution
ux,y,z=x4sinycosz.  
The equation is subject to the following boundary conditions, again derived from the exact solution given above.

10.1 Program Text

Program Text (d03fafe.f90)

10.2 Program Data

Program Data (d03fafe.d)

10.3 Program Results

Program Results (d03fafe.r)