. . 6 , , . .
. 6. π .
1 , . ; zero (. . 4) ; , . 2 normsq : .
3 . , [7]. , . v -
. ( , ).
, . (. . 7) . , .
. 7. in,
x k 𝛼.
. 7 . 11 (. . 5.). , .
.
.
( ). , . , . , .
: , , .
++ OpenMP, ( 2). . . , . .
|
|
1). .. // : . . . (30 - 2 2000 ., . ). .: - . 2000. C. 71-72.
2). .., .., .. // . No. 2. 2004. . 16-32.
3). .., .. . . : , 2005. 210 .
4). .., .. . : - , 1994. 103 .
5). .., .. . .: -, 2002. 600 .
6). .., .., .., .. // . 2009. . 49. 8. . 1369-1384.
7). .., .. // . " ". 2011. No. 37(254), . 10. C. 12-21.
8). A.., . ., . . : . .: , 1985. 161 .
. , . .: , 1966. 600 .
9). . . : "", 2001. 624 .
10). .. . : - , 2009. 200 .
1
( 𝛼); ; A | |
( ) | |
( 𝛼); b | |
, . | |
c | |
C | |
dataChange(t) | t (t 1) |
σ | xk |
ε | ( : ) |
g | |
G | g |
γ | |
init | |
factor | |
k | |
kγ | |
k𝛼 | 𝛼, ; k𝛼 = MinInt , |
K | |
L | |
m | |
MinFloat | , |
MinInt | , |
n | |
normsq | : |
p | , |
P | MPI- |
π(z, k) | z , k; xk (-1,), |
q | |
q' | |
r | |
R | r, |
rank() | , MPI- |
s | : |
T | z |
w | |
k- ; , k- | |
y | 𝛼 |
z | |
zero(k) | -, k |
|
|
2
#include "include.h"
using namespace std;
int main(int argc, char *argv[]) {
LoadData(); // input.txt
Init(); //
Pi(/*axis*/ 0,/*index*/ 0); // _x _z index axis
cout << "x="; for (int j = 0; j<_n; j++) cout << _x[j] << "\t"; cout << endl; fflush(stdout);
printf("\nPress ENTER"); fflush(stdout); getchar();
return 0;
};
static void Init() {//
_mc = _m + _n * 2; //
assert(_MAX_MC >= _mc);
#ifndef _NO_MPI
MPI_Init(&argc, &argv); // MPI
MPI_Comm_rank(MPI_COMM_WORLD, &_rank); // _rank MPI
MPI_Comm_size(MPI_COMM_WORLD, &_size); // _size MPI
#else
_rank = 0;
#endif
//
for (int i = 0; i<_n; i++)
_A[i + _m][i] = -1;
for (int i = 0; i<_n; i++)
_A[i + _m + _n][i] = 1;
_r = _R; //
_s = _r / ((_K-1)/2); //
for (int i = 0; i < _n; i++) //
_z[i] = _c[i] * _T;
};
static void LoadData() {// input.txt
FILE *stream;
fpos_t pos;
errno_t errNo = fopen_s(&stream, _INPUT, "r");
if (errNo!= 0) {
cout << "File '" << _INPUT << "' not found!" << endl;
#ifndef _NO_MPI
MPI_Finalize();
#endif
printf("\nPress ENTER"); fflush(stdout); getchar();
abort();
};
SkipComments(stream, &pos);
fscanf_s(stream, "%d", &_n); //
assert(_n <= _MAX_N);
|
|
assert(_n>0);
SkipComments(stream, &pos);
fscanf_s(stream, "%d", &_m);// - ( -x_j<=0)
assert(_m <= _MAX_M);
assert(_m>0);
SkipComments(stream, &pos);
for (int i = 0; i<_m; i++) // ( -x_j<=0)
for (int j = 0; j<_n; j++)
fscanf_s(stream, "%lf", &_A[i][j]);
SkipComments(stream, &pos);
for (int i = 0; i<_m; i++) // b ( -x_j<=0)
fscanf_s(stream, "%lf", &_b[i]);
SkipComments(stream, &pos);
for (int j = 0; j<_n; j++) // c
fscanf_s(stream, "%lf", &_c[j]);
SkipComments(stream, &pos);
for (int j = 0; j<_n; j++) //
fscanf_s(stream, "%lf", &_g[j]);
SkipComments(stream, &pos);
fscanf_s(stream, "%lf", &_R); // C ,
SkipComments(stream, &pos);
fscanf_s(stream, "%d", &_K); // C
assert(_K>0);
assert(_K % 2 == 1);
SkipComments(stream, &pos);
fscanf_s(stream, "%d", &_L); // C
SkipComments(stream, &pos);
fscanf_s(stream, "%lf", &_T); // C _z
fclose(stream);
};
static void SkipComments(FILE *stream, fpos_t *pos) {
/* (, '%') */
fscanf_s(stream, "\n");
fgetpos(stream, pos);
while (getc(stream) == '%') {
while (getc(stream)!= 10);
fscanf_s(stream, "\n");
fgetpos(stream, pos);
};
fsetpos(stream, pos);
};
static void Pi(int axis, int index) { // _x _z index axis
#define b_alpha(i) _b[i+_m] // b'
Point x_stroke; //
double S[_MAX_N]; //
double _Scal_ai_x[_MAX_MC]; // A' Scal _x
double factor; //
double dist; //
//#pragma omp parallel num_threads(4)
//#pragma omp parallel for simd
for (int j = 0; j < _n; j++)
_x[j] = _z[j];
Zero(axis, index);
for (int i = 0; i < _n; i++)
b_alpha(i) = -_y[i];
for (int i = 0; i < _n; i++)
b_alpha(i + _n) = _y[i] + _s;
CalcNormsq();
//
do {
//#pragma simd
for (int l = 0; l < _L; l++) {
for (int j = 0; j < _n; j++) {
x_stroke[j] = _x[j];
S[j] = 0;
};
for (int i = 0; i<_mc; i++) {
_Scal_ai_x[i] = 0;
for (int j = 0; j<_n; j++) // A' Scal _x
_Scal_ai_x[i] += _A[i][j] * _x[j];
factor= __max(_Scal_ai_x[i] - _b[i], 0) / _normsq[i];
for (int j = 0; j < _n; j++)
S[j] += factor*_A[i][j];
};
for (int j = 0; j < _n; j++)
_x[j] = x_stroke[j] - (_lambda / _m)*S[j];
};
DataChange(); //
dist = Dist(_x, x_stroke);
}while (dist>_EPS_F);
};
static void CalcNormsq() { // a_i
for (int i = 0; i<_mc; i++) {
double s = 0;
for (int j = 0; j<_n; j++)
|
|
s += _A[i][j] * _A[i][j];
_normsq[i] = s;
};
};
static void Zero(int axis, int index) { // y[*] index axis
assert(axis >= 0 && axis<_n);
assert(abs(index) < (_K - 1) / 2);
for (int j = 0; j < _n; j++)
_y[j] = _g[j];
_y[axis] = _g[axis] + index*_s;
};
static void DataChange() {};
static double Dist(Point x1, Point x2) { //
double x, dist;
dist = 0;
for (int j = 0; j<_n; j++) {
x = x1[j] - x2[j];
dist += x*x;
};
dist = sqrt(dist);
return dist;
};
[1]) .