Поиск:

Ответ в темуСоздание новой темы Создание опроса
> [PETSc] Типы вещественных данных, single или double precision 
:(
    Опции темы
volodyap
Дата 23.6.2013, 08:52 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 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.

PM MAIL   Вверх
Фантом
Дата 23.6.2013, 11:17 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Вы это прекратите!
***


Профиль
Группа: Участник Клуба
Сообщений: 1516
Регистрация: 23.3.2008

Репутация: 5
Всего: 49



Цитата(volodyap @  23.6.2013,  09:52 Найти цитируемый пост)
1. В документации по системе PETSc не указано: какого типа вещественные данные (single или double precision) используются в процедурах?

По умолчанию - двойной точности.

Цитата(volodyap @  23.6.2013,  09:52 Найти цитируемый пост)

Но тогда непонятно. Дело в том, что в решаемой задаче точное и приближенное решения равны между собой и равны единице.

Откуда сие следует?
PM   Вверх
volodyap
Дата 23.6.2013, 12:09 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 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
PM MAIL   Вверх
Фантом
Дата 23.6.2013, 13:06 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Вы это прекратите!
***


Профиль
Группа: Участник Клуба
Сообщений: 1516
Регистрация: 23.3.2008

Репутация: 5
Всего: 49



Цитата(volodyap @  23.6.2013,  13:09 Найти цитируемый пост)
Точное решение (вектор u) и приближенное решение (вектор x) в данной программе совпадают, так как
перед тем, как решается система уравнений Ax = b, вектор правых частей системы b вычисляется по формуле b = Au.

Накопление ошибок при этом все равно возможно, и для плохо обусловленной матрицы можно получить и не такую потерю точности.
PM   Вверх
volodyap
Дата 23.6.2013, 20:14 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 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 слишком грубая.

PM MAIL   Вверх
  
Ответ в темуСоздание новой темы Создание опроса
0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Fortran | Следующая тема »


 




[ Время генерации скрипта: 0.1194 ]   [ Использовано запросов: 22 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.