{
/* Переопределение размера */
integer dim ;
#pragma omp single
{dim = IMVL_dim;}
/* Локальные данные */
integer NumIter = 0;
real NormRight, NormR, Residual,
alpha, alphach, alphazn,
beta, betach, betazn;
/* Необходимые вектора */
real *x = IMVL_CreateVector(IMVL_dim),
*q = IMVL_CreateVector(IMVL_dim),
*r = IMVL_CreateVector(IMVL_dim),
*f = IMVL_right,
*z = IMVL_CreateVector(IMVL_dim);
#pragma omp single
{x = IMVL_CreateVector(IMVL_dim),
q = IMVL_CreateVector(IMVL_dim),
r = IMVL_CreateVector(IMVL_dim),
f = IMVL_right,
z = IMVL_CreateVector(IMVL_dim);}
/* Реализация метода */
IMVL_PreConditioner();
#pragma omp single
{ NormRight = sqrt(IMVL_ScalarMultiply(f,f)); } // |f|
IMVL_SetStartVector(x); // x = x0
IMVL_MultMatrixVector(x,q); // q = Ax
IMVL_SummVectorMultConstant(f,q,-1,r); // r = f-q
IMVL_MultDiagMatrixVector(r,q); // q = D^(-1)r
IMVL_CopyVector(z,q); // z = q
#pragma omp single
{NormR = sqrt(IMVL_ScalarMultiply(r,r));
// |r|
Residual = NormR/NormRight; // |r|/|f|
alphach = IMVL_ScalarMultiply(q,r); } // alch = (q,r)
/* Пока не конец по итерациям или по невязке */
while(Residual>IMVL_EPS && NumIter<MaxIter)
{
IMVL_MultMatrixVector(z,q); // q = Az
#pragma omp single
{alphazn = IMVL_ScalarMultiply(q,z); // alzn = (q,z)
alpha = alphach/alphazn;} // al = alch/alzn
IMVL_SummVectorMultConstant(x,z,alpha,x); // x = x + al * z
IMVL_SummVectorMultConstant(r,q,-alpha,r); // r = r - al * q
IMVL_MultDiagMatrixVector(r,q); // q = D^(-1)r
#pragma omp single
{ betach = IMVL_ScalarMultiply(q,r); // btch = (q,r)
betazn = alphach; // btzn = alch
alphach = betach; // alch = betach
beta = betach/betazn; } // bt = btch/btzn
IMVL_SummVectorMultConstant(q,z,beta,z); // z = q + bt * z
#pragma omp single
{NormR = sqrt(IMVL_ScalarMultiply(r,r)); // |r|
Residual = NormR/NormRight; } // |r|/|f|
NumIter++;
}
/* Указатель ответа настраиваем на решение */
#pragma omp single
{IMVL_Solve = x;}
/* Если выход по невязке */
if(NumIter<MaxIter) return IMVL_RUN_TO_EPSILON;
else return IMVL_RUN_TO_MAX_ITER;
}
/***** Сброс значений матрицы и вектора правой части *****/
void IMVL_Reset()
{
/* Сброс матрицы */
memset(IMVL_val,0,IMVL_row_ptr[IMVL_dim]*sizereal);
/* Сброс вектора правой части */
memset(IMVL_right,0,IMVL_dim*sizereal);
/* Сброс предобуславливающей диагонали */
memset(IMVL_di,0,IMVL_dim*sizereal);
/* Переопределение указателя */
IMVL_start_ptr = IMVL_next_start;
}
Результат работы программы
( время выполнения (сек.))
1. Исходная программа – 0,345
2. Параллельный вариант программы (с OpenMP) – 0.344
Таким образом, время выполнения программы почти не отличается от времени выполнения параллельного варианта, при этом теоретическое ускорение линейно зависит от числа нитей.
Уважаемый посетитель!
Чтобы распечатать файл, скачайте его (в формате Word).
Ссылка на скачивание - внизу страницы.