An optimized parallel algorithm is constructed for solving the nonlinear inverse gravimetry problem of nding several boundary surfaces between layers in multilayered medium. The algorithm is based on the modi ed nonlinear conjugate gradient method with weighting factors. The e cient implementation for GPU was developed. A model problem with synthetic gravitational data was solved.
The problem considered in this paper is nding several interfaces between layers
in a multilayered medium using known gravitational data [
In works [
Here, we construct a parallel algorithm on the basis of this method and implement it for GPUs using the CUDA technology.
We compare the CPU and GPU implementations in terms of computation time in solving the model problem with synthetic gravitational data. ? This work was nancially supported by the Ministry of Education and Science of the Republic of Kazakhstan (project AP 05133873).
We assume that the lower half-space is composed of several layers with constant densities, which are separated by the sought surfaces Sl; l = 1::L, where L is the number of boundary surfaces (see, Fig. 1).
The gravitational eld due to this half-space is equal to the sum of the
gravitational elds due to each surface. Let the boundary surfaces be speci ed
by the functions z = l(x; y), let the density contrasts on them be l, and
let the surfaces have the horizontal asymptotic planes z = Hl. The eld g
produced by the superposition of the boundary surfaces and measured on the
Earth's surface z = 0 is found by the following equation (with accuracy up to a
constant term of summation) [
Z1 Z1 where f is the gravitational constant.
This equation is the Fredholm nonlinear equation of the rst kind of functions l and is ill-posed: it has a nonunique solution which unstably depends on the initial data.
After the discretization of equation (1) on the rectangular grid n = M N , where the right-hand side g(x; y; 0) is given, and the approximation of the left-side integral operator A( 1; :::; l) by the quadrature rules, we obtain the vector of the right-hand side F with the length n, the combined surfaces vector z = [ 1(x1; y1); :::; 1(xM ; yN ); :::; 2(x1; y1); :::; 2(xM ; yN ); :::; L(xM ; yN )] with the length Ln, and the system of nonlinear equations having the form A(z) = F: (2) (3) (4) 3
In this work, to solve problem (2), we use the approach based on the modi ed
linearized conjugate gradient method with weighting factors [
This method has the following form: zk+1 = zk pk = vk +
kpk 1; p0 = v0; k = max hvk; vk vk 1 kvk 1k
i ; 0 v =
S(z); S(z) = A0(z) (A(z)
F ); hpk; S(zk) k kA0(zk) pkk p ; where zk is the solution estimate at kth iteration, is the damping factor, is the vector of weighting factors, A0(z) is the Jacobean matrix of the discretized integral operator A(z), is operation of componentwise vector multiplication.
In this work, we propose the following rules for selection of the weighting factors: [F1; :::; FL] ! [f1; :::; fLn] ! [ 1; :::; Ln]; Fl ! [ n(l 1); n(l 1)+1; :::; nl]; l =
pfi2 + im=1a::xnfpfi2 + g ; 0 < < 1; where is the smoothing parameter.
The elds Fl are extracted from the total eld F using heightwise
transformation technique from [
The condition kA(z) F k = kF k < " for su ciently small " is used as the termination criterion for the iterative process.
In works [
Most expensive part is calculation of the Jacobian matrix A0 and vector A(zk) at each iteration. The matrix and vector are divided into a number of fragments, and each fragment is processed by its own thread. The elements of the Jacobian matrix are calculated on-the- y, which means that the value of an element is computed when calling this element, without storing it in memory.
The adjustment of the kernel execution parameters for the grid size is an
important problem. In previous works [
This imposes some constraints on the input data and GPUs con guration: grid size should be divisible by 128 (128, 256, 512, 1024 ...) and GPUs number should be a power of 2 (1, 2, 4, 8, ...).
To test the constructed modi ed algorithm and to compare it with the unmodied one in terms of execution time, we use the model problem of reconstructing three boundary surfaces on a large grid (512 512 nodes) using the quasi-real data.
The model gravitational eld shown in Fig. 3 was obtained by solving the
forward problem using three surfaces S1; S2; S3 (see, Fig. 4a) with asymprotic
planes H1 = 10km; H2 = 20km; H3 = 30km from work [
The inverse problem was solved using the eight-core Intel Xeon E5-2650 processor and NVIDIA Tesla M2090 GPUs) incorporated in Uran parallel computing system. The Jacobian matrix for this problem is 262144 786432. The threshold parameter = 0:1 was used, which gives the bandwidth parameters 1 = 0:078, 2 = 0:094, 3 = 0:1. The smoothing parameter was = 0:2. For the stopping criterion, " = 0:1 was used. The algorithm took 90 iterations.
The reconstructed surfaces S1; S2; S3 are shown in Fig. 4b. The relative errors
Sl = kSlk are lower than 1%.
The table contains the execution times Tm for various number m of Tesla M2090 GPUs, relative speedup Sm = T1=Tm and e ciency Em = Sm=m. It also contains computation time for eight-core CPU.
The implementation for multiple GPU demonstrates an excellent scaling; the e ciency is more than 90% for eight GPUs. For solving the inverse gravimetry problem of nding several boundary surfaces in multilayered medium, the parallel algorithm was constructed and implemented for multiple GPUs using the CUDA technology. The algorithm is based on the nonlinear conjugate gradient method and on the approach with weighting factors previously proposed by authors. The numerical implementation uses the modi cation based on approximation the Jacobian matrix by dropping out the less signi cant elements and to replacing the matrix by a block-band one.
The model problem of reconstructing three surfaces using the quasi-real gravitational data was solved on a large grid. The GPU implementation reduces the computation time by two orders of magnitude in comparison with CPU. The multi-GPU implementation demonstrates an excellent scaling and nearly 90% e ciency.