Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2002/ps/art1_1.ps Äàòà èçìåíåíèÿ: Mon Dec 16 17:39:35 2002 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:59:35 2012 Êîäèðîâêà:

Numerical Methods and Programming, 2002, Vol. 3 1
UDC 519.6
SOME FEATURES OF SOLVING THE PROBLEMS OF MAPPING FOR ALLOCATION OF
CHEMICAL ELEMENTS ON STELLAR SURFACES AS ILL­POSED PROBLEMS WITH
THE USE OF MULTIPROCESSOR SYSTEMS
A. G. Yagola 1 , V. N. Titarenko 1 , M. P. Vasiliev 1 , and E. V. Shimanovskaya 1
The problem of mapping for allocation of chemical elements on stellar surfaces is considered as
an ill­posed one. The minimization of Tikhonov's functional with the choice of a regularization
parameter according to the finite­dimensional generalized discrepancy principle is used to construct a
regularizing algorithm for solving the problem under consideration. The method based on projections
of conjugate gradients onto a set of non­negative vectors is chosen as a method of minimization. An
approach to numerical solution of the minimization problem is considered for the above mapping
problem. A multiprocessor computer is used to solve several model problems numerically. Some
schemes of parallelization are offered and special features of their realization are discussed.
1. Introduction. Many stars manifest periodic variations in the absorption lines profiles of some chemical
elements. Such variations may be explained by the fact that chemical elements possess nonuniform distributions
on surfaces of rotating stars and, hence, by the difference in the absorption lines profiles in various regions on the
stellar surfaces. Spectral structure of continuous radiation for a star and its brightness, which are mainly defined
by physical conditions, are constant or have very small variations. So the changes of profiles are produced by
heterogeneous chemical compositions rather than by different physical conditions on the stellar surface. This
statement may be checked by the fact that absorption lines profiles for an element with different ionization
degrees vary synchronously and profiles for different elements with close excitation and ionization potentials
vary asynchronously.
Problems to reconstruct a magnetic field distribution on various stars arise frequently in modern astro­
physics. Quite often, the stars described above possess essential magnetic fields varying with the same periods
as the line profiles. Thereby the heterogeneity is most likely produced by a strong magnetic field of a star. This
property adds interest to the reconstruction of the stellar chemical composition distribution for a subsequent
comparison with the distribution of the magnetic field.
2. The problem statement. Let's consider the problem of the absorption lines profiles in the stellar
spectrum observed from the Earth. The problem is set as in [1, 2]. By definition the energy dE being transferred
by radiation across a unit surface area at an angle ` to its normal in the solid angle 2ú sin ` d` in the wavelength
band [–; – + d–] during the period [t; t + dt] is
dE = 2úI(M; `; –) sin ` d` cos ` d– dt; (2:1)
where I(M; `; –) is a radiation intensity at a point M on the stellar surface. The energy coming to the Earth
from the entire visible surface of a star is proportional to the quantity F +
– =
RR
I(M; `; –) cos ` dM , where dM
is a surface element on the sphere, ` is the angle between the normal to the stellar surface at the point M and
the viewing line. The integration is accomplished over the entire visible surface of the star. Now we introduce
the spherical system of coordinates (`; ') whose axis is directed along the viewing line (see Fig. 1). So the unit
surface area is dM = sin ` d` d', and the visible surface is chosen so that the condition cos ` ? 0 holds.
It should be taken into account that for stars rotating rapidly enough a Doppler's shift of the radiation may
differ at the different points on the stellar surface. This shift may attain values comparable with the line width.
Thereby the argument – of the function I(M; `; –) in the integrand for F +
– should be replaced by –+4–D (M ),
where 4–D (M ) is the Doppler's shift of radiation.
We should distinguish the radiation intensity in the continuum I 0 (M; `; –) and the radiation intensity in
the line. Because the considered line widths are rather small, the continuum intensity is practically constant.
In addition, we may assume that the continuum intensity is irrelevant to a particular point M on the stellar
surface, since physical conditions on the star are approximately homogeneous. In accordance with what has
been said, we might agree with I 0 (M; `; –) = I 0 (`).
1 M.V. Lomonosov Moscow State University, Faculty of Physics, Vorobyevy Gory, Moscow, 119899; e­mail:
yagola@inverse.phys.msu.su; ill­posed@mail.ru
c
fl Research Computing Center, Moscow State University, 119992 Moscow, Russian Federation

2 Numerical Methods and Programming, 2002, Vol. 3
M
N
O
Axis of
gyration
Line of view
Equator
'


Fig. 1. The system of coordinates on the rotating star
The dependence of the continuum intensity on the angle (known as the limb darkening) takes the form
I 0 (`) = I 0 u 1 (`), where u 1 (`) is a linear function of cos `. We assume that u 1 (`) = 1 \Gamma u 1 + u 1 cos `. The
dependence of the coefficient u 1 upon the temperature is carefully analyzed, the stellar temperature may be
recovered from its spectral class.
Investigations revealed that function
r(M; –) = 1 \Gamma I(M; 0; –)
I 0 (M; 0) (2:2)
can always be expressed more precisely by appeal to Minnaerth's formula
r(M; –) =
` 1
x(M;–) + 1
R c (M )
' \Gamma1
; (2:3)
where R c (M ) is the central depth of the line, x(M;–) = z(M )k(–), k(–) is so­called Foygt profile. The coefficient
z(M ) is defined by the relative content of the absorbing element in the atmosphere.
The dependence of the radiation intensity in the line I(M; `; –) upon the angle ` becomes rather compli­
cated. In the paper we assume that the appropriate dependence may be written as
1 \Gamma I(M; `; –)
I 0 (M; `)
= r(M; –)u 2 (`); (2:4)
where u 2 (`) is a linear function with respect to cos `: u 2 (`) = 1 \Gamma u 2 + u 2 cos `. The parameter u 2 may be
found, for instance, by numerical integrating of the transfer equations for the atmosphere model of the star.
The value of the parameter u 2 is often taken as zero.
When physical conditions happen to be homogeneous, the parameters R c and z involved in Minnaerth's
formula (2.3) become dependent. The variables R c and z correspond to the central depth and width of the line.
It has been shown that the dependence of the central depth on the value z may be approximated as follows:
R c = h(z) = C 1
\Gamma 1 \Gamma e \GammaC 2 z
\Delta ; (2:5)
where the parameters C 1 , C 2 are suitable chosen on the basis of the investigation of the atmosphere model.
Finally, the absorption lines profiles in the spectrum of the radiation for the whole stellar disk observed
from the Earth may be determined by the two­dimensional integral equation which is nonlinear relative to z(M ):
r(–) = 1
A
ZZ
cos `?0
f
\Gamma
z(M ); – + 4–D (M )
\Delta
u 1 (`)u 2 (`) cos ` dM; (2:6)
where
A = 2ú
ú=2
Z
0
cos ` u 1 (`) sin ` d`; f(z; –) =
` 1
h(z)
+ 1
zk(–)
' \Gamma1
: (2:7)

Numerical Methods and Programming, 2002, Vol. 3 3
In the process of the star's rotation the profiles R(–) may vary. Thereby, if t is the time and ! is the
angular rotation velocity, equation (2.6) reduces to
R(–; !t) =
ZZ
f
\Gamma z(M ); – + 4–D (M; !t)
\Delta H(M;!t)dM: (2:8)
If we regard the latitude ' and the longitude / as the coordinates on the stellar surface (see Fig. 1), then a
surface element is equal to dM = cos 'd'd/,
4–D ('; /; !t) = – 0
c
v eq sin fi cos ' cos(/ \Gamma !t); (2:9)
H('; /; !t) =
8 !
:
1
A u 1 (`)u 2 (`) cos ` cos '; cos ` ? 0;
0; cos ` 6 0;
(2:10)
cos ` = cos fi sin ' + cos ' sin fi cos(/ \Gamma !t): (2:11)
Here v eq is the equatorial velocity of the star's rotation and fi is the angle of the slop of the star's rotation axis
to the viewing line.
The Foygt profile k(–) characterizing the dependence of the profile upon the wavelength is a convolution
of two functions describing the profile of the line widened by the pressure (collisions) and a Doppler's shift in
the process of microturbulent motions in the stellar atmosphere:
k(–) = ae l
úae d
+1 Z
\Gamma1
e \Gammat 2
dt
ae 2
l
ae 2
d
+
`
t \Gamma –
ae d
' 2 : (2:12)
Here ae l is the quantity of the Lorentz widening by the pressure, the parameter ae d characterizes the line widening
by microturbulent motions.
We will recover the function z(M ) from the profiles of the line measured at different phases which are
represented by the function R(–; !t). The restriction z(M ) ? 0 on the function z(M ) arises naturally. The
problem was formulated in such a form for the first time in [3].
3. A regularization technique. Nonlinear integral equation (2.8) is the first­kind Fredholm equation,
so it is an ill­posed problem. Therefore, regularization algorithms should be constructed to solve it. Consider
algorithms based on minimization of Tikhonov's smoothing functional. Tikhonov's functional for equation (2.8)
is:
M ff [z] = I[z] + ff\Omega\Gamma z] j
ZZ
d– d(!t)
ae
R(–; !t) \Gamma
ZZ
f
\Gamma
z(M ); – + 4–D (M; !t)
\Delta
H(M;!t)dM
oe 2
+
+ ff
ZZ (
z 2 +
` @z
@'
' 2
+
` @z
@/
' 2
)
d'd/:
(3:1)
Assume that the function z(M ) belongs to the space W 1
2 on the stellar surface and it's required to minimize
the functional on the set of all nonnegative functions z(M ) ? 0. The function z(M ) is sought on the set
('; /) 2 [\Gammafi; ú
2 ] \Theta [0; 2ú], because the star's regions with ' 2 [\Gamma ú
2 ; \Gammafi] are not visible from the Earth. The
function R(–; !t) representing the experimental data is supposed to be known for – 2
h
– 0 \Gamma 4–
2 ; – 0 + 4–
2
i
and
t 2 [0; T ].
It's necessary to consider the finite difference approximation of the problem to solve it numerically. Let's
introduce a rectangular grid f' i ; / j g on the stellar surface that is the Cartesian product of grids f' i g n
1 , f/ j g m
1
so that \Gammafi 6 ' 1 ! ' 2 ! : : : ! 'n 6 ú
2 , 0 6 / 1 ! / 2 ! : : : ! /m 6 2ú. The grids f' i g n
1 , f/ j g m
1 are equidistant
and they are defined by formulas ' i = \Gammafi + h'
2 + h' (i \Gamma 1), i = 1; n; / j = h/ (j \Gamma 1), j = 1; m, where steps
h' = ú=2+fi
n , h/ = 2ú
m . We will use the notation z ij = z(' i ; / j ), “
z is the vector z ij , i = 1; n, j = 1; m. We
introduce equidistant grids with respect to the wavelength – and the phase !t: f– k g N–
1 and f(!t) l g N!t
1 , so that
h – = 4–
N– \Gamma 1 , – k = – 0 \Gamma 4–
2 + h– (k \Gamma 1), k = 1; N– ; h!t = !T
N!t \Gamma 1 , (!t) l = h!t (l \Gamma 1), l = 1; N!t . The
integrals in (3.1) are approximated by the rectangles formula. Thus, one of the finite difference approximations
to the smoothing functional is:
c
M ff [“z] = “
I [“z] + ff
b\Omega [“z]; (3:2)

4 Numerical Methods and Programming, 2002, Vol. 3
where “
I[“z] j g( “
J [“z]), g(x) j x and
“
J [“z] =
N–
X
k=1
N!t X
l=1
W kl
n
R kl \Gamma
n
X
i=1
m
X
j=1
h'h/H
\Gamma ' i ; / j ; (!t) l
\Delta f
\Gamma z ij ; – k + 4–D (' i ; / j ; (!t) l )
\Delta o 2
; (3:3)
b
\Omega [“z] = h'h/
n
X
i=1
m
X
j=1
z 2
ij + h/
h'
n
X
i=2
m
X
j=1
(z ij \Gamma z i\Gamma1;j ) 2 + h'
h/
n
X
i=1
n m
X
j=2
(z ij \Gamma z i;j \Gamma1 ) 2 + (z i1 \Gamma z im ) 2
o
: (3:4)
Here W kl ? 0 are the weights introduced for the following reasons. The accuracy in measuring R(–; !t) and the
number of measurements of R(–; !t) may differ in different phases. The profiles corresponding to each phase
are usually known with a high degree of accuracy in the center of the line and with a lower degree of accuracy
on the wings. For the concordance with infinite dimensional problem we assume W kl ¸ h–h!t (in all models
considered below it's assumed that W kl = h–h!t ). So, the initial problem leads to the minimization problem
on the N­dimensional set b
D of nonnegative vectors z ij , where N j n \Theta m, followed by a proper choice of the
regularization parameter ff in agreement with some principle.
The choice of functionals b
\Omega\Gamma “
z], “
J [“z], g(x) in such a form ensures the agreement with all conditions (see [2])
to construct regularization algorithms for an ill­posed problem.
4. The choice of the regularization parameter. We assume that the input data inaccuracy ffi depends
on the inaccuracy of determination of stellar absorption lines profiles only.
ffi 2 j
ZZ
d– d(!t)
\Gamma R(–; !t) \Gamma R ffi (–; !t)
\Delta 2
; (4:1)
where R(–; !t) is an exact line profile, R ffi (–; !t) is an experimentally measured profile. Then instead of the
exact functional J [z] we have an approximate functional J ffi [z] which satisfies the condition of the approximation:
fi fi J [z] \Gamma J ffi [z]
fi fi 6 \Psi 0
\Gamma
ffi; \Omega\Gamma z]
\Delta
: (4:2)
Instead of the initial functional J ffi [z] we deal with its finite­dimensional approximation “
J ffi [“z] which satisfies
the following condition of the approximation:
fi fi “
J ffi [“z] \Gamma J ffi [z]
fi fi 6 b
\Psi
\Gamma j; \Omega\Gamma z]
\Delta : (4:3)
The vector of the approximation j is of the form j j (ffi; 1=N; 1= e
N), where N = n \Theta m, e
N = N– \Theta N!t . It's
clear that j ! 0 as ffi ! 0, N; e
N !1.
Note that approximation measures \Psi 0 (ffi; \Omega\Gamma z]) and b
\Psi(j; z]) should obey specific properties from [2], which
are required to construct a regularization algorithm. The corollary to conditions of approximation (4.2) and
(4.3) is
fi fi J [z] \Gamma “
J j [“z]
fi fi 6 \Psi 0
\Gamma
ffi; \Omega\Gamma z]
\Delta
+ b
\Psi
\Gamma
j; \Omega\Gamma z]
\Delta
j \Psi
\Gamma
j; z]
\Delta
: (4:4)
The function \Psi(j; \Omega\Gamma z]) is the total measure of approximation for the exact functional J [z] by the finite dimen­
sional functional “
J j [“z]. Fig. 2 and Fig. 3 illustrate dependencies of the function \Psi(j; \Omega\Gamma z]) upon the problem
dimension N and e
N , inaccuracy ffi and the value of the stabilizing
functional\Omega\Gamma z]. To express these dependencies
a model problem was considered.
Problem 1: T = 10 5 , fi = 40 ffi , ! = 10 \Gamma4 , v eq = 10 5 , – 0 = 10 \Gamma6 , 4– = 10 \Gamma9 , C 1 = 1, C 2 = 10 \Gamma3 ,
u 1 = 0:5, u 2 = 0:5, ae d = 1, ae l = 10 \Gamma3 , c = 3 \Delta 10 8 . The integral equation (3.1) is nonlinear with respect to
the function z(M ). So it's hard to find the analytical expression of the approximation measure \Psi(j; \Omega\Gamma z]). We
considered several model functions z(M ) and chose maximum values of the function \Psi(j; z]) as values of the
approximation measure under conditions of fixed values of j
and\Omega\Gamma z].
Fig. 2 a) demonstrates a saturation of the approximation measure \Psi(j; \Omega\Gamma z]) at relatively low length N of
the vector “
z (N ú 600). The dependence \Psi(j; \Omega\Gamma z]) on the length e
N of the input data vector represented in
Fig. 2 b) is clear: \Psi(j; \Omega\Gamma z]) ¸ e
N \Gamma2 .

Numerical Methods and Programming, 2002, Vol. 3 5
0 200 400 600 800 1000 1200
43
42
41
40
39
38
a)
0 1000 2000 3000
44
42
40
38
36
b)
Fig. 2. a) The dependence of ln
\Gamma
\Psi(j; \Omega\Gamma z])
\Delta
on N for e
N =
900,\Omega\Gamma z] = 1, ffi = 0. b) The dependence
of ln
\Gamma \Psi(j; \Omega\Gamma z])
\Delta on e
N for N =
900,\Omega\Gamma z] = 1, ffi = 0
Fig. 3 a) illustrates the dependence of ln
\Gamma
\Psi(j; \Omega\Gamma z])
\Delta
on ln ffi . This dependence is characterized by the
saturation at low values of ffi , i.e., there is a horizontal asymptote for \Psi(j; \Omega\Gamma z]) at ffi ! 0. For large ffi the
dependence of the approximation measure \Psi(j; \Omega\Gamma z]) becomes linear with respect to ffi : \Psi(j; z]) ú ffi .
Fig. 3 b) shows the dependence of ln
\Gamma
\Psi(j; z])
\Delta
on ln
\Gamma
\Omega\Gamma z]
\Delta
for ffi = 0. This dependence is linear with
respect to ln
\Gamma
\Omega\Gamma z]
\Delta
. For ffi 6= 0 the form of the function ln
\Gamma
\Psi(j; \Omega\Gamma z])
\Delta
of ln
\Gamma
\Omega\Gamma z]
\Delta
becomes similar to the
dependence of ln
\Gamma
\Psi(j; \Omega\Gamma z])
\Delta
on ln ffi , i.e., there are horizontal and incline asymptotes. The angle of inclination
for the incline asymptote increases with decreasing of ffi .
48 46 44 42 40 38 36 34
42
40
38
36
a)
14 12 10 8 6 4 2 0
54
52
50
48
46
44
42
b)
Fig. 3. a) The dependence of ln
\Gamma \Psi(j; \Omega\Gamma z])
\Delta on ln ffi for N = e
N =
900,\Omega\Gamma z] = 1. b) The dependence
of ln
\Gamma \Psi(j; \Omega\Gamma z])
\Delta on ln
\Gamma
\Omega\Gamma z]
\Delta for N = e
N = 900, ffi = 0
The choice of the regularization parameter ff is crucial for solving of an ill­posed problem. The regularization
parameter should depend on input data, their errors and a method of approximation of the initial problem.
The finite­dimensional generalized principles of discrepancy, quasisolutions and smoothing functional described
in [2] may be used as algorithms to find the regularization parameter ff.
For fixed j define the following variables:
“ – j “ – j = inf
n
“
J j [“z] + \Psi(j; b
\Omega\Gamma “ z]) : “
z 2 b
D
o
; (4:5)
“ ú(ff) = f
\Gamma “ – + \Psi(j; b
\Omega\Gamma “
z ff ])
\Delta
; (4:6)
“
ae(ff) = “
I[“z ff ] \Gamma “ ú(ff); (4:7)
where “
z ff 2 b
Z ff j Arg inf
\Phi c
M ff [“z] : “
z 2 b
D
\Psi
. The choice of the regularization parameter ff according to the finite­
dimensional generalized discrepancy principle consists in finding a solution of the equation with the monotone
function
“
ae(ff) = 0: (4:8)
Then z j ¯
P “
z ! ¯ z at j ! 0, where ¯
P : “
z ! z.
5. The minimization of the smoothing functional. To minimize the smoothing functional c
M ff [“z] we
use the method of projection of conjugate gradients. The iteration process is constructed as follows: the initial
approximation is set as a vector “ z (0) with nonnegative components. At a given point “
z (i) of the minimizing

6 Numerical Methods and Programming, 2002, Vol. 3
sequence the gradient g (i) = grad c
M ff [“z (i) ] = grad “
I[“z (i) ] + ff grad b
\Omega\Gamma “
z (i) ] is calculated. Analytical expressions
for gradients of the functionals “
I [“z], b
\Omega\Gamma “
z] may be obtained from formulas (3.3), (3.4):
@ “
I
@z ij
= \Gamma2h 'h/
N!t X
l=1
H
\Gamma
' i ; / j ; (!t) l
\Delta N–
X
k=1
h
W kl
n
: : :
o @f
@z ij
\Gamma
z ij ; – k + 4–D (' i ; / j ; (!t) l )
\Delta i
; (5:1)
@
b\Omega @z ij
= 2h'h/ z ij + 2 h/
h'\Omega
0
1 + 2 h'
h/\Omega
0
2 ; (5:2)
\Omega 0
1 = 2z ij \Gamma
8 ? !
? :
z ij + z i+1j ; i = 1;
z i\Gamma1j + z i+1j ; i 2 2; n \Gamma 1;
z i\Gamma1j + z ij ; i = n;
\Omega 0
2 = 2z ij \Gamma
8 ? !
? :
z ij+1 + z im ; j = 1;
z ij \Gamma1 + z ij+1 ; j 2 2; m \Gamma 1;
z i1 + z ij \Gamma1 ; j = m:
(5:3)
From these formulas we see that to find grad “
I [“z] it's necessary to calculate @
@z f(z; –):
@
@z f(z; –) = k(–) z 2 k(–)h 0 (z) + h 2 (z)
\Gamma zk(–) + h(z)
\Delta 2 ; @
@z h(z) = C 1 C 2 e \GammaC 2 z : (5:4)
Note that functions f(z; –),
@
@z f(z; –) are undefined at the point z = 0 but they are right­continuous. Thereby
we determine them at this point according to continuity:
f(0; –) = 0; @
@z
f(0; –) = C 1 C 2 k(–)
k(–) + C 1 C 2
: (5:5)
With a knowledge of the gradient g (i\Gamma1) and the direction of descent h (i\Gamma1) at the previous step, it's possible
to find a new direction of descent by appeal to the formula
h (i) = g (i) + fl i h (i\Gamma1) ; fl i =
\Gamma
g (i) \Gamma g (i\Gamma1) ; g (i)
\Delta
\Gamma
g (i\Gamma1) ; g (i\Gamma1)
\Delta (5:6)
(this is the Polak­Ribiere version of the method of conjugate gradients). The next step is to construct the
one­parameter set of the elements
“
z – = P
\Gamma
“ z (i) \Gamma –h (i)
\Delta
; – ? 0; (5:7)
where P is the projector of N­vector onto the nonnegative octant b
D. Then the problem of the one­dimensional
minimization of the functional c
M ff [“z] on the set b
D should be solved. The point, at which the functional c
M ff [“z]
has its minimum on the set b
D, is taken as the next element “ z (i+1) of the minimizing sequence. To solve the
one­dimensional minimization problem it's possible to use the quadratic approximation of the function
oe(–) j c
M ff
\Theta
P (“z (i) \Gamma –h (i) )
\Lambda
\Gamma c
M ff [“z (i) ] (5:8)
at three points 0, – step , 2– step . Then an approximate solution of the one­dimensional minimization problem is
taken in the form:
– min = – step \Delta 1
2
oe 2 \Gamma 4oe 1
oe 2 \Gamma 2oe 1
; oe 1 j oe(– step ); oe 2 j oe(2– step ): (5:9)
One should refine this point if it is necessary. After every N steps the method must be ``revised'' by equating fl i to
zero in formula (5.6). Moreover, the algorithm is to be changed if the inaccuracy in solving the one­dimensional
minimization problem yields
\Gamma h (i) ; g (i)
\Delta
6 0 at the current step.
6. Some features of the numerical implementation. Let's consider some features of the numerical
solution of the problem stated above. According to formula (2.12), Foygt profile k(–) is a convolution of
two functions. The evaluation of the convolution is accompanied by cumbersome calculations, so it would be
erroneous to find it every time. The value of the function k(–) may be obtained by its linear interpolation on
some grid. It's also possible to calculate the function k(–) using the approximation by Matveev's formula
k(–) =
p
ln 2
p
úae v
(1 \Gamma ¸)2 \Gammaj 2
+ 1
úae v
1
1 + j 2
\Gamma ¸(1 \Gamma ¸)
úae v
` 3
2 ln 2 + 1 + ¸
' Ÿ
0:066e \Gamma0:4j 2
\Gamma 1
40 \Gamma 5:5j 2 + j 4
–
; (6:1)

Numerical Methods and Programming, 2002, Vol. 3 7
0 1000 2000 3000 4000 5000
0
1000
2000
3000
a)
0 500 1000 1500 2000 2500
0
100
200
300
b)
Fig. 4. a) The dependence of computing times ratio for grad c
M ff [“z] and c
M ff [“z] on N for e
N = 100.
b) The dependence of computing times ratio for grad c
M ff [“z] and c
M ff [“z] on e
N for N = 400
ae v = 1
2
`
ae l +
q
ae 2
l + 4ae 2
d
'
+ 0:05ae l
/
1 \Gamma 2ae l
ae l +
p
ae 2
l + ae 2
d
!
; j = – \Gamma – 0
ae v
; ¸ = ae l
ae v
: (6:2)
From formulas (3.3) and (5.1) for the discrepancy functional “
I [“z] and its gradient @ “
I
@z ij
one could see that
summing often involves functions defined on one­, two­ or threedimensional grids, for instance, the function
H
\Gamma
' i ; / j ; (!t) l
\Delta
. Such functions should be computed previously, i.e., one can construct a matrix H ijl , which
can be used. Some calculations, e.g. / \Gamma !t to find the Doppler's shift 4–D (' i ; / j ; (!t) l ) or cos ', can be
``put outside brackets'', i.e., an expression, which is common for a certain sum (for summation by k in (3.3),
(5.1)), can be computed before computing the sum and then can be substituted. Obviously, it complicates the
program to solve the problem. But the computing time decreases substantially (by times) because all functions
have a simple structure. However, the summation is performed sometimes over six variables (indexes), as for
@ “
I
@z ij
, that leads to a large computing time. Therefore a correct organization of a summation and a computing
order leads to a substantial reduction of the program runtime.
A minimization of a convex function on a certain set, using first­order methods with computation of function
and gradient values at a point, leads to the question about the required accuracy of solving the one­dimensional
minimization problem. The answer depends on the structure of the function to be minimized. If the time
to compute a function's gradient doesn't exceed the time to compute a function value at a point, then the
one­dimensional minimization problem can be solved with quite low accuracy. It follows, since the program can
make several iterations for an N­dimensional minimization during the time spent on to find the exact minimum
of the one­dimensional function. In this case the value of the function at the found point is often less than
the minimum value of the function on the one­dimensional set considered above. If computing of the gradient
requires a lot of time, then the one­ dimensional minimization should be carried out more precisely. For the
considered functional c
M ff [“z], Fig. 4 shows ratios of computing times for the gradient and for the functional as
functions of lengths N and e
N of the vector “
z and input data respectively. One can see that in the first case the
dependence is linear with respect to N . In the second case the dependence on e
N may be thought of as absent.
Correlations obtained for computing times for the functional and its gradient leads to the conclusion that the
problem of the one­dimensional minimization for the considered multidimensional problem should be solved as
precisely as possible.
The choice of the step – step for the one­dimensional minimization of the convex function oe(–) is an important
problem. If – step is small enough, then values oe 1 and oe 2 may be close and errors of rounding begin to affect
on computation results. Thus, computed values may be not values of a convex function. If the step – step is
too large, the value – min couldn't also be found precisely because of inaccuracies of rounding. To resolve this
issue we can make the step – step dynamic. Let there be a step – step . Find – min using formula (5.9). Then if
– min
– step
2 (0; &] [ [2 \Gamma &; +1), where & is a small number set by the user, then set – step = – min and compute – min .
Shurely, one should take into account that the point “
z corresponding to 2– step should belong to the set of a
priori constraints, i.e., all its coordinates must be nonnegative. Therefore – step should be chosen in such a way
to satisfy the last condition.
Note that when a point “
z (i) of the minimizing sequence is close to the minimum of the functional c
M ff [“z],

8 Numerical Methods and Programming, 2002, Vol. 3
errors in computed functional's values can make a subsequent minimization useless, since the functional's values
at subsequent points “
z (p) , p ? i can exceed c
M ff [“z (i) ]. Thereby, while minimizing the function oe(–), the return
from the subprogram if calculation errors impact on the function's convexity should exist along with standard
conditions of return by the number of iterations or by the gradient norm.
The regularization parameter ff determines a rate of convergence of a minimizing sequence. The functional
b
\Omega\Gamma “
z] is a quadratic function, the functional “
I [“z] hasn't such a property. Therefore, it's obvious that a minimizing
sequence for b
\Omega\Gamma “ z] converges faster than for “
I[“z]. The rate of convergence for “
I [“z] and b
\Omega\Gamma “
z] depending on the
number of iterations p (the number of gradient calculations) is shown in Fig. 5.
0 100 200 300
62
60
58
56
54
52
50
N=64
N=100
N=144
a)
0 100 200 300 400
60
40
20
0
N=900
N=400
N=225
b)
Fig. 5. a) The dependence of ln( “
I[“z]) on the number of minimization iterations. b) The dependence
of ln( b
\Omega\Gamma “ z]) on the number of minimization iterations
To minimize the functional “
I [“z] the exact function is ¯
z = 1 + 0:3 cos ' sin / with parameters from the
Problem 1 and the initial estimate is z j 1. To minimize the functional b
\Omega\Gamma “
z] the exact function is ¯
z j 0, the
initial estimate is z = 1+0:3 cos ' sin /. From Fig. 5 one can see that rates of convergence are close to exponential
ones for considered model problems: Ae \Gammaap and Be \Gammabp (a; b ? 0) for “
I[“z] and “
\Omega\Gamma “ z] respectively. Certainly, the
dependence ln
\Gamma “
I [“z]
\Delta on the number of iterations p showed in Fig. 5 a) isn't linear for sequential minimization
iterations, but a tendency of linear decreasing appears for a great number of minimization iterations, i.e., a right
line \Gammabp + ln B is an asymptote in some sense. From Fig. 5 one can see that the conjugate gradient projection
method quickly jumps into a neighborhood of the minimum. It is demonstrated by characteristic oscillations
near straight lines for small p (the function's values for p ! 20 haven't shown in the figure, a sharp descent
dependent on initial estimate is observed for that p).
It should be noted that there is an observed dependence of constants a, b on the dimension N : a ¸ N
for “
I[“z], b ¸ 1=N for b
\Omega\Gamma “
z]. Thus, during minimization of the functional c
M ff [“z] an exponential decrease of the
functional value with increasing number of iterations will be observed: Ce \Gammacp , c ? 0. The dependence of the
constant c on the dimension N is determined by the regularization parameter ff: c ¸ 1=N for small ff and c ¸ N
for large ff.
7. Parallelization of the minimization problem. The process to solve a minimization problem by
means of the conjugate gradients projection method on the set of nonnegative vectors involves computation
of the functional c
M ff [“z] and its gradient grad c
M ff [“z]. Formulas (3.3), (3.4) for the functional c
M ff [“z] and (5.1),
(5.2) for the gradient @
@z ij
c
M ff lead to the conclusion about a possibility of using multiprocessor systems. The
problem could be parallelized, i.e., one could change the program so that independent parts of the program run
on different processors.
Model problems were solved on the computation cluster in Scientific Research Computer Center of Moscow
State University. This cluster is built on dual Pentium III­processor nodes. The nodes are consolidated into
two­dimensional lattice using high­speed network SCI. Fast Ethernet is used as a service network. See [4] for
more details about the cluster.
The library MPI (Message Passing Interface) [5] was used in programs to solve model chemical elements
distribution mapping problems on the multiprocessor system. The programming language was Fortran 90.
The initial problem is solved using nsize (nsize ? 2) parallel processes with numbers 0; : : : ; nsize \Gamma 1
that are performed on different processors and interact with each other if necessary. A zero process (a process
with number 0) performs all nonparallelized operations: reading from a file, writing into a file, array generation.
The parallelization is implemented only for minimization of the functional c
M ff [“z]: the nonzero processes aren't
involved before and after minimization.
All processes are contained in the group with predefined identifier mpi—comm—world. Each process has
the number irank. Fig. 6 demonstrates the diagram of parallelization for the minimization problem. The

Numerical Methods and Programming, 2002, Vol. 3 9
Sequential part I
istate=1
Evaluation
of the function
Sequential part II
istate=-1
Evaluation
of the gradient
Sequential part III
irank6=0
istate
istate6=0
istate>0
Parallel evaluations
for the gradient
Parallel evaluations
for the function
6=0
6=0
0
Fig. 6. A parallelization scheme for the minimization problem
zero process has three associated blocks, in which a minimization of the functional c
M ff [“z] is performed. If it's
necessary to calculate values of the functional c
M ff [“z], its gradient grad c
M ff [“z] or to return from the subprogram,
the zero process sends a control variable istate to all processes by means of the subprogram mpi—bcast from
MPI:
call—mpi—bcast(istate,1,mpi—integer,0,mpi—comm—world,ierr)
If a value of the variable istate is equal to zero, all processes finish their works and return from the
minimization subprogram. If a value of the variable istate is positive, a calculation of the functional c
M ff [“z] is
performing. If a value of the variable istate is negative, a calculation of the gradient grad c
M ff [“z] is performing.
Each non­zero process calculates values and sends them to the zero process for the subsequent formation of
values of the functional or its gradient by the zero process.

10 Numerical Methods and Programming, 2002, Vol. 3
z
i=1
fun=alpha
omega(: : :)
i e
N
s
fun=fun+s
i=i+1
6=0
any
a)
z
i=irank
i e
N
s=tihfun(i,z,: : :)
s
i=i+nsize-1
0
0
b)
Fig. 7. a) A functional value calculation scheme for the zero process. b) A non­zero process work
scheme for a functional value calculation
Let's consider details of calculating of a value of the functional c
M ff [“z] at a particular point “ z. By reason of
the simple form of the functional b
\Omega\Gamma “
z], its value is calculated by means of one process (zero) using the function
omega(: : : ). Calculation of the functional “
J [“z] may be parallelized by forcing each process to compute the square
of the expression under double­summation sign in (3.3). Thus, the zero process sends a vector “
z (the variable z
in the scheme on Fig. 7 a)) to all other processes. After that every non­zero process calculates the square of the
expression for specified variables k and l. All multidimensional arrays could be expressed as one­dimensional
ones, so in the scheme the indexes k and l from (3.3) are substituted for variable i ranging from 1 to e N. The
square of the expression is computed using the function tihfun(i; z; : : :), each process calculates values for its
indexes i. Surely, it may appear that some processes compute all their functions tihfun(i; z; : : : ) at a certain
point of time while the others compute only the half. Then it's necessary to implement a different multiprocessing
scheme, in which processes involve into the calculations as they become unallocated. For considered modeling
problems calculation times for obtaining values of functions tihfun(i; z; : : : ) were independent on i, i.e., the
suggested scheme appears to be sufficiently fast. Computed values of the function tihfun(i; z; : : :) are sent to
the zero process using the function mpi—send, the zero process performs summation:
call mpi—send(s,1,mpi—double—precision,0,i,mpi—comm—world,ierr)
The zero process receive variables s by means of the function mpi—recv:
call—mpi—recv(s,1,mpi—double—precision,mpi—any—source,mpi—any—tag,mpi—comm—world,istatus,ierr)
As one may notice, the zero process sums variables s in the order of the receipt. It is the property of the
calculation for a value of the functional c
M ff [“z]. The variables s are doubleprecision numbers of the same
order. However, in summation they are rounded. Since the variables s come from different processes in
different order, the rounding leads to differences in last and next to last digits in values of the functional
c
M ff [“z] in repetitive computations. Increasing a number of iterations, one can notice that when the number of
minimization iterations approaches 100 the difference reaches several percents. Moreover, it may appear that
for a first repetitive calculation a minimum has been found (return according to a gradient norm), for a second
one a minimum is been looked for, for a third one a return happens because of the violation of the convexity of
the function to be minimized.
When computing model problems, minimum and maximum values of c
M ff [“z] have been found for each iter­
ation in repetitive computations. Fig. 8 demonstrates dependencies of ratios between function values scattering
width and average value in repetitive calculations on a number of iterations and dimension N . 15 repetitive
calculations have been carried out. These dependencies are characterized by growth with increasing of the

Numerical Methods and Programming, 2002, Vol. 3 11
0 100 200 300
0
0.05
0.1
0.15
0.2
3
6
4
a)
0 100 200 300
0
0.1
0.2
0.3
N=100
N=150
N=400
b)
Fig. 8. The dependence of ratios between function values scattering width and average value in
repetitive calculations on the number of iterations a) for a different number of processes (N = 100)
and b) for a different number of dimension (8 processes)
number of iterations. But by virtue of the fact that in every repetitive computation of minimizing consequences
converges to the one point, sharp ``wells'' may appear in dependences (especially on Fig. 8 b) for N = 400).
It should be pointed out that curves behavior doesn't depend on dimension of a problem or on the number of
processes.
Thus, when solving a minimization problem on multiprocessor computers, a key property of the ordinary
(consecutive) program --- conservation of a result for repetitive computations --- becomes violated. If this
property compliance is obligatory, the solution is simple: the zero process doesn't sum received values but
rather record them in a certain array; once all values are received, a subprogram for summation of the array
elements is called, e.g. using preparatory sorting of elements. In most cases the result conservation in repetitive
computations isn't necessary: a return from the subprogram according to the gradient norm or the convexity
property violation happens only near the minimum point. Thereby, in repetitive computations a minimumpoint
is found anyway, changes apply only to a number of iterations (a time of solving the minimization problem).
A scheme of multiprocessing of the gradient calculation problem grad c
M ff [“z] is shown in Fig. 9. The
parallelization algorithm is the same as for computing the functional c
M ff [“z]. But there is a difference: each
non­zero process calculates k element of a gradient vector grad c
M ff [“z] (the array grad in the scheme in Fig. 9 a))
and sends the calculated element s with its index k in the array to the zero process. Thus, the order of receipt
of elements by the zero process isn't important in formation of the array grad, because indexes k determine
their position in the array grad. It's apparent that the value of the gradient grad c
M ff [“z] that was found using
the suggested scheme is conserved in repetitive computations.
The zero process on schemes presented in the paper is considered as a process that send and receive data
from the other processes. It's possible to use the zero process in basic parallel calculations to calculate values of
the functional c
M ff [“z] and its gradient grad c
M ff [“z]. To reduce computational costs caused by receipt and transfer
messages between processes, and to have the zero process working effectively, i.e., its work wouldn't include
only a message exchange, one can do the following. Each process performs all calculations necessary to find a
value of the functional or its gradient and then send data to the zero process.
In the program for chemical elements distribution mapping problem parallel computations were used only
for calculating values of the Tikhonov's functional c
M ff [“z] and its gradient grad c
M ff [“z]. All rest calculations
performed in series. Thereby using q processors for solving the problem didn't reduce the runtime in q \Gamma 1
times exactly in comparison with the unparalleled program (according to the adopted parallelization scheme,
the zero process doesn't participate in primary parallel calculations, it just send and receive data). It's one of
the expression of the Amdal's empirical law that allows for an unparalleled code in the program. But because
a runtime for a serial code for large numbers N and e
N is negligibly small in comparison with the unparalleled
code, one can assume that the runtime decrease in q \Gamma 1 times on q processors.
Surely, if one consider parallelized blocks of the program, e.g. calculation of the functional and its gradient,
he or she will notice that amounts of calculations for each process should perform differ. Thereby for a moment
of time all processes will finish their work while some processes will just begin last calculations of functions
tihfun(i; z; : : :), tihgrad(k; z; : : :).
For model problems the expression Ü q = t q
t 2
(q \Gamma 1) \Gamma 1 was calculated, where q is the number of processors
and t q is the time of calculation of a gradient on q processors. As would be expected, numbers Ü q approach zero

12 Numerical Methods and Programming, 2002, Vol. 3
z
i=1
grad=0.0
iN
s, k
grad(k)=s
i=i+1
6=0
any
a)
z
k=irank
kN
s=tihgrad(k,z,: : :)
s, k
k=k+nsize-1
0
0
b)
Fig. 9. a) Scheme of the calculation of a gradient for the zero process. b) Scheme of the operation
of a non­zero process during gradient calculation
0 1000 2000 3000 4000 5000
0
1000
2000
3000
a)
0 500 1000 1500 2000 2500
0
100
200
300
400
500
b)
Fig. 10. a) The dependence of a gradient computing time (seconds) grad c
M ff [“z] on N for e
N = 100.
b) The dependence of a gradient computing time (seconds) grad c
M ff [“z] on e
N for N = 400

Numerical Methods and Programming, 2002, Vol. 3 13
with increasing N : Ü 8 = 3:26 for N = 100, Ü 8 = 0:71 for N = 900, Ü 8 = 0:46 for N = 3600.
Fig. 10 illustrates the dependences of a gradient computing time grad c
M ff [“z] on N and e
N of the problem to
be solved for q = 2 (nonparallel program) and allows to evaluate intervals of time required to solve the problem.
The function on Fig. 10 a) is well­approximated by the function 1:38 \Delta 10 \Gamma4 N 2 , the function on Fig. 10 b is
well­approximated by the function 0:221 e
N. Thus, a general dependence of the gradient computing time on N
and e
N is of the form tgrad = 1:387 \Delta 10 \Gamma6 N 2 e
N (seconds). For the computing time of the function c
M ff [“z] we
have t = 2:128 \Delta 10 \Gamma6 N e
N (seconds).
Using multiprocessor systems to solve problems of mapping of the chemical elements distribution on the
stellar surface reduces computation time substantionally. A program with multiprocessing runs for several hours
on computational cluster while an ordinary program (using one processor) runs for several days.
The authors would like to thank Moscow State University Interdisciplinary Scientific Project # 31 ``Mag­
netic Fields and Turbulence in Space''.
###### ##########
1. Goncharskii A.V., Cherepashchuk A.M., Yagola A.G. Ill­Posed Problems in Astrophysics, Nauka, Moscow, 1985 (in
Russian).
2. Tikhonov A.N., Leonov A.S., Yagola A.G. Nonlinear Ill­Posed Problems. Chapman and Hall, 1998.
3. Goncharskii A.V., Stepanov V.V., Khokhlova V.L., Yagola A.G. Restoration of local profiles of lines from observations
in the spectrum of Ap­stars // Letters to Astron. J., 1977, 3, # 6, 278­282 (in Russian).
4. MSU SRCC Computational Cluster (http://parallel.ru/cluster).
5. MPI: A Message­Passing Interface standard // The Message Passing Interface Forum, Version 1.1, June 12, 1995
(http://www.mpi­forum.org).
Received 15 January 2002