Последовательный и параллельный варианты программы. Теоретическое и практическое ускорение параллельной программы, страница 2

        B[1][0] = u1; B[1][1] = ud; B[1][2] = u3; B[1][3] = u2;

        B[2][0] = u2; B[2][1] = u3; B[2][2] = ud; B[2][3] = u1;

        B[3][0] = u3; B[3][1] = u2; B[3][2] = u1; B[3][3] = ud;

        /* Формирование матрицы масс */

        ud = hx*hy/9.0;

        u1 = hx*hy/18.0;

        u2 = u1;

        u3 = hx*hy/36.0;

        C[0][0] = ud; C[0][1] = u1; C[0][2] = u2; C[0][3] = u3;

        C[1][0] = u1; C[1][1] = ud; C[1][2] = u3; C[1][3] = u2;

        C[2][0] = u2; C[2][1] = u3; C[2][2] = ud; C[2][3] = u1;

        C[3][0] = u3; C[3][1] = u2; C[3][2] = u1; C[3][3] = ud;

        /* Запуск элементной постановки */

        for(i=0 ; i<nx ; i++)

          for(k=0 ; k<ny ; k++)

            {

               /* Нахождение коэффициента теплопроводности */

               mu = GetMu(i*hx+0.5*hx, k*hy+0.5*hy);

               /* Формирование матрицы правой части */

               f1 = GetF(i*hx, k*hy);

               f2 = GetF((i+1)*hx, k*hy);

               f3 = GetF(i*hx, (k+1)*hy);

               f4 = GetF((i+1)*hx, (k+1)*hy);

               F[0] = C[0][0]*f1+C[0][1]*f2+C[0][2]*f3+C[0][3]*f4;

               F[1] = C[1][0]*f1+C[1][1]*f2+C[1][2]*f3+C[1][3]*f4;

               F[2] = C[2][0]*f1+C[2][1]*f2+C[2][2]*f3+C[2][3]*f4;

               F[3] = C[3][0]*f1+C[3][1]*f2+C[3][2]*f3+C[2][3]*f4;

               /* Нахождение опорной точки */

               n1 = k*(nx+1)+i;

               n2 = k*(nx+1)+i+1;

               n3 = (k+1)*(nx+1)+i;

               n4 = (k+1)*(nx+1)+i+1;

               /* Добавление элемента в вектор правой части */

               IMVL_AddElementVector(n1,F[0]);

               IMVL_AddElementVector(n2,F[1]);

               IMVL_AddElementVector(n3,F[2]);

               IMVL_AddElementVector(n4,F[3]);

               /* Добавление элементов в матрицу */

               IMVL_AddElement(n1,n1,mu*B[0][0]);

               IMVL_AddElement(n1,n2,mu*B[0][1]);

               IMVL_AddElement(n1,n3,mu*B[0][2]);

               IMVL_AddElement(n1,n4,mu*B[0][3]);

               IMVL_AddElement(n2,n1,mu*B[1][0]);

               IMVL_AddElement(n2,n2,mu*B[1][1]);

               IMVL_AddElement(n2,n3,mu*B[1][2]);

               IMVL_AddElement(n2,n4,mu*B[1][3]);

               IMVL_AddElement(n3,n1,mu*B[2][0]);

               IMVL_AddElement(n3,n2,mu*B[2][1]);

               IMVL_AddElement(n3,n3,mu*B[2][2]);

               IMVL_AddElement(n3,n4,mu*B[2][3]);

               IMVL_AddElement(n4,n1,mu*B[3][0]);

               IMVL_AddElement(n4,n2,mu*B[3][1]);

               IMVL_AddElement(n4,n3,mu*B[3][2]);

               IMVL_AddElement(n4,n4,mu*B[3][3]);

            }

        /* Учёт краевых условий */

        for(k=0;k<=ny;k++)

        {

            IMVL_AddElement(k*(nx+1),k*(nx+1),1.0e+20);

            IMVL_AddElementVector(k*(nx+1),1.0e+20);

            IMVL_AddElement(k*(nx+1)+nx,k*(nx+1)+nx,1.0e+20);

            IMVL_AddElementVector(k*(nx+1)+nx,0.0);

        }

        /* Решение СЛАУ */

            printf("Max Num threads in region is = %d\n",omp_get_max_threads());

        #pragma omp parallel

                {

                   #pragma omp barrier

                        {

                       IMVL_ConjGradientMethod(5*nx*ny);

                         #pragma omp single

                 printf("Num threads in region is = %d\n",omp_get_num_threads());

                    }

                }

        /* Сохранение решения */

        FSolve = fopen("solve.dat","w");

        for(k=0;k<ny;k++)

          for(i=0;i<=nx;i++)

            fprintf(FSolve,"%lf\t%lf\t%lf\n",

                    i*hx,k*hy,IMVL_Solve[k*(nx+1)+i]);

        fclose(FSolve);

            /*for(i=0;i<100000000;i++)

            {