Параллельная программа для метода холецкого с помощью openMp и mpi - MPI - Ответ 1441406

/*  Решает систему уравнений для квадратной
 *  симметричной положительно определённой
 *  (с диагональным преобладанием) матрицы
 *  коэффициентов А:
 *        A*x=b;
 *  где A размеров nodesxnodes. Матрица А
 *  предполагается не разреженной. Нумерация 
 *  элементов везде начинается с нуля.
 *  Процедура представляет собой разложение Холесского:
 *        A=L*transpose(L),
 *  после которого выполняются прямое исключение и 
 *  обратная подстановка. A и b не сохраняются. 
*/
void eqsolv_simple_holesskii(Real **A, int nodes, Real *b, Real *x) {
    // Разложение Холесского: замена A верхним и нижним 
    // треугольными множителями.
    A[0][0]=sqrt(A[0][0]);
    A[1][0]/=A[0][0];
    A[0][1]=A[1][0];
    A[1][1]=sqrt(A[1][1]-A[1][0]*A[1][0]);
 
    int irow,irow1;
    int icol, icol1;
    Real sum;
    int k;
    for (irow=2; irow<nodes; irow++) {
        irow1=irow-1;
        A[irow][0]/=A[0][0];
        A[0][irow]=A[irow][0];
        #pragma omp parallel for shared(irow1,A) private(icol, icol1, sum, k)
        for (icol=1; icol<=irow1; icol++) {
            icol1=icol-1;
            sum=0.0;   
            for (k=0; k<=icol1; k++) sum+=A[irow][k]*A[icol][k];
            A[irow][icol]=(A[irow][icol]-sum)/A[icol][icol];
            A[icol][irow]=A[irow][icol];
        }
        sum=0.0;
        #pragma omp parallel for shared(A,irow,irow1) private(k) reduction (+: sum)
        for (k=0; k<=irow1; k++) sum+=A[irow][k]*A[irow][k];
        A[irow][irow]=sqrt(A[irow][irow]-sum);
    }
    
    // Прямое исключение. Происходит разрушение правой части
    b[0]/=A[0][0];
 
    for (irow=1; irow<nodes; irow++) {
        irow1=irow-1;
        sum=0.0;
        #pragma omp parallel for shared(A,b,irow,irow1) private(icol) reduction (+: sum)
        for (icol=0; icol<=irow1; icol++) sum+=A[irow][icol]*b[icol];
        b[irow]=(b[irow]-sum)/A[irow][irow];
    }
 
    // Обратная подстановка используется верхний треугольный множитель
    x[nodes-1]=b[nodes-1]/A[nodes-1][nodes-1];
    for (k=1; k<=nodes; k++) {
        irow=nodes+1-k-1;
        irow1=irow+1;
        sum=0.0;
        #pragma omp parallel for shared(A,x,irow,irow1,nodes) private(icol) reduction (+: sum)
        for (icol=irow1; icol<nodes; icol++) sum+=A[irow][icol]*x[icol];
        x[irow]=(b[irow]-sum)/A[irow][irow];
    }
 
} // eqsolv_simple_holesskii