Форум программистов, компьютерный форум, киберфорум
C/С++ под Linux
Войти
Регистрация
Восстановить пароль
Блоги Сообщество Поиск Заказать работу  
 
Рейтинг 4.92/12: Рейтинг темы: голосов - 12, средняя оценка - 4.92
4 / 4 / 0
Регистрация: 25.03.2011
Сообщений: 28

вторая производная с использованием SSE

12.10.2012, 23:06. Показов 2691. Ответов 10
Метки нет (Все метки)

Студворк — интернет-сервис помощи студентам
Задача: Таблично задана функция, необходимо вычислить вторую производную методом конечных разностей.
Выбранный подход к решению: Считываем в регистры 3 массива по 4 элемента, один со сдвигом влево (копируем с i-1 элемента: строка 65), второй с нулевым сдвигом (строка 66) и третий со сдвигом вправо (копируем с i+1 элемента: строка 67). Делаем вычисления. Копируем обратно. Программа работает, но sse версия работает много медленнее чем обычное решение.
Подозреваю, что проблема в использовании _mm_set_ps для копирования в регистры (строки 65-67) и организации обратного копирования из регистров (строки 69-71).

Вопрос: Как можно оптимизировать/сделать вменяемыми операции копирования?


C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
// 1st_SSE.cpp
//
 
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include <iostream>
using std::cout;
using std::cin;
using std::endl;
 
 
int main()
{
const int N=16*1024+1;
float *d2Ydx2=(float*)malloc(4*sizeof(float));
__m128 Ysse,YRsse,YMsse,YLsse;
__m128 rrs=_mm_set1_ps(-2.0f);
__m128 dX2s=_mm_set1_ps((N*N)/(6.28319f*6.28319f));
__m128 *d2Ydx2SSE=(__m128*) d2Ydx2;
float *Y1=(float*)malloc(N*sizeof(float));
float *Y=(float*)malloc(N*sizeof(float));
float *Y3=(float*)malloc(N*sizeof(float));
__m128 *Y2=(__m128*) Y;
timespec ts_beg, ts_end;
 
float rr,dX2;
const int SSEN=(N-1)/4;
 
dX2=(N*N)/(6.28319f*6.28319f);
rr=-2.0f;
//задаем исходную функцию
for (int i=0;i<=N;i++)
{
Y[i]=sin(i*(6.28319f/N));
Y3[i]=0.f;
Y1[i]=0.f;
}
 
 
 
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
// решение в лоб
for (int i=1;i<=N-1;i++)
{
Y1[i]=(Y[i-1]+rr*Y[i]+Y[i+1])*dX2;
 
};
 
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
float time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
 
 
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_beg);
 
for (int i=1;i<=N-1;i=i+4)
{
Ysse=_mm_set1_ps(0.f);
//копируем в регистры
YLsse=_mm_set_ps(Y[i+2],Y[i+1],Y[i],Y[i-1]);
YMsse=_mm_set_ps(Y[i+3],Y[i+2],Y[i+1],Y[i]);
YRsse=_mm_set_ps(Y[i+4],Y[i+3],Y[i+2],Y[i+1]);
//вычисляем производную
Ysse=_mm_mul_ps(YMsse,rrs);
Ysse=_mm_add_ps(Ysse,YLsse);
Ysse=_mm_add_ps(Ysse,YRsse);
_mm_store_ps(d2Ydx2, _mm_mul_ps(Ysse,dX2s));
//копируем обратно
for (int j=0;j<4;j=j++)
{
Y3[i+j]=d2Ydx2[j];
}
};
 
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts_end);
time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec)/1e9;
cout << "Time required = " << time <<endl;
 
for (int i=0;i<=N;i=i++)
{
if (Y1[i]!=Y3[i]) cout<<i<<"\t"<<Y1[i]<<"\t"<<Y3[i]<<endl;
 
};
free(d2Ydx2);
free(Y1);
free(Y);
free(Y3);
return 0;
}
0
cpp_developer
Эксперт
20123 / 5690 / 1417
Регистрация: 09.04.2010
Сообщений: 22,546
Блог
12.10.2012, 23:06
Ответы с готовыми решениями:

вторая производная
Подскажите, плиз, как найти вторую производную функции в точке?

вторая производная в С
у меня определенная функцияю y=(cos(x)+6)/sinx мне нужно найти вторую производнуюб не пойму как мне это реализовать в С?

