|
|
|
volodyap |
|
|||
Новичок Профиль Группа: Участник Сообщений: 3 Регистрация: 22.6.2013 Репутация: нет Всего: нет |
1. В документации по системе PETSc не указано: какого типа вещественные данные (single или double precision) используются в процедурах?
2. Более того. Если запустить следующую программу на Фортране из учебного пособия Беликова Д.А. и др. ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ НА МНОГОПРОЦЕССОРНЫХ СИСТЕМАХ program main implicit none #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscksp.h" #include "include/finclude/petscpc.h" Vec x, b, u ! приближенное решение, вектор правых частей, точное решение Mat A ! матрица, определяющая линейную систему KSP ksp ! линейный солвер PC pc ! предобуславливатель Double precision norm ! норма ошибки INTEGER i, j, Ii, Jj, Istart, Iend, m, n, its, ierr, rank Double precision v, one Call PetscInitialize(PETSC_NULL_CHARACTER, ierr) Call MPI_COMM_RANK(PETSC_COMM_WORLD, Rank, Ierr) c Зададим размеры матрицы m = 8 n = 7 c Вычисляем матрицу и вектор правых частей, c которые определяют линейную систему Ax = b c Создаём матрицу, распределенную по активированным процессам Call MatCreate(PETSC_COMM_WORLD, A, ierr) Call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr) Call MatSetFromOptions(A, ierr) c Определяем номера строк, распределенных по локальным процессам Call MatGetOwnershipRange(A, Istart, Iend, ierr) c Заполним матрицу do Ii=istart,iend-1 v = -1.0 i = Ii/n j = Ii - i*n if (i.gt.0) then JJ = Ii - n Call MatSetValues(A, 1, Ii, 1, JJ, v, INSERT_VALUES, ierr) end if if (i.lt.(m-1)) then JJ = Ii + n Call MatSetValues(A, 1, Ii, 1, JJ, v, INSERT_VALUES, ierr) end if if (j.gt.0) then JJ = Ii - 1 Call MatSetValues(A, 1, Ii, 1, JJ, v, INSERT_VALUES, ierr) end if if (j.lt.(n-1)) then JJ = Ii + 1 Call MatSetValues(A, 1, Ii, 1, JJ, v, INSERT_VALUES, ierr) end if v = 4.0 Call MatSetValues(A, 1, Ii, 1, Ii, v, INSERT_VALUES, ierr) end do c Обрабатываем матрицу Call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) Call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) Call MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr) c Создаём параллельные векторы Call VecCreate(PETSC_COMM_WORLD, u, ierr) Call VecSetSizes(u, PETSC_DECIDE, m*n, ierr) Call VecSetFromOptions(u, ierr) Call VecDuplicate(u, b, ierr) Call VecDuplicate(b, x, ierr) c Устанавливаем точное решение (u=1.0), а затем вычисляем c вектор правых частей one=1.0 Call VecSet(u, one, ierr) Call MatMult(A, u, b, ierr) c Cоздадим линейный солвер Call KSPCreate(PETSC_COMM_WORLD, ksp, ierr) c Устанавливаем операторы. Матрица, определяющая линейную систему, c будет также предобуславливающей матрицей Call KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN, ierr) c Метод решения будем устанавливать в командной строке c -ksp_type <type> Call KSPSetFromOptions(ksp, ierr) c Решаем СЛАУ Call KSPSolve(ksp, b, x, ierr) c Проверка решения Call VecAXPY(x, -one, u, ierr) Call VecNorm(x, NORM_2, norm, ierr) Call KSPGetIterationNumber(ksp, its, ierr) if (rank .eq. 0) then if (norm .gt. 1.e-12) then write(6,100) norm, its else write(6,110) its endif endif 100 format('Norm of error ',e10.4,',iterations ',i5) 110 format('Norm of error < 1.e-12,iterations ',i5) Call VecView(u, PETSC_VIEWER_STDOUT_WORLD, ierr) c Уничтожаем все использовавшиееся PETSc объекты Call KSPDestroy(ksp, ierr) Call VecDestroy(u, ierr) Call VecDestroy(x, ierr) Call VecDestroy(b, ierr) Call MatDestroy(A, ierr) c Завершаем работу PETSc Call PetscFinalize(ierr) end то ошибка приближенного решения, посчитанная по среднеквадратичной норме, имеет порядок 1.0e-4. Вот какие результаты приводит автор Запустим программу на 4 процессорах, используя разные методы: mpirun -np 4 laplas -ksp_type gmres Norm of error 0.3097E-04,iterations 10 mpirun -np 4 laplas -ksp_type bcgs Norm of error 0.5074E-04,iterations 6 mpirun -np 4 laplas -ksp_type cg Norm of error 0.3140E-04,iterations 10 Я проверял. Результаты действительно имеют приведенный вид. Но тогда непонятно. Дело в том, что в решаемой задаче точное и приближенное решения равны между собой и равны единице. Поэтому в лучшем случае ошибка должна получаться нулем или на пределе точности: порядка 1.0e-7 для вещественных данных single и порядка 1.0e-13 для вещественных данных double precision. |
|||
|
||||
Фантом |
|
||||
Вы это прекратите! Профиль Группа: Участник Клуба Сообщений: 1516 Регистрация: 23.3.2008 Репутация: 5 Всего: 49 |
По умолчанию - двойной точности.
Откуда сие следует? |
||||
|
|||||
volodyap |
|
|||
Новичок Профиль Группа: Участник Сообщений: 3 Регистрация: 22.6.2013 Репутация: нет Всего: нет |
Точное решение (вектор u) и приближенное решение (вектор x) в данной программе совпадают, так как
перед тем, как решается система уравнений Ax = b, вектор правых частей системы b вычисляется по формуле b = Au. Добавлено через 4 минуты и 3 секунды Вы правы. Я внимательней прочитал документацию и нашел, что действительно вещественные числа используются с двойной точностью. Но тогда программа должна давать порядо ошибки не 1.0e-4 и даже не 1.0e-6, а примерно 1.0e-12 -- 1.0e-13 |
|||
|
||||
Фантом |
|
|||
Вы это прекратите! Профиль Группа: Участник Клуба Сообщений: 1516 Регистрация: 23.3.2008 Репутация: 5 Всего: 49 |
Накопление ошибок при этом все равно возможно, и для плохо обусловленной матрицы можно получить и не такую потерю точности. |
|||
|
||||
volodyap |
|
|||
Новичок Профиль Группа: Участник Сообщений: 3 Регистрация: 22.6.2013 Репутация: нет Всего: нет |
Какое накопление ошибок?
Матрица размерности 56 с пятью ненулевыми диагоналями. Всего ненулевых элементов порядка 250. Да при решении подобной задачи с двойной точностью ошибка должна быть порядка 1.0e-13. Число обусловленности этой матрицы равняется отношению максимального собственного значения к минимальному: lambda_max/lambda_min=7.7271/0.2729=28.3148. Поэтому матрица не является плохо обусловленной. Я думаю дело не в этом, а в том, что для получения решения производится 6 или 10 итераций, что для такой матрицы маловато. Только я не знаю, как проверить с какой точностью решает задачу процедура Call KSPSolve(ksp, b, x, ierr). Обычно в нормальных системах имеется возможность задания необходимой точности epsilon. А в PETSc'e я такой возможности не вижу. Если не разберусь с точностью решения системы линейных алгебраических уравнений в PETSc'e, то не имеет смысла работать с таким пакетом. Точность 1.0e-4 слишком грубая. |
|||
|
||||
0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей) | |
0 Пользователей: | |
« Предыдущая тема | Fortran | Следующая тема » |
|
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности Powered by Invision Power Board(R) 1.3 © 2003 IPS, Inc. |