# Trigonometric Transforms Interface Code Examples

Code presented in this section computes solutions of three simple 1D Helmholtz problems with different boundary conditions: DD, NN and ND cases, where “D” denotes a Dirichlet boundary condition and “N” stands for a Neumann boundary condition. Example “C Example to Solve a Set of 1D Helmholtz Problems” implements the computations in C and Example “Fortran Example to Solve a Set of 1D Helmholtz Problems” provides Fortran 90 code.

The algorithm of computing the solution uses Trigonometric Transform routines. In the DD case, the sine transform is computed, the NN case uses the cosine transform and the ND case corresponds to the staggered cosine transform.

Other details of the Helmholtz problems being solved are printed out along with the computed solutions.

Upon successful execution of Example “C Example to Solve a Set of 1D Helmholtz Problems” the following text is printed out (Example “Fortran Example to Solve a Set of 1D Helmholtz Problems” generates similar output):

`Example of use of MKL Trigonometric Transforms`
`**********************************************`
``
`This example gives the the solutions of the 1D differential problems`
`with the equation -u"+u=f(x), 0<x<1,`
`and with 3 types of boundary conditions:`
`DD case: u(0)=u(1)=0,`
`NN case: u'(0)=u'(1)=0,`
`ND case: u'(0)=u(1)=0.`
`-----------------------------------------------------------------------`
`In general, the error should be of order O(1.0/n**2)`
`For this example, the value of n is 8`
`The approximation error should be of order 5.0e-002 if everything is OK`
`-----------------------------------------------------------------------`
`Note that n should be even to use Trigonometric Transforms !`
`-----------------------------------------------------------------------`
`DOUBLE PRECISION COMPUTATIONS`
`=====================================================================`
``
`The computed solution of DD problem is`
``
`u[0]= 0.000`
`u[1]= 0.153`
`u[2]= 0.524`
`u[3]= 0.895`
`u[4]= 1.049`
`u[5]= 0.895`
`u[6]= 0.524`
`u[7]= 0.153`
`u[8]= 0.000`
``
`Error=4.873e-002`
``
`The computed solution of NN problem is`
``
`u[0]=-0.026`
`u[1]= 0.128`
`u[2]= 0.500`
`u[3]= 0.872`
`u[4]= 1.026`
`u[5]= 0.872`
`u[6]= 0.500`
`u[7]= 0.128`
`u[8]=-0.026`
``
`Error=2.583e-002`
``
`The computed solution of ND problem is`
``
`u[0]=-0.009`
`u[1]= 0.145`
`u[2]= 0.517`
`u[3]= 0.890`
`u[4]= 1.045`
`u[5]= 0.892`
`u[6]= 0.522`
`u[7]= 0.152`
`u[8]= 0.000`
``
`Error=4.470e-002`

C code of the computations is presented below:

## C Example to Solve a Set of 1D Helmholtz Problems

/*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2005-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

!

!*******************************************************************************

! Content:

! Double precision C test example for trigonometric transforms

!*******************************************************************************

!

! This example gives the solution of the 1D differential problems

! with the equation -u"+u=f(x), 0<x<1, and with 3 types of boundary conditions:

! u(0)=u(1)=0 (DD case), or u'(0)=u'(1)=0 (NN case), or u'(0)=u(1)=0 (ND case)

*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include "mkl_trig_transforms.h"

int main(void)

{

MKL_INT n=8, i, k, tt_type;

MKL_INT ir, ipar[128];

double pi=3.14159265358979324, xi, c;

double c1, c2, c3, c4, c5, c6;

double *u, *f, *dpar, *lambda;

DFTI_DESCRIPTOR_HANDLE handle = 0;

/* Printing the header for the example */

printf("\n Example of use of MKL Trigonometric Transforms\n");

printf(" **********************************************\n\n");

printf(" This example gives the the solutions of the 1D differential problems\n");