Интерполяция и вторая производная
Нужна программа для интерполяции табличных данных и последующего вычисление 2 производной по полученным результатам Данные в архиве

10
1264 / 978 / 384
Регистрация: 02.09.2012
Сообщений: 3,024
13.10.2012, 22:53
фигню написал... удалено
0
4 / 4 / 0
Регистрация: 25.03.2011
Сообщений: 28
14.10.2012, 10:38  [ТС]
Как можно организовать копирование начиная с любого i-го (не кратного 4) элемента массива в регистр и наоборот? Без использования _mm_set_ps и костыля из строк 74-77.
0
1264 / 978 / 384
Регистрация: 02.09.2012
Сообщений: 3,024
14.10.2012, 11:31
попытка номер 2 вот код

Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include <iostream>
using std::cout;
using std::cin;
using std::endl;
 
 
int main(int arg, char* argv[])
{
    const int N = 16 * 1024;
    const float RR = -2.0f;
    const float DX2 = (N * N) / (6.28319f * 6.28319f);
 
    float *YL __attribute__ ((aligned(16))) = (float*)_mm_malloc(N * sizeof(float), 16);
    float *Y __attribute__ ((aligned(16))) = (float*)_mm_malloc(N * sizeof(float), 16);
    float *YR __attribute__ ((aligned(16))) = (float*)_mm_malloc(N * sizeof(float), 16);
 
    float *Y1 __attribute__ ((aligned(16))) = (float*)_mm_malloc(N * sizeof(float), 16);
    float *Y3 __attribute__ ((aligned(16))) = (float*)_mm_malloc(N * sizeof(float), 16);
 
    //задаем исходную функцию
    for (int i = 0; i < N; i++)
    {
        Y[i] = sin(i * (6.28319f / N));
        Y1[i] = 0.f;
        Y3[i] = 0.f;
    }
    for (int i = 0; i < N; i++)
    {
        if (i > 0) YL[i] = Y[i - 1]; else YL[i] = 0.f;
        if (i < (N - 1)) YR[i] = Y[i + 1]; else YR[i] = 0.f;
    }
 
    __m128 RRsse = _mm_set1_ps(RR);
    __m128 DX2sse = _mm_set1_ps(DX2);
 
    timespec ts_beg, ts_end;
 
    // решение в лоб
    clock_gettime(CLOCK_REALTIME, &ts_beg);
    for (int i = 0; i < N; i++)
    {
        Y1[i] = (YL[i] + RR * Y[i] + YR[i]) * DX2;
    }
    clock_gettime(CLOCK_REALTIME, &ts_end);
    float time = (ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec) / 1.0e9;
    cout << "Time required = " << time << endl;
 
    clock_gettime(CLOCK_REALTIME, &ts_beg);
    for (int i = 0; i < N; i = i + 4)
    {
        __m128 y_i = _mm_load_ps(&Y[i]);
        __m128 y_l = _mm_load_ps(&YL[i]);
        __m128 y_r = _mm_load_ps(&YR[i]);
        __m128 result = _mm_mul_ps(y_i, RRsse);
        result = _mm_add_ps(result, y_l);
        result = _mm_add_ps(result, y_r);
        result = _mm_mul_ps(result, DX2sse);
        _mm_store_ps(&Y3[i], result);
    }
    clock_gettime(CLOCK_REALTIME, &ts_end);
    time=(ts_end.tv_sec - ts_beg.tv_sec) + (ts_end.tv_nsec - ts_beg.tv_nsec) / 1.0e9;
    cout << "Time required = " << time << endl;
 
    for (int i = 0; i < N; i++)
    {
        if (Y1[i] != Y3[i])
            cout << i << ";\t" << YL[i] << ";\t" << Y[i] << ";\t" << YR[i] << ";\t" <<
              Y1[i] << ";\t" << Y3[i] << endl;
    }
 
    _mm_free(Y);
    _mm_free(Y1);
    _mm_free(Y3);
 
    return 0;
}


Вот результат (core 2 duo 2,5GHz):
Code
1
2
Time required = 3.7296e-05
Time required = 1.8439e-05
1. Чтобы пользоваться преимуществами быстрой загрузки данных, нужен load, для которого требуется 16-байтовой вырванивание; то есть каждая четверка float-ов в матрице должна лежать по адресу кратному 16.
2. N также придется сделать кратным 16, чтобы массив ровно ложился. У вас было 16*1024+1... +1 - плохо, убрал.
3. Чтобы "выравнено" положить i-1 и i+1 элементы пришлось для них завести отдельные массивы, иначе сдвиг на единицу нарушает выравнивание.

