Талтыкина (1235576), страница 5
Текст из файла (страница 5)
Конец рекурсии определяется максимальным уровнемдерева. Построенное дерево сохраняется в оперативной памяти и являетсявходными данными для следующего этапа метода. Кластер представляет собой структуру данных, в которой описываются какие и сколько точек сеткисодержит кластер, его размеры и положение. Определение этой структуру31приведено на следующем листинге.type cuber e a l ( 8 ) : : x1 ( 3 ) , x2 ( 3 ) , y1 ( 3 ) , y2 ( 3 )r e a l ( 8 ) : : s i d e ! l e n g t h o f cube s i d ereal (8) : : center (3)integer : : geninteger : : nmbr e a l ( 8 ) , a l l o c a t a b l e : : xt ( : , : )type ( cube ) , pointer : : u1 , u2 , u3 , u4 , u5 , u6 , u7 , u8end typeПолями структуры cube являются координаты вершин куба x1, x2, y1, y2;длина его стороны side; центр куба center; уровень дерева gen, на которомнаходится данный куб; число точек сетки nmb, попавших в куб; сами точкисетки xt, а также указатели u1, u2, u3, u4, u5, u6, u7,u8 на дочерние кубы.Каждый блок отвечает взаимодействию пары кластеров, если точки кластеров геометрически удалены друг от друга, то блок попадает в «дальнюю»зону, иначе в «ближнюю».
Расстояние считается по формуле || − ||∞ =max || − ||, где = (1 , 2 , 3 ) и = (1 , 2 , 3 ) – центры кубов(кластеров).По формулам (2.12) определяется, в какую зону попадет блок. Размеры блока определяются количеством точек в кластерах. Список блоков в программеорганизован в виде набора структур типа block.
Описание структуры приведено в следующем листинге.type blockcomplex ( 8 ) , a l l o c a t a b l e : : u ( : , : )complex ( 8 ) , a l l o c a t a b l e : : v ( : , : )integer , a l l o c a t a b l e : : nmb_col ( : )integer , a l l o c a t a b l e : : nmb_row ( : )integer : : nmb1 , nmb2l o g i c a l : : checkendtypetype v e c t o rtype ( block ) , dimension ( : ) , pointer : : data=> null ( )end type v e c t o rЗдесь u и v соответственно строки и столбцы матрицы ; nmb_row и nmb_col32номера точек строк и столбцов; nmb1 и nmb2 количество точек в строках истолбцах соответственно; check - переменная, принимающая значение true, если блок принадлежит «ближней» зоне, и значение false в противном случае;старые и новые номера точек сетки nmb_points и new_nmb_points соответственно.
Старые и новые номера точек сетки нужны на втором этапе метода,так как там происходит локальная перенумерация точек. Предполагается,что точки, попавшие в один блок, нумеруются последовательно, и при переходе к следующему блоку нумерация продолжается. Все построенные блокихранятся в виде массива структур в переменной vector_block.На последнем этапе мозаично-скелетонного метода блоки, попавшие в«дальнюю» зону, аппроксимируются суммами одноранговых матриц. Согласно алгоритму 1 ищутся строки u и столбцы v и рассчитывается погрешностьпо формуле (2.16). Критерий остановки определяется формулой (2.14). Такжесуммируется мозаичный ранг по формуле (2.11).ep s = 1step = 1s 2=0d0nmb=min(nmb1 , nmb2)a l l o c a t e ( u1 (nmb1 ) ,u_(nmb1 , nmb) ,v_(nmb2 , nmb ) )i j (2) = 1do while ( eps>epsilon_app ) ! ! ! i n c o m p l e t e c r o s s ap p ro xi m at io ntmp_vector = 0d0do l = 1 , s t ep −1c a l l zaxpy (nmb2 , u_( i j ( 2 ) , l ) ,v_ ( : , l ) , 1 , tmp_vector ( 1 : nmb2 ) , 1 )enddoi f ( nomer == 1 ) thendo l l = 1 , nmb2c a l l yadro1_element ( np , dinteg_mod ( tmp%nmb_row( l l ) ) ,d i n t e g ( tmp%nmb_col ( i j ( 2 ) ) ) , xt , ksr , tmp%nmb_row( l l ) ,tmp%nmb_col ( i j ( 2 ) ) , v_( l l , s t e p ) )enddoelsedo l l = 1 , nmb2c a l l yadro2_element ( np , d i n t e g ( tmp%nmb_col ( i j ( 2 ) ) ), xt , cosnor , ksr , tmp%nmb_row( l l ) , tmp%nmb_col ( i j ( 2 ) ) ,vn , v_( l l , s t e p ) )enddo33endifv_ ( : , s t e p ) = v_ ( : , s t e p ) − tmp_vector ( 1 : nmb2)i j ( 1 ) = maxloc ( cdabs (v_ ( : , s t e p ) ) , 1 )e l = v_( i j ( 1 ) , s t e p )tmp_vector = 0d0do l = 1 , s t ep −1c a l l zaxpy (nmb1 , v_( i j ( 1 ) , l ) ,u_ ( : , l ) , 1 , tmp_vector ( 1 : nmb1 ) , 1 )enddoi f ( nomer == 1 ) thendo l l = 1 , nmb1c a l l yadro1_element ( np , d i n t e g ( tmp%nmb_row( i j ( 1 ) ) ) ,d i n t e g ( tmp%nmb_col ( l l ) ) , xt , ksr , tmp%nmb_row( i j ( 1 ) ) ,tmp%nmb_col ( l l ) ,u_( l l , s t e p ) )enddoelsedo l l = 1 , nmb1c a l l yadro2_element ( np , d i n t e g ( tmp%nmb_col ( l l ) ), xt , cosnor , ksr , tmp%nmb_row( i j ( 1 ) ) , tmp%nmb_col ( l l ) ,vn , u_( l l , s t e p ) )enddoendifu_ ( : , s t e p ) = u_ ( : , s t e p ) − tmp_vector ( 1 : nmb1)u1 ( : ) = cdabs (u_ ( : , s t e p ) )u1 ( i j ( 2 ) ) = 0d0i j ( 2 ) = maxloc ( u1 , 1 )u_ ( : , s t e p ) = u_ ( : , s t e p ) / e lcom = com + d s q r t ( eps1 )eps1 = dot_product (u_ ( : , s t e p ) ,u_ ( : , s t e p ) ) *dot_product (v_ ( : , s t e p ) ,v_ ( : , s t e p ) )s 2 = s 2 + 2d0* d s q r t ( eps1 )* com+eps1ep s = (nmb−s t e p )* d s q r t ( eps1 / s2 ) !step = step + 1enddomosaic_rank = mosaic_rank + min ( ( s te p −1)*(nmb1+nmb2 ) , nmb1*nmb2)Здесь eps – погрешность аппроксимации блока на текущем шаге step;34epsilon_app задает относительную погрешность аппроксимации, одинаковуюдля всех блоков «дальней» зоны; eps1, s2 и com – вспомогательные переменные; ij – массив из двух чисел, содержащий номера максимальных элементов строки и столбца на текущем шаге; u и v – строки и столбцы матрицы;tmp_vector и u_1 – вспомогательные массивы; zaxpy - процедура библиотеки Intel MKL, предназначенная для умножения константы на вектор; nomer– переменная, которая определяет, какое уравнение (1.15) или (1.17) используется; dinteg – массив, в котором хранятся численные значения интегралов (2.2); tmp - экземпляр структуры block; ksr - волновое число; процедураyadro1_element рассчитывает элемент матрицы, аппроксимирующей оператор () в формулах (1.13), а процедура yadro2_element в зависимости отпараметра vn вычисляет элемент матрицы, приближающей операторы ()*или ()в формулах (1.13); cosnor – массив направляющих косинусов; el –элемент матрицы, стоящий на пересечении строки u и столбца v; функцияcdabs возвращает модуль комплексного числа; maxloc – индекс максимального элемента массива; dot_product возвращает результат скалярного произведения векторов; dsqrt рассчитывает корень числа, заданного в аргументефункции; mosaic_rank – переменная, накапливающая мозаичный ранг всейматрицы.Блоки «ближней» зоны рассчитываются для задач 1 и 2 по формулам(2.4), а для задачи 3 по формуле (1.13).
Такие блоки хранятся в памяти ввиде матриц. Их ранг, равный произведению количества строк и количествастолбцов, также добавляется к общему мозаичному рангу. Алгоритм расчетаблоков «ближней» зоны приведен на следующем листинге.a l l o c a t e ( matr (nmb2 , nmb1 ) )i f ( nomer == 1 ) thendo i i = 1 ,nmb2do j j = 1 ,nmb1c a l l yadro1_element ( np , dinteg_mod ( tmp%nmb_row( i i ) ) ,dinteg_mod ( tmp%nmb_col ( j j ) ) , xt , ksr , tmp%nmb_row( i i ) ,tmp%nmb_col ( j j ) , matr ( i i , j j ) )enddoenddoelsedo i i = 1 ,nmb235do j j = 1 ,nmb1i f ( tmp%nmb_col ( j j ) .
ne . tmp%nmb_row( i i ) ) thenc a l l yadro2_element ( np , dinteg_mod ( tmp%nmb_col ( j j ) ) , xt ,cosnor , ksr , tmp%nmb_row( i i ) , tmp%nmb_col ( j j ) , vn , matr ( i i , j j ) )elsei f ( nomer==2) matr ( i i , j j ) = d i a g ( tmp%nmb_row( i i ) )i f ( nomer==3) matr ( i i , j j ) = diag_norm ( tmp%nmb_row( i i ) )endifenddoenddoendifdeallocate ( matr )mosaic_rank = mosaic_rank + nmb1*nmb2В данном листинге за matr обозначена матрица, в которой хранится блок«ближней» зоны.В GMRES умножение блока на вектор зависит от зоны, в которую попадает блок. Блоки «ближней» зоны используют функцию из библиотекиIntel MKL, которая требует 2 операций для умножения квадратной матрицы × на вектор размера .
Для умножения блоков «дальней» зонынаписана специальная функция block_vector. Она позволяет умножить блокразмера × на вектор за линейное число операций. Ниже, в листинге,приведено определение этой функции.subroutine b l o c k _ v e c t o r ( k , m1, m2, u , v , vec , r e s )include "mkl.fi"integer : : k , i , m1, m2complex ( 8 ) : : u (m1, k ) , v (m2, k )complex ( 8 ) : : vec (m1) , r e s (m2)complex ( 8 ) : : comr e s = ( 0 d0 , 0 d0 )do i = 1 , kcom = dot_product ( u , vec )c a l l zaxpy (m2, com , v ( : , i ) , 1 , r e s , 1 )enddoend subroutineЗдесь res – результат умножения скелетона * на вектор vec, m1 и m2 раз36мерности соответственно векторов vec и res. Результатом работы программырешения краевых задач и задач дифракции является файл с плотностямисоответствующих интегральных уравнений. Решения исходных задач восстанавливаются в другой программе расчёта погрешностей решений.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