printf(" with the equation -u\"+u=f(x), 0<x<1, \n");

printf(" and with 3 types of boundary conditions:\n");

printf(" DD case: u(0)=u(1)=0,\n");

printf(" NN case: u'(0)=u'(1)=0,\n");

printf(" ND case: u'(0)=u(1)=0.\n");

printf(" -----------------------------------------------------------------------\n");

printf(" In general, the error should be of order O(1.0/n**2)\n");

printf(" For this example, the value of n is %d\n", n);

printf(" The approximation error should be of order 5.0e-002 if everything is OK\n");

printf(" -----------------------------------------------------------------------\n");

printf(" DOUBLE PRECISION COMPUTATIONS \n");

printf(" =======================================================================\n\n");

u=(double*)malloc((n+1)*sizeof(double));

f=(double*)malloc((n+1)*sizeof(double));

/* NOTE: This example uses shorter dpar array of size 3n/2+2 instead of 5n/2+2

as only sine, cosine, and staggered cosine transforms are used. More details

can be found in section "Partial Differential Equations Support" of the Intel MKL Reference Manual. */

dpar=(double*)malloc((3*n/2+2)*sizeof(double));

lambda=(double*)malloc((n+1)*sizeof(double));

for(i=0;i<=2;i++)

{

/* Varying the type of the transform */

tt_type=i;

/* Computing test solutions u(x) */

for(k=0;k<=n;k++)

{

xi=1.0E0*k/n;

u[k]=pow(sin(pi*xi),2.0E0);

}

/* Computing the right-hand side f(x) */

for(k=0;k<=n;k++)

{

f[k]=(4.0E0*(pi*pi)+1.0E0)*u[k]-2.0E0*(pi*pi);

}

/* Computing the right-hand side for the algebraic system */

for(k=0;k<=n;k++)

{

f[k]=f[k]/(n*n);

}

if (tt_type==0)

{

/* The Dirichlet boundary conditions */

f[0]=0.0E0;

f[n]=0.0E0;

}

if (tt_type==2)

{

/* The mixed Neumann-Dirichlet boundary conditions */

f[n]=0.0E0;

}

/* Computing the eigenvalues for the three-point finite-difference problem */

if (tt_type==0||tt_type==1)

{

for(k=0;k<=n;k++)

{

lambda[k]=pow(2.0E0*sin(0.5E0*pi*k/n),2.0E0)+1.0E0/(n*n);

}

}

if (tt_type==2)

{

for(k=0;k<=n;k++)

{

lambda[k]=pow(2.0E0*sin(0.25E0*pi*(2*k+1)/n),2.0E0)+1.0E0/(n*n);

}

}

/* Computing the solution of 1D problem using trigonometric transforms

First we initialize the transform */

d_init_trig_transform(&n,&tt_type,ipar,dpar,&ir);

if (ir!=0) goto FAILURE;

/* Then we commit the transform. Note that the data in f will be changed

at this stage !

If you want to keep them, save them in some other array

before the call to the routine */

d_commit_trig_transform(f,&handle,ipar,dpar,&ir);

if (ir!=0) goto FAILURE;

/* Now we can apply trigonometric transform */

d_forward_trig_transform(f,&handle,ipar,dpar,&ir);

if (ir!=0) goto FAILURE;

/* Scaling the solution by the eigenvalues */

for(k=0;k<=n;k++)

{

f[k]=f[k]/lambda[k];

}

/* Now we can apply trigonometric transform once again

as ONLY input vector f has changed */

d_backward_trig_transform(f,&handle,ipar,dpar,&ir);

if (ir!=0) goto FAILURE;

/* Cleaning the memory used by handle

Now we can use handle for other kind of trigonometric transform */

free_trig_transform(&handle,ipar,&ir);

if (ir!=0) goto FAILURE;

/* Performing the error analysis */

c1=0.0E0;

c2=0.0E0;

c3=0.0E0;

for(k=0;k<=n;k++)

{

/* Computing the absolute value of the exact solution */

c4=fabs(u[k]);

/* Computing the absolute value of the computed solution

Note that the solution is now in place of the former right-hand side ! */

c5=fabs(f[k]);

/* Computing the absolute error */

c6=fabs(f[k]-u[k]);

/* Computing the maximum among the above 3 values c4-c6 */

if (c4>c1) c1=c4;

if (c5>c2) c2=c5;

if (c6>c3) c3=c6;

}

/* Printing the results */

if (tt_type==0)

{

printf("The computed solution of DD problem is\n\n");

}

if (tt_type==1)

{

printf("The computed solution of NN problem is\n\n");

}

if (tt_type==2)

{

printf("The computed solution of ND problem is\n\n");

}

for(k=0;k<=n;k++)

{

printf("u[%d]=",k);

printf("%6.3f\n",f[k]);

}

printf("\nError=%6.3e\n\n",c3/c1);

if (c3/c1>0.05)

{

printf("The computed solution seems to be inaccurate. ");

goto FAILURE;

}

/* End of the loop over the different kind of transforms and problems */

}