Подозреваю, что оптимизация здесь еще не закончена, будут улучшения, напишите.
1
4 / 4 / 0
Регистрация: 25.03.2011
Сообщений: 28
14.10.2012, 11:48  [ТС]
grgdvo, большое спасибо за развернутый ответ
на самом деле данная задачка просто кусочек основной - трехмерная задача теплопроводности, то есть в ней полный массив данных будет размером Nx*Ny*Nz, причем Nx и Ny (количество элементов сетки вдоль осей x и y) будут не всегда кратные 4 (не переклинит ли load в один момент?). Ну и их количество тоже будет не всегда маленькое, не хотелось бы память забивать еще двумя массивами (для левой и правой точки при вычислении производных), опять же пере присваивание на каждом шаге будет съедать некоторое время.
0
1264 / 978 / 384
Регистрация: 02.09.2012
Сообщений: 3,024
14.10.2012, 15:33
Цитата Сообщение от Mike_Texnik Посмотреть сообщение
не переклинит ли load в один момент?
Переклинит! Если адрес не кратен 16, то будет segmentation fault.

Цитата Сообщение от Mike_Texnik Посмотреть сообщение
Ну и их количество тоже будет не всегда маленькое, не хотелось бы память забивать еще двумя массивами (для левой и правой точки при вычислении производных), опять же пере присваивание на каждом шаге будет съедать некоторое время.).
Это крест параллельных вычислений, который надо нести каждому, кто в это погрузился. Идеального решения не будет! Либо вы подбираете (подгоняете) задачу под конкретную машину и ее конкрентные особенности, либо у вас падает произвиодтельность в разы и больше. Кстати с памятью проблема не очень большая. Сейчас памяти в компьютерах достаточно, она дешевая и можно без надобности ее НЕ экономить.

Если у вас планируется размерности по координатам превышать 100 единиц и выше, то сразу переходите на MPI и другой принцип программирования. С SIMD вы не получите должного ускорения (идеальный потолок для SSE : ускорение в 4 раза).
0
4 / 4 / 0
Регистрация: 25.03.2011
Сообщений: 28
14.10.2012, 16:21  [ТС]
SIMD планирую использовать совместно с MPI (но для начала с OpenMP скрестить хочу), верно ли понимаю, что корректного использования SSE придется отказаться от развертывания 3D массива в 1D ? иначе при умножении на 4 (в начале считывания) будет получаться черти что?
0
 Аватар для Kastaneda
5232 / 3205 / 362
Регистрация: 12.12.2009
Сообщений: 8,143
Записей в блоге: 2
14.10.2012, 21:00
Для чистоты эксперемента советую разделить вычисления разными способами. И запускать тот или иной, например, в зависимости от аргумента командной строки. Т.к. есть вероятность, что при вычислении "в лоб" кэш-мисов больше, чем при дальнейшем вычислении с SSE.
Еще для замера производительности (и времени в т.ч.) могу посоветовать утилиту perf, умеет показывать много всего, в т.ч. и кэш-мисы.
0
1264 / 978 / 384
Регистрация: 02.09.2012
Сообщений: 3,024
15.10.2012, 15:21
Цитата Сообщение от Mike_Texnik Посмотреть сообщение
SIMD планирую использовать совместно с MPI (но для начала с OpenMP скрестить хочу), верно ли понимаю, что корректного использования SSE придется отказаться от развертывания 3D массива в 1D ? иначе при умножении на 4 (в начале считывания) будет получаться черти что?
Что-то вы все в кучу собрали Не сложновато ли?

Мое мнение. От каждой технологии, использованной в программе, должен быть выхлоп. У вас газогидродинамика - сложные уравнения и большие матрицы. Честно говоря SSE не для больших матриц, да и потолок ускорения небольшой. SSE имеет смысл использовать для небольших внутренних вычислений вспомагательных величин, когда на лицо явно "векторные" вычисления (то есть однотипные операции над списком данных).

