A memory e cient algorithm is constructed and implemented for solving the inverse problem of nding a variable density in a curvilinear layer using gravitational data. The algorithm is based on the stabilized biconjugated gradient method. The modi cation based on approximation of the SLAE matrix by Toeplitz-block-Toeplitz one signi cantly reduces memory requirements. The parallel algoritms were implemented for the Uran supercomputer using the hybrid MPI+OpenMP technology. A model problem with synthetic data was solved.
The problem of nding a density in a horizontal or curvilinear layer between given depths was studied in works [1{7]. Solution process consists in two separate steps. The rst one is to extract approximately anomalous e ect of the considered layer from the observed gravitational eld using special technique [3]. This technique was implemented for various multiprocessor systems [5, 6]. The extracted eld is the right part of an integral equation of the desired density. After discretization of the area and approximation of the integral operator, both of these steps are reduced to solving ill-conditioned systems of linear algebraic equations.
In work [7], a parallel algoritm for solving density reconstruction problem was created on the basis of the stabilized biconjugate gradient (BiCGSTAB) method and implemented for multicore processor using the OpenMP technology.
In this work, this algorithm is implemented using the hybrid MPI+OpenMP technology for the Uran supercomputer installed at the Krasovskii Institute of Mathematics and Mechanics.
The memory e cient algoritm for solving this problem is constructed and implemented for multicore processor. The modi cation based on approximation of the SLAE matrix by Toeplitz-block-Toeplitz one signi cantly reduces memory requirements.
Let us assume a Cartesian coordinate system where plane x0y coincides with the Earth's surface, and axis z is directed downwards.
The vertical attraction of gravity g(x; y; 0) generated by the layer = f(x; y; z) : H1(x; y) 6 z 6 H2(x; y)g is given by the following equation [7]: f
Z1 Z1 1 1 p(x p(x x0)2 x0)2
1 (y
1 (y y0)2
H12(x; y) y0)2
H22(x; y)
(x0; y0)dx0dy0 =
g(x; y; 0);
(
Let us assume that the boundary surfaces H1(x; y) and H2(x; y) satisfy the conditions H1(x; y) < H2(x; y), 8(x; y) 2 R2, lim Hi(x; y) = const, and x;y! 1 the density is constant in areas (x; y; z) : H1 6 z 6 H1(x; y) and (x; y; z) : H2(x; y) 6 z 6 H2 , where H1 = min H1(x; y) and H2 = max H2(x; y).
x;y x;y
Problem (
After discretization on a grid n = M N and approximation of the integral
operator using quadrature formulae, expression (
N M f X X u 2 [1::M ]; v 2 [1::N ]: is the regularization parameter.
BiCGSTAB algoritm. For solving system (
k = qk = rk 1 + k(qk 1 !k 1vk 1) k
k 1 k 1 !k
1 0 = 0 = !0 = 1; v0 = q0 = 0 vk = Aqk; sk = rk 1 tk = Ask; zk = zk 1 + !ksk + enddo k = kvk
k hr; vki ! = htk; ski htk; tki k k q ; rk = sk
k k ! t
The initial estimate is z0 = 0. The condition Azk b = kbk < " for some su ciently small " is taken as a termination criterion.
For large ne grids, storing of coe cient matrix A can take large memory space. E.g., for the 29 29 grid, the matrix would be of 218 218 dimension. Full storage of this matrix with double precision requires 512 GB of memory, which is more than one node of the Uran supercomputer has.
For solving the density reconstruction problem using the BiCGSTAB method, a parallel algoritm was developed utilizing the MPI technology. The coe cient matrix A is divided by a number of horizontal bars with respect to the number of used nodes, and each bar is stored and processed on its own node. The hybrid MPI+OpenMP technology was used to parallelize the algoritm for nodes' multicore processors.
Memory e cient algoritm. In work [6], the memory-e cient algoritm for solving the density reconstrustion problem in a horizontal layer was proposed and implemented for multicore processor. This algoritm is based on exploiting the Toeplitz-block-Toeplitz structure of the coe cient matrix.
The elements of matrix can be written in the following form:
ak;p;l;q =
x y
p(xk
p(xk
1
1
xl)2 + (yp
xl)2 + (yp
yq)2 + H12
yq)2 + H22
;
(
Apparently, the matrix elements depend only on the terms (xk xl)2 + (yp yq)2. Note that ak;p;l;q = ak;p+1;l;q+1 and ak;p;l;q = ak+1;p;l+1;q. The former equation means that in each block each descending diagonal from left to right is constant. The latter one means that each block diagonal is constant as well. In other words, the matrix is symmetric Toeplitz-block-Toeplitz.
The matrix structure for a 6 4 grid is shown in Fig. 2. The matrix has 4 4 blocks of 6 6 dimension. Equal elements are marked by the same colors.
The obvious way of storing this matrix is to store the rst row only. Each subsequent row is obtained by double cycle shifting of elements. This method requires storing M N elements, whilst full storage requires M 2 N 2.
In the case of a curvilinear layer, the matrix elements will have the form ak;p;l;q = x y q(xk q(xk 1 1 xl)2 + (yp
yq)2 + H12;(p 1)M+k
xl)2 + (yp
yq)2 + H22;(p 1)M+k
;
(
In this work, we propose the following way to store this matrix approximately.
Note that the farther the matrix element is from the main diagonal, the lower
its value is. We compute the elements of the main block diagonal (i.e., elements
that has p = q) using formula (
MN H1 and H2. The mean values Hw = P Hw;i=M N can be used.
i=1
The number of stored elements will be M (N 2 + 1). E.g, to store the 218 218 matrix, we need only 256 MB instead of 512 GB needed for full storage. Thus, to solve the problem using this algorithm, we need only one node with multicore processor. The algorithm was implemented utilizing the OpenMP technology 4
Let us demonstrate application of the constructed algorithms to solving the model density reconstruction problem on a large grid.
The model problem of reconstructing the density distribution in a curvilinear layer between the depths H1 = 9:5 km and H2 = 11:5 km is considered.
The boundary surfaces are shown in Fig. 3.
Figure 4 shows the model gravitational eld. This eld was obtained by
solving the forward problem (
The inverse problem was solved on 29 29 grid by the BiCGSTAB algorithm and the memory e cient algorithm as described above.
Figure 5b shows the reconstructed density distribution ^. The regularization parameter was = 0:1, " = 0:005 was used as termination criterion.
Solution by the memory e cient algorithm has a relative error k^ k = k k not greater than the error of solution by the unmodi ed BiCGSTAB. Both errors are lower than 0:05.
Table 1 shows the computing times for solving the problem by the BiCGSTAB method using the MPI technology with the Uran supercomputer. Each of supercomputer's nodes has Intel E5-2650 CPU and 192 GB of RAM. One core was used on each node. The coe cient matrix requires 512 GB for full storage. Thus, to solve the problem, we need at least 4 nodes.
The last table line contains the time to solve the problem utilizing the hybrid MPI+OpenMP technology with 8 cores at each node. The computing time reduces to 20 seconds and the e ciency of this technology was 70 %.
The memory e cient algorithm requires 256 MB to store matrix; thus, to solve the problem, we need only one node with multicore processor.
Table 2 shows computing times, parallel speedup, and e ciency for solving the problem by the memory e cient algorithm using the node with two eightcore Intel E5-2650 CPUs. The parallel algoritm show good scaling. 1.5 { 1.9 3.7 7.4
The memory e cient parallel algoritm was constructed on the basis of the BiCGSTAB method for solving the inverse problem of density reconstruction in a curvilinear layer. The modi cation is based on approximation of matrix blocks by the blocks of Toeplitz-block-Toeplitz matrix. Optimized method for matrix storage signi cantly reduces the memory requirements.
The parallel programs were developed for solving the inverse problem of density reconstruction in a curvilinear layer on the Uran supercomputer.
The rst program is based on BiCGSTAB algorithm with full matrix storage using the hybrid MPI+OpenMP technology. The second program is based on the memory e cient algorithm with optimized matrix storage utilizing the OpenMP technology.
The model problem with synthetic data was solved, the speedup and e ciency of parallel algorithms was studied.
This work was partly supported by the Center of Excellence \Geoinformation technologies and geophisical data complex interpretation" of the Ural Federal University Program and by Ural Branch of RAS, project 18-1-1-8.