/* Success message to print if everything is OK */

printf("This example has successfully PASSED through all steps of computation!\n");

/* Jumping over failure message */

goto SUCCESS;

/* Failure message to print if something went wrong */

FAILURE: printf("Failed to compute the solution(s)...\n");

SUCCESS: return 0;

/* End of the example code */

}

Fortran 90 code of the computations is presented below:

## Fortran Example to Solve a Set of 1D Helmholtz Problems

!*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2005-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

!

!*******************************************************************************

! Content:

! Double precision Fortran90 test example for trigonometric transforms

!*******************************************************************************

!

! This example gives the solution of the 1D differential problems

! with the equation -u"+u=f(x), 0<x<1, and with 3 types of boundary conditions:

! u(0)=u(1)=0 (DD case), or u'(0)=u'(1)=0 (NN case), or u'(0)=u(1)=0 (ND case)

program d_tt_example_bvp

use mkl_dfti

use mkl_trig_transforms

implicit none

integer n, i, k,j, tt_type

integer ir, ipar(128)

parameter (n=8)

double precision pi, xi

double precision c1, c2, c3, c4, c5, c6

! NOTE: This example uses shorter dpar array of size 3n/2+2 instead of 5n/2+2

! as only sine, cosine, and staggered cosine transforms are used. More details

! can be found in section "Partial Differential Equations Support" of the Intel MKL Reference Manual. */

double precision u(n+1), f(n+1), dpar(3*n/2+2), lambda(n+1)

parameter (pi=3.14159265358979324D0)

type(dfti_descriptor), pointer :: handle

! Printing the header for the example

print *, ''

print *, ' Example of use of MKL Trigonometric Transforms'

print *, ' **********************************************'

print *, ''

print *, ' This example gives the solution of the 1D differential problems'

print *, ' with the equation -u"+u=f(x), 0<x<1, '

print *, ' and with 3 types of boundary conditions:'

print *, ' DD case: u(0)=u(1)=0,'

print *, ' NN case: u''(0)=u''(1)=0,'

print *, ' ND case: u''(0)=u(1)=0.'

print *, ' -----------------------------------------------------------------------'

print *, ' In general, the error should be of order O(1.0/n**2)'

print *, ' For this example, the value of n is', n

print *, ' The approximation error should be of order 0.5E-01, if everything is OK'

print *, ' -----------------------------------------------------------------------'

print *, ' DOUBLE PRECISION COMPUTATIONS '

print *, ' ======================================================================='

print *, ''

do i=0,2

! Varying the type of the transform

tt_type=i

! Computing test solution u(x)

do k=1,n+1

xi=1.0D0*(k-1)/n

u(k)=dsin(pi*xi)**2

end do