По поводу матриц: любая, даже 5555-мерная , матрица в памяти компьютера - это одномерный вектор. Все зависит от того как вы двигаетесь по размерностям. Бывают ситуации, когда исходное представление единой трехмерной матрице (в соответствии с физической областью расчета) оказывается невыгодным. Поэтому изначально не зацикливайте себя на этом, выбирайте тот вариант представления, который будет удобен с точки зрения скорости вычислений.
0
4 / 4 / 0
Регистрация: 25.03.2011
Сообщений: 28
26.10.2012, 10:36  [ТС]
Относительно успешно завершились игры с SSE, ниже привожу решение задачки тепловодности для кубика с изолированными стенками. Применил OpenMP и SSE, использование SSE дало ускорение в два раза, что ожидаемо. Однако есть проблема: Visual Studio прекрасно понял OpenMP и на 2 ядрах программка пошла в 2 раза быстрее, а вот gcc нее понял шутку и программа работать даже чуть медленнее.
Опции компиляции: g++ proga.cpp -fopenmp -msse -o out.out

Программка:
Кликните здесь для просмотра всего текста

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
#include <stdio.h>
#include <malloc.h>
#include <xmmintrin.h>
#include <iostream>
#include <time.h>
#include <fstream>
#include <omp.h>
 
using namespace std;
 
const int Nx=128,Ny=128,Nz=128;
int Nth=4; // number of OpenMP threads
 
//Allocation for Visual Studio
//__declspec(align(16)) float U[Nx][Ny][Nz];
//__declspec(align(16)) float U0[Nx][Ny][Nz];
//Allocation for GCC
__attribute__((aligned(16))) float U[Nx][Ny][Nz];
__attribute__((aligned(16))) float U0[Nx][Ny][Nz];
 
 
const float alpha=0.1172f,dt=0.001f; //thermal diffusion coefficient, timestep
const float Lx=4.f,Ly=4.f,Lz=4.0f; //sizes of beam in X, Y and Z directions
const float dX=Lx/float(Nx-1),dY=Ly/float(Ny-1), dZ=Lz/float(Nz-1);
const float dX2=1.0f/(dX*dX),dY2=1.0f/(dY*dY),dZ2=1.0f/(dZ*dZ);
const float E1=1.0f/3.0f, D2=alpha*dt;
float Xi,Yj,Zk;
int k,i,j;
 
    float d2Udx2,d2Udy2,d2Udz2;
    __m128 dX2_sse = _mm_set1_ps(dX2*D2);
    __m128 dY2_sse  = _mm_set1_ps(dY2*D2);
 
    __m128 dZ2_sse  = _mm_set1_ps(dZ2*D2);
    __m128 Gam_sse  = _mm_set1_ps(1.0f-2.0f*D2*(dX2+dY2+dZ2));
 
    __m128 U_Fwd,U_Back,U_Up,U_Down,U_Left,U_Right,U_Cent;
    __m128 U_Res,U_ResTemp;
    
    __m128 E1_sse  = _mm_set1_ps(E1);
    __m128 Chet_sse  = _mm_set1_ps(4.0f);
 
    int Kmin_Weld=int(0.45f*Lz/dZ);
    int Kmax_Weld=int(0.55f*Lz/dZ);
    int Jmin_Weld=int(0.45f*Ly/dY);
    int Jmax_Weld=int(0.55f*Ly/dY);
 
int Init()
{
    for (k=0;k<Nz;k++)
    {
        Zk=dZ*((float)k);
        for (int i=0;i<Nx;i++)
        {
            Xi=dX*((float)i);
            for (int j=0;j<Ny;j++)
            {
                Yj=dY*((float)j);
                U0[i][j][k]=300.0f;
                U[i][j][k]=300.0f;
                if ((i==(Nx-10))&&(Yj>=0.45f*Ly)&&(Yj<=0.55f*Ly)&&(Zk>=0.45f*Lz)&&(Zk<=0.55f*Lz)) U0[i][j][k]=3000.0f;
            }
        }
    }
    return 0;
}
 
