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

#pragma omp for private(i,IMVL_di) shared(IMVL_EPS) schedule(static,2)

        /* Преобразование диагонали */

            { for(i=0 ; i<IMVL_dim ; i++)

        {

            if(fabs(IMVL_di[i])<IMVL_EPS)

                return IMVL_NULL_DIAGONAL;

            else

                IMVL_di[i] = 1.0/IMVL_di[i];

        }

            }

        return IMVL_OK;

}

/***** Уничтожение матрицы и векторов *****/

void IMVL_Destroy()

{

        #ifdef __cplusplus

           delete IMVL_memory;

        #else

           free(IMVL_memory);

        #endif

}

/***** Умножение матрицы на вектор *****/

void IMVL_MultMatrixVector( real* InVector,  // входящий вектор

                            real* OutVector  // полученный вектор

                          )

{

        /* Переопределение указателей и размера матрицы */

      integer *col_ind ,*row_ptr ;

        integer  dim ;

            real    *val;

       #pragma omp single    

      {     col_ind = IMVL_col_ind,

        row_ptr = IMVL_row_ptr;

        dim = IMVL_dim;

            val = IMVL_val;}

        /* Локальные данные */

        integer i, k, end;

        real    valelem;

        /* Умножение матрицы на вектор */

#pragma omp for private(i,k,end,valelem) shared(row_ptr,col_ind,InVector,OutVector,val) schedule(static,2)       

            {for(i=0 ; i<dim ; i++)

        {

            valelem = 0.0;

            end = row_ptr[i+1];

            /* умножение только не нулевых элементов */

            for(k=row_ptr[i] ; k<end ; k++)

            {

               valelem += val[k] * InVector[col_ind[k]];

            }

            OutVector[i] = valelem;

            }}

}

/***** Умножение вектора на обратную диагональ матрицы *****/

void IMVL_MultDiagMatrixVector( real* InVector,  // входящий вектор

                                real* OutVector  // полученный вектор

                              )

{

        /* Переопределение размера и указателей матрицы */

      integer dim ;

      real    *di ;

#pragma omp single

      {dim = IMVL_dim;

      di = IMVL_di;}

        /* Локальные перемнные */

        integer i;

        /* Умножение диагонали */

#pragma omp for private(i) shared(OutVector,InVector, di) schedule(static,2)    

    {for(i=0 ; i<dim ; i++)

                  OutVector[i] = di[i] * InVector[i];}

}

/***** Сложение векторов с домножением X=Y+aZ *****/

void IMVL_SummVectorMultConstant( real* y, real* z,

                                  real al, real* x)

{

        /* Переопределение размера матрицы */

integer dim ;   

#pragma omp single     

{dim = IMVL_dim;}

        /* Локальные данные */

        integer i;

        /* Вычислительный цикл */

#pragma omp for  shared(x,y,z,al)private(i) schedule(static,2)

            {for(i=0 ; i<dim ; i++)

                  x[i] = y[i] + al * z[i];}

}

/***** Скалярное произведение *****/

real IMVL_ScalarMultiply( real* x, real* y)

{

        /* Переопределение размера матрицы */

integer dim ;   

#pragma omp single     

{dim = IMVL_dim;}

        /* Локальные данные */

        integer i;

        real dot = 0.0;

        /* Вычислительный цикл */

#pragma omp for private(i) shared(x,y) schedule(static,2) reduction(+:dot)

            { for(i=0 ; i<dim ; i++) dot += x[i] * y[i];}

        /* Скалярное проивздение */

        return dot;

}

/***** Установка начального приближения *****/

void IMVL_SetStartVector(real *x)

{

       /* установка нулевого начального приближения */

    #pragma omp single  

      {memset(x,0,IMVL_dim*sizereal); }

}

/***** Копирование векторов *****/

void IMVL_CopyVector(real *dest, real *src)

{

    #pragma omp single   

      {memcpy(dest,src,IMVL_dim*sizereal);}

}

/***** Метод сопряжённых градиентов *****/

integer IMVL_ConjGradientMethod(integer MaxIter)