! Computing the right-hand side f(x)

do k=1,n+1

f(k)=(4.0D0*(pi**2)+1.0D0)*u(k)-2.0D0*(pi**2)

end do

! Computing the right-hand side for the algebraic system

do k=1,n+1

f(k)=f(k)/(n**2)

end do

if (tt_type.eq.0) then

! The Dirichlet boundary conditions

f(1)=0.0D0

f(n+1)=0.0D0

end if

if (tt_type.eq.2) then

! The mixed Neumann-Dirichlet boundary conditions

f(n+1)=0.0D0

end if

! Computing the eigenvalues for the three-point finite-difference problem

if (tt_type.eq.0.or.tt_type.eq.1) then

do k=1,n+1

lambda(k)=(2.0D0*dsin(0.5D0*pi*(k-1)/n))**2+1.0D0/(n**2)

end do

end if

if (tt_type.eq.2) then

do k=1,n+1

lambda(k)=(2.0D0*dsin(0.25D0*pi*(2*k-1)/n))**2+1.0D0/(n**2)

end do

end if

! Computing the solution of 1D problem using trigonometric transforms

! First we initialize the transform

CALL D_INIT_TRIG_TRANSFORM(n,tt_type,ipar,dpar,ir)

if (ir.ne.0) goto 99

! Then we commit the transform. Note that the data in f will be changed at this stage !

! If you want to keep them, save them in some other array before the call to the routine

CALL D_COMMIT_TRIG_TRANSFORM(f,handle,ipar,dpar,ir)

if (ir.ne.0) goto 99

! Now we can apply trigonometric transform

CALL D_FORWARD_TRIG_TRANSFORM(f,handle,ipar,dpar,ir)

if (ir.ne.0) goto 99

! Scaling the solution by the eigenvalues

do k=1,n+1

f(k)=f(k)/lambda(k)

end do

! Now we can apply trigonometric transform once again as ONLY input vector f has changed

CALL D_BACKWARD_TRIG_TRANSFORM(f,handle,ipar,dpar,ir)

if (ir.ne.0) goto 99

! Cleaning the memory used by handle

! Now we can use handle for other KIND of trigonometric transform

CALL FREE_TRIG_TRANSFORM(handle,ipar,ir)

if (ir.ne.0) goto 99

! Performing the error analysis

c1=0.0D0

c2=0.0D0

c3=0.0D0

do k=1,n+1

! Computing the absolute value of the exact solution

c4=dabs(u(k))

! Computing the absolute value of the computed solution

! Note that the solution is now in place of the former right-hand side !

c5=dabs(f(k))

! Computing the absolute error

c6=dabs(f(k)-u(k))

! Computing the maximum among the above 3 values c4-c6

if (c4.gt.c1) c1=c4

if (c5.gt.c2) c2=c5

if (c6.gt.c3) c3=c6

end do

! Printing the results

if (tt_type.eq.0) then

print *, 'The computed solution of DD problem is'

end if

if (tt_type.eq.1) then

print *, 'The computed solution of NN problem is'

end if

if (tt_type.eq.2) then

print *, 'The computed solution of ND problem is'

endif

print *, ''

do k=1,n+1

write(*,11) k,f(k)

end do

print *, ''

write(*,12) c3/c1

print *, ''

if (c3/c1.ge.5.0D-2) then

print *, 'The computed solution seems to be inaccurate.'

goto 99

endif

! End of the loop over the different kind of transforms and problems

end do

! Success message to print if everything is OK

print *, 'This example has successfully PASSED through all steps of computation!'

! Jumping over failure message

goto 1

! Failure message to print if something went wrong

99 continue

print *, 'FAILED to compute the solution(s)...'

1 continue

! Print formats

11 format(1x,'u(',I1,')=',F6.3)

12 format(1x,'Relative error =',E10.3)

! End of the example code

end

Parent topic: PDE Support Code Examples