float Calc()
{
    
    int jYmin,jYmax,ijk,jMin,jMax;
    #pragma omp parallel private(i,j,k,U_Fwd,U_Back,U_Up,U_Down,U_Left,U_Right,U_Cent,U_Res,U_ResTemp)
    {
    #pragma omp for
    for (i=1;i<Nx-1;i++)
    {
        for (j=1;j<Ny-1;j++)
        {
            for (k=1;k<Nz-1;k=k+4)
            {
                U_Cent = _mm_loadu_ps(&U0[i][j][k]);
                U_Cent = _mm_mul_ps(U_Cent, Gam_sse );
 
                U_Back = _mm_loadu_ps(&U0[i][j][k-1]);
                U_Fwd = _mm_loadu_ps(&U0[i][j][k+1]);
                
                U_ResTemp = _mm_add_ps(U_Back , U_Fwd);
                U_Res = _mm_mul_ps(U_ResTemp, dZ2_sse );
 
                U_Left = _mm_loadu_ps(&U0[i][j-1][k]);
                U_Right = _mm_loadu_ps(&U0[i][j+1][k]);
 
                U_ResTemp = _mm_add_ps(U_Left , U_Right);
                U_ResTemp = _mm_mul_ps(U_ResTemp, dY2_sse );
                U_Res = _mm_add_ps(U_Res, U_ResTemp );
 
                U_Down = _mm_loadu_ps(&U0[i-1][j][k]);
                U_Up = _mm_loadu_ps(&U0[i+1][j][k]);
 
                U_ResTemp = _mm_add_ps(U_Down , U_Up);
                U_ResTemp = _mm_mul_ps(U_ResTemp, dX2_sse );
                U_Res = _mm_add_ps(U_Res, U_ResTemp );
                U_Res = _mm_add_ps(U_Res, U_Cent );
                _mm_storeu_ps(&U[i][j][k], U_Res);
 
                
            }
        }
    }
    }
    
    // BC on XY surfaces and set U to U0
    //begin of the second paralle section
    #pragma omp parallel private(jYmin,jYmax,i,j,k,U_Fwd,U_Back,U_Res,jMin,jMax)
    {
    #pragma omp for
    for (i=0;i<Nx;i++)
    {
        jMin=0;
        jMax=Ny;
        
        jYmin=0;
        jYmax=Ny-1;
        for (j=jMin;j<jMax;j++)
        {
            U[i][j][0]=E1*(4.0f*U[i][j][1]-U[i][j][2]);
            U[i][j][Nz-1]=E1*(4.0f*U[i][j][Nz-2]-U[i][j][Nz-3]);
            for (k=0;k<Nz;k=k+4)
            {
 
                // BC on YZ surfaces
                if (i==0)
                {
                    U_Back = _mm_load_ps(&U[1][j][k]);
                    U_Res = _mm_mul_ps(U_Back , Chet_sse);
                    U_Fwd = _mm_load_ps(&U[2][j][k]);
                    U_Res = _mm_sub_ps(U_Res, U_Fwd );
                    U_Res = _mm_mul_ps(U_Res, E1_sse );
                    _mm_store_ps(&U[0][j][k], U_Res);
                }
 
                // BC on XZ surfaces
                if (i==Nx-1)
                {
                    U_Back = _mm_load_ps(&U[Nx-2][j][k]);
                    U_Res = _mm_mul_ps(U_Back , Chet_sse);
                    U_Fwd = _mm_load_ps(&U[Nx-3][j][k]);
                    U_Res = _mm_sub_ps(U_Res, U_Fwd );
                    U_Res = _mm_mul_ps(U_Res, E1_sse );
                    _mm_store_ps(&U[Nx-1][j][k], U_Res);
                }
 
 
 
 
                if (j==jYmin)
                {
                    U_Back = _mm_load_ps(&U[i][jYmin+1][k]);
                    U_Res = _mm_mul_ps(U_Back , Chet_sse);
                    U_Fwd = _mm_load_ps(&U[i][jYmin+2][k]);
                    U_Res = _mm_sub_ps(U_Res, U_Fwd );
                    U_Res = _mm_mul_ps(U_Res, E1_sse );
                    _mm_store_ps(&U[i][jYmin][k], U_Res);
                }
 
 
                if (j==jYmax)
                {
                    U_Back = _mm_load_ps(&U[i][jYmax-1][k]);
                    U_Res = _mm_mul_ps(U_Back , Chet_sse);
                    U_Fwd = _mm_load_ps(&U[i][jYmax-2][k]);
                    U_Res = _mm_sub_ps(U_Res, U_Fwd );
                    U_Res = _mm_mul_ps(U_Res, E1_sse );
                    _mm_store_ps(&U[i][jYmax][k], U_Res);
                }
                _mm_store_ps(&U0[i][j][k], _mm_load_ps(&U[i][j][k]));
            }
        }
    }
    }
    i=Nx-10;
 
    #pragma omp parallel private(j,k)
    {
    #pragma omp for
    for (int k=Kmin_Weld;k<Kmax_Weld;k++)
    {   
 
        {
            for (j=Jmin_Weld;j<Jmax_Weld;j++)
            {
                U0[i][j][k]=3000.0;
                U[i][j][k]=3000.0;
            }
        }
    }
    }
    return 0;
}
 
 
int main()
{
    ofstream resultTemp;
    ofstream resultNyMnMx;
    ofstream resultU0;
    ofstream Info;
 
    double tstart, tstop, ttime;
    omp_set_num_threads(Nth);
    Init();
    resultU0.width(15);
    resultU0.precision(7);
    resultU0.open("U0.dat");
 
    for (k=0;k<Nz;k++)
    {   
        Zk=dZ*float(k);
 
        i=Nx-10;
        {
            Xi=dX*float(i);
            for (j=0;j<Ny;j++)
            {
                Yj=dY*float(j);
                resultU0<<fixed<<Xi<<"\t"<<Yj<<"\t"<<Zk<<"\t"<< U0[i][j][k] << endl;
            }
        }
    }
 
 
    resultU0.close();
    
 
    cout << "initialization finished" << endl;
    tstart = (double)clock()/CLOCKS_PER_SEC;
    for (int n=1;n<500;n++)
    {
        Calc();
    }
    tstop = (double)clock()/CLOCKS_PER_SEC;
    cout<<"Time required for solution is "<<tstop-tstart<<" Seconds"<<endl;
    resultU0.open("U1.dat");
    for (k=0;k<Nz;k++)
    {   
        Zk=dZ*float(k);
        //for (i=Nx-1;i<Nx;i++)
        i=Nx-10;
        {
            Xi=dX*float(i);
            for (j=0;j<Ny;j++)
            {
                Yj=dY*float(j);
                resultU0<<fixed<<Xi<<"\t"<<Yj<<"\t"<<Zk<<"\t"<< U[i][j][k] << endl;
            }
        }
    }
    resultU0.close();
    cout << "DONE!!!" << endl;
    cin>>Xi;
    return 0;
}



