Fast Helmholtz Solver : ipar documentation

Lars_S_ · ‎04-11-2017

Hello,

I am calling the MKL shared libraries (2017.2.187) fom the julia language in order to solve a 2D cartesian helmholtz equation (a poisson equation to be more specific). And I think I found some discrepancies in the documentation for the ipar Array (https://software.intel.com/en-us/node/522080).

Here is the julia code:

# julia code: www.julialang.org

# call [...]\mkl\bin\mklvars.bat intel64 ilp64  first to set the environment
  
  MKL_INT = Int32;
  
  r = linspace(0,1,132);
  z = linspace(0,2.5,136);

  nx = length(r);
  ny = length(z);

  in_ax = r[1];
  in_bx = r[nx];
  in_ay = z[1];
  in_by = z[ny];

  in_BCtype    = Array{UInt8}(4);
  in_BCtype[1] = UInt8('D');
  in_BCtype[2] = UInt8('N');
  in_BCtype[3] = UInt8('N');
  in_BCtype[4] = UInt8('N');

  in_q = 0.0; 

  inout_ipar       = Array{MKL_INT}(128);
  inout_ipar[:]    = 123123123; #dummy value in order to identify unaltered entries
  inout_dpar       = Array{Float64}(Int32(5*nx/2+7));
  inout_dpar[:]    = 0.0;
  inout_stat       = 0;
  
  inout_f          = Array{Float64}(Int32((nx+1)*(ny+1)));
  inout_f[:]       = 0.0;
  in_bd_ax         = Array{Float64}(Int32(ny));
  in_bd_ax[:]      = 0.0;
  in_bd_bx         = copy(in_bd_ax);
  in_bd_ay         = Array{Float64}(Int32(nx));
  in_bd_ay[:]      = 0.0;
  in_bd_by         = copy(in_bd_ay);
  inout_xhandle    = Ptr{Any};

  ccall((:D_INIT_HELMHOLTZ_2D, "mkl_rt"),  # function and library
        Ptr{Void},              # ReturnType Void --> no return value
        (Ptr{Float64},  #ax     # ArgumentTypes as a Tuple
        Ptr{Float64},   #bx
        Ptr{Float64},   #ay
        Ptr{Float64},   #by
        Ptr{MKL_INT},   #nx
        Ptr{MKL_INT},   #ny
        Ptr{UInt8},     #BCtype
        Ptr{Float64},   #q
        Ptr{MKL_INT},   #ipar
        Ptr{Float64},   #dpar
        Ptr{MKL_INT},   #stat
        ),
        &in_ax,         #ax     # Arguments passed by reference
        &in_bx,         #bx
        &in_ay,         #ay
        &in_by,         #by
        &nx,            #nx
        &ny,            #ny
        in_BCtype,      #BCtype  
        &in_q,          #q
        inout_ipar,     #ipar  
        inout_dpar,     #dpar
        &inout_stat,    #stat 
        );

  # print ipar array elements in a tabular manner with 0-based indices
  # as presented in the documentation
  println("Values of ipar after INIT call:")
  println("idx\tvalue")
  for i=1:length(inout_ipar)
	if inout_ipar==123123123
		println("$(i-1)\tunset");
	else
		println("$(i-1)\t$(inout_ipar)");
    end
  end
  
  ccall((:D_COMMIT_HELMHOLTZ_2D, "mkl_rt"),
        Ptr{Void},              # ReturnType Void 
        (Ptr{Float64},  #f      # ArgumentTypes as a Tuple
        Ptr{Float64},   #bd_ax
        Ptr{Float64},   #bd_bx
        Ptr{Float64},   #bd_ay
        Ptr{Float64},   #bd_by
        Ptr{Any},       #xhandle
        Ptr{MKL_INT},   #ipar
        Ptr{Float64},   #dpar
        Ptr{MKL_INT},   #stat
        ),
        inout_f,        # Arguments passed by reference
        in_bd_ax,
        in_bd_bx,
        in_bd_ay,
        in_bd_by,
        &inout_xhandle,
        inout_ipar,
        inout_dpar,
        &inout_stat,
        );
  
  
  # print ipar array elements in a tabular manner with 0-based indices
  # as presented in the documentation
  println("Values of ipar after COMMIT call:")
  println("idx\tvalue")
  for i=1:length(inout_ipar)
	if inout_ipar==123123123
		println("$(i-1)\tunset");
	else
		println("$(i-1)\t$(inout_ipar)");
    end
  end

And this is the output:

Values of ipar after INIT call:
idx     value
0       0
1       1
2       1
3       unset
4       unset
5       unset
6       0
7       1
8       1
9       1
10      unset
11      unset
12      132
13      136
14      unset
15      unset
16      unset
17      unset
18      unset
19      unset
20      unset
21      unset
22      unset
23      1
24      1
25      unset
26      unset
27      unset
28      unset
29      unset
30      unset
31      unset
32      unset
33      unset
34      unset
35      unset
36      unset
37      unset
38      unset
39      unset
40      unset
41      unset
42      unset
43      unset
44      unset
45      unset
46      unset
47      unset
48      unset
49      unset
50      unset
51      unset
52      unset
53      unset
54      unset
55      unset
56      unset
57      unset
58      unset
59      unset
60      unset
61      unset
62      unset
63      unset
64      unset
65      unset
66      unset
67      unset
68      unset
69      unset
70      unset
71      unset
72      unset
73      unset
74      unset
75      unset
76      unset
77      unset
78      unset
79      unset
80      unset
81      unset
82      unset
83      unset
84      unset
85      unset
86      unset
87      unset
88      unset
89      unset
90      unset
91      unset
92      unset
93      unset
94      unset
95      unset
96      unset
97      unset
98      unset
99      unset
100     unset
101     unset
102     unset
103     unset
104     unset
105     unset
106     unset
107     unset
108     unset
109     unset
110     unset
111     unset
112     unset
113     unset
114     unset
115     unset
116     unset
117     unset
118     unset
119     unset
120     unset
121     unset
122     unset
123     unset
124     unset
125     unset
126     unset
127     unset

Values of ipar after COMMIT call:
idx     value
0       0
1       1
2       1
3       unset
4       unset
5       unset
6       0
7       1
8       1
9       1
10      unset
11      unset
12      132
13      136
14      unset
15      6
16      138
17      unset
18      unset
19      139
20      337
21      unset
22      unset
23      1
24      1
25      unset
26      unset
27      unset
28      unset
29      unset
30      unset
31      unset
32      unset
33      unset
34      unset
35      unset
36      unset
37      unset
38      unset
39      unset
40      132
41      1
42      1
43      unset
44      unset
45      2
46      0
47      1
48      1
49      1
50      0
51      unset
52      unset
53      unset
54      unset
55      unset
56      unset
57      unset
58      unset
59      unset
60      unset
61      unset
62      unset
63      unset
64      unset
65      unset
66      unset
67      unset
68      unset
69      unset
70      unset
71      unset
72      unset
73      unset
74      unset
75      unset
76      unset
77      unset
78      unset
79      unset
80      unset
81      unset
82      unset
83      unset
84      unset
85      unset
86      unset
87      unset
88      unset
89      unset
90      unset
91      unset
92      unset
93      unset
94      unset
95      unset
96      unset
97      unset
98      unset
99      unset
100     unset
101     unset
102     unset
103     unset
104     unset
105     unset
106     unset
107     unset
108     unset
109     unset
110     unset
111     unset
112     unset
113     unset
114     unset
115     unset
116     unset
117     unset
118     unset
119     unset
120     unset
121     unset
122     unset
123     unset
124     unset
125     unset
126     unset
127     unset

And it seems to me that the following issues arise:

Not only is ipar[3] unused (as stated in the docs) but also ipar[4] and ipar[5], because the values for my selected boundary conditions are stored in ipar[6-9] and not ipar[4-7] as stated in the docs. This increases the index of all following information by 2. From now on, I will use the "new" indices.
The value in ipar[16] equals ipar[15]+ipar[12]+0 and not ipar[15]+ipar[12]+1 as stated in the docs.
In the tabular description of indices 17 to 22 for the 2D Cartesian case:
- ipar[21] is unset, but the docs state that it should be ipar[20]+1. But maybe this will only be set in the 2D_HELMHOLTZ procedure and might not be an issue
- the docs state, that ipar[22] should be ipar[21]+3*ipar[14]/4 , but ipar[14] is unset because nz is not used for the 2D case. So the table entry for the 2D Cartesian case should state "unused" for ipar[22]

Greetings,

Lars

Zhen_Z_Intel · ‎04-11-2017

Hi Lars,

You are using 32bit integer, the interface layer lib linked should be lp64 not ilp4. The ilp64 is used for 64bits integer.

Best regards,
Fiona