Вопрос: где напакостил, что gcc не хочет работать?
0
1264 / 978 / 384
Регистрация: 02.09.2012
Сообщений: 3,024
26.10.2012, 23:38
Может стоит еще добавить флаг оптимизации?

Code
1
g++ proga.cpp -fopenmp -msse -O3 -o out.out
"-O3" бывает вреден для математики, если увидите, что программа перестала считать верный результат, понизьте уровень до "-O2"
0
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
raxper
Эксперт
30234 / 6612 / 1498
Регистрация: 28.12.2010
Сообщений: 21,154
Блог
26.10.2012, 23:38
Помогаю со студенческими работами здесь

Вторая производная методом полинома Лагранжа
По формулам которые ниже на картинке надо написать функцию, которая на входе получает массив со значениями x, массив со значениями функции...

Вторая производная
Подскажите, пожалуйста, правильно ли понимаю. Пускай f(x1, ..., xn) - функция. Тогда её вторая производная равна : \sum_{i =...

Вторая производная
Проблема такая: У меня есть функция y=f(a,b,c) нужно вычислить вторую производную в точке d: y''=d2(f(a,b,c))/da2|d Как это...

Вторая производная.
Прошу помощи в нахождении производной. Производная вида:...

Вторая производная функции EXP(x)
Все знаки (10 штук) верные REM y''(1) = 2.718281828 DEFDBL A-Z CLS dx = .000005 x = 1 y1 = (EXP(x + 2 * dx) - EXP(x -...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
11
Ответ Создать тему
Новые блоги и статьи
SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий
8Observer8 02.03.2026
Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11680&amp;d=1772460536 Одним из. . .
Реалии
Hrethgir 01.03.2026
Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .
Ритм жизни
kumehtar 27.02.2026
Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .
SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten
8Observer8 27.02.2026
Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .
SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами
8Observer8 20.02.2026
Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .
Конвертировать закладки radiotray-ng в m3u-плейлист
damix 19.02.2026
Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .
Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI
Eddy_Em 18.02.2026
Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .
Камера Toupcam IUA500KMA
Eddy_Em 12.02.2026
Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2026, CyberForum.ru