Macierze

...
Mat A;
...

należy utworzyć, podobnie jak w przypadku wektorów, obiekt typu macierz. Robi się to przy pomocy funkcji: MatCreate(MPI_Comm comm,int m,int n,Mat *A).Funkcja tworzy macierz A o rozmiarze , gdzie m to liczba wierszy, a n kolumn. Wstawianie wartości do macierzy odbywa się analogicznie jak w przypadku wektorów, przy pomocy funkcji MatSetValue(Mat A,int row,int col,Scalar val,InsertMode mode). Ta funkcja wpisuje w macierz A wartość val w miejsce row,col. Ostatni argument określa, podobnie jak w przypadku wektorów, sposób wstawiania wartości do macierzy i może przyjąć wartości INSERT_VALUES lub ADD_VALUES. Do zwolnienia pamięci zajmowanej przez macierz służy funkcja MatDestroy(Mat A), która działa tak samo jak jej odpowiednik dla wektorów.

Przykład

#include "sles.h"

static char help[]="przykład_3\n";
int main(int argc,char **args)
{
  Mat     A; 
  int i,n=5;
  /* inicjalizacja */
  PetscInitialize(&argc,&args,PETSC_NULL,help);
  /* wygenerowanie macierzy trójdiagonalnej*/
  MatCreate(MPI_COMM_WORLD,n,n,&A);

  /* wypełniamy macierz wierszami */
  for(i=0;i<n;i++)
  {
   MatSetValue(A,i,i,2,INSERT_VALUES);
   if(i<n-1)MatSetValue(A,i,i+1,1,INSERT_VALUES);
   if(i>0)MatSetValue(A,i,i-1,1,INSERT_VALUES);
  }
  MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
  
  /* wypisuje macierz na konsolę */
  MatView(A,VIEWER_STDOUT_SELF);
  
  /* zakończenie */
  MatDestroy(A);
  PetscFinalize();
  return(0);
 }

hydra:/usr/home/bin/> przykład
row 0: 0 2  1 -1 
row 1: 0 -1  1 2  2 -1 
row 2: 1 -1  2 2  3 -1 
row 3: 2 -1  3 2  4 -1 
row 4: 3 1  4 2 

hydra:/usr/home/bin/> _

Przyjrzyjmy się, jak działa przeglądarka macierzy MatView() w trybie tekstowym: macierz jest wypisywana wierszami. W każdym wierszu podawane są tylko niezerowe elementy (bardzo słusznie!); dlatego format wydruku jest następujący:

...
row wiersz: kolumna wartość kolumna wartość kolumna wartość
....

Co jeszcze warto wiedzieć

MatSetValues(Mat A,int m,int *row,int n,int *col,Scalar *val,InsertMode mode)

Ta funkcja rozprasza w macierzy A macierz o rozmiarze mxn zapisaną w tablicy val w miejsca zapisane w tablicach row i col. Sposób działania tej funkcji wyjaśnia rysunek 1.

Aby wartości zostały fizycznie wpisane do macierzy należy, zawsze wywołać funkcje: MatAssemblyBegin(Mat A,MatAssemblyType type) i MatAssemblyEnd(Mat A, MatAssemblyType type). Zmienna type może przyjąć dwie wartości: MAT_FLUSH_ASSEMBLY bądź MAT_FINAL_ASSEMBLY. Pierwsza wartość jest używana, gdy zmieniamy tryb wstawiania wartości (INSERT_VALUES na ADD_VALUES, bądź na odwrót), druga - gdy kończymy wstawianie wartości i macierz ma zostać użyta. Jeżeli na macierzy nie została wykonana operacja FINAL_ASSEMBLY to próba jej wykorzystania w programie skończy się komunikatem o błędzie.

Funkcja MatCreate() tworzy szkielet macierzy rozrzedzonej pewnego (domyślnego) rodzaju. Użytkownik może też samemu wpłynąć na sposób generowania macierzy. Na przykład, można zaznaczyć jak duże kawałki macierzy mają być przechowywane w poszczególnych procesach, czy też, jak macierze mają być składowane: wierszami czy kolumnami. Szczegółowy opis można znaleźć w manualu [2] i man pages [3]. Dla potrzeb początkującego użytkownika przydatna jest jeszcze funkcja MatCompress(Mat A), która poddaje macierz kompresji, tak aby zajmowała jak najmniej miejsca w pamięci, gdyż funkcja MatCreate alokuje 10 razy ilość wierszy macierzy, co z reguły jest za dużo. Warto zatem po zakończeniu wpisywania wartości i wywołaniu funkcji MatAssemblyEnd(), wywołać funkcję MatCompress(). W przypadku gdy wiemy, że macierz zajmie więcej miejsca niż domyślne "10 razy ilość wierszy" trzeba użyć funkcji MatCreateSeqAIJ(), lub podobnej, ponieważ alokowanie dodatkowego miejsca w czasie wstawiania wartości do macierzy może spowodować nawet 50 krotne (!) spowolnienie działania funkcji.

Jak już wspomnieliśmy, dla macierzy istnieje funkcja MatView(Mat A,Viewer view), która pokazuje macierz w przeglądarce. Przeglądarka VIEWER_DRAWX_WORLD pokazuje, w okienku graficznym, strukturę niezerowych miejsc macierzy. Rysunek 2 pokazuje efekt przeglądania w trybie graficznym macierzy z poprzedniego przykładu.

Podstawowe operacje macierzowe

Nazwa funkcji	Operacja
MatEqual(Mat A,Mat B)	sprawdza czy A=B
MatMult(Mat A, Vec x, Vec y)
MatMultAdd(Mat A, Vec x,Vec y, Vec z)
MatMultTrans(Mat A, Vec x,Vec y)
MatMultTransAdd(Mat A, Vec x,Vec y, Vec z)

Podsumowanie

• Zadeklarować zmienną typu Mat
• Utworzyć strukturę macierzy funkcją MatCreate() (lub inną funkcją, pozwalającą na podanie przypuszczalnego współczynnika zapełnienia macierzy).
• Wstawić do macierzy wartości funkcją MatSetValues()
• Wywołać funkcje MatAssemblyBegin() i MatAssemblyEnd() (z opcją FINAL_ASSEMBLY)
• Gdy macierz nie jest już potrzebna, zwolnić pamięć funkcją MatDestroy()

Przykład podsumowujący

po czym mnożony jest przez nią wektor złożony z samych jedynek, a efekt zostaje wypisany na konsolę. Uważny czytelnik zapewne zauważy, że poprawna macierz dyskretyzacji laplasjanu

#include "sles.h"

int main(int argc,char **args)
{
  Vec     x,y;
  Mat     A; 
  int     i,k,l,n=4,m=4,flg,col[5];
  Scalar  val[5];

/* inicjalizacja */
  PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);

/* wektory przydatne przy demonstracji */
  VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,n*m,&x);
  VecDuplicate(x,&y);
  val[0]=1.0;
  VecSet(val,x);

/* Tworzenie macierzy laplasjanu na prostokacie mxn */
  MatCreate(MPI_COMM_WORLD,n*m,n*m,&A);
  val[0] = -1.0; val[1] = -1.0; val[2] = 4.0;
  val[3]=-1.0;val[4]=-1.0;
  for (i=0; i<n*m; i++ ) 
  {
    col[0]=i+1; col[1]=i-1; col[2]=i;
    col[3]=i+m; col[4]=i-m;
    l=0;k=5;
    if(i<m)k=4;
    if((i+1)>(n-1)*m){k=4;col[3]=i-m;}
    if((i%m)==0){l=1;k=k-1;col[1]=i+1;}
    if((i+1)%m==0)if(i!=0){l=1;k=k-1;}
    MatSetValues(A,1,&i,k,&col[l],&val[l],INSERT_VALUES);
  }

/* Składanie macierzy*/
  MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

/*Mnożenie wektora x przez macierz A (y=A*x)*/
  MatMult(A,x,y);

/*Wypisanie wektora y na konsolę*/
  VecView(y,VIEWER_STDOUT_SELF);

/* zakończenie */
  VecDestroy(x);
  MatDestroy(A);
  PetscFinalize();
  return(0);
 }

hydra:/usr/home/>p5.6
2
1
1
2
1
0
0
1
1
0
0
1
2
1
1
2
hydra:/usr/home/> _

Macierze powłokowe

jest macierzową reprezentacją (w standardowej bazie R³) następującego operatora:

Taka definicja jest szczególnie ważna w przypadku macierzy rozrzedzonych, w których powtarzają się elementy, gdyż dzięki temu nie musimy tracić miejsca w pamięci na powtarzające się wartości. Na przykład macierz powstała z dyskretyzacji jednowymiarowego laplasjanu ma (z dokładnością do stałej) postać:

Działanie tą macierzą na wektor można zapisać następująco:

,

gdzie

Autorzy pakietu przewidzieli możliwość definiowania macierzy w ten sposób. Użytkownik informuje pakiet jedynie o rozmiarze macierzy i podaje poszczególne funkcje realizujące operacje na tej macierzy (czy tą macierzą). Istotne jest to, że taka macierz jest traktowana tak samo jak inne rodzaje macierzy. Dopóki program nie usiłuje wykonać na nich operacji niewykonalnej, to znaczy nie zdefiniowanej przez użytkownika, to nie ma absolutnie żadnej różnicy w posługiwaniu się tą macierzą i macierzą "zwykłą". Macierze definiowane w duchu operatorowym zwane są powłokowymi (ang. shell matrix).

Mając zadeklarowaną zmienną typu Mat, możemy utworzyć ją w postaci macierzy powłokowej. Do tego służy funkcja:
MatShellCreate(MPI_Comm com,int m,int n,int M,int N,void *ctx,Mat *mat)

Pierwszy argument to kanał komunikacyjny MPI (np. MPI_COMM_WORLD), drugi i trzeci argument (m i n) to rozmiar lokalny macierzy, to znaczy rozmiar części macierzy zapisanej w procesorze wywołującym tę funkcję. Czwarty i piąty argument (M i N) to globalny rozmiar macierzy. Oczywiście, w wersji sekwencyjnej będziemy mieli n=N i m=M; ponadto (zazwyczaj) nasze macierze będą kwadratowe więc wtedy dodatkowo m=n=M=N. Szósty argument (ctx) to wskaźnik do kontekstu związanego z macierzą, lub PETSC_NULL, gdy nie korzystamy z kontekstu (idea kontekstu zostanie szczegółowo wytłumaczona w rozdziale "Koncept kontekstu"). Ostatni argument to wskaźnik do zmiennej typu Mat z którą będzie powiązana macierz powłokowa.

Do definiowania operacji, którą można wykonać na macierzy powłokowej służy funkcja:

MatShellSetOperation(Mat mat,MatOperation op, void *f)

Pierwszy argument tej funkcji to zmienna powiązana z obiektem typu Mat, dla którego definiujemy operację. Drugi argument to stała okreslająca rodzaj operacji, a trzeci to wskaźnik do funkcji wykonującej tą operację. Brzmi to groźnie, ale w praktyce jako f przekazujemy po prostu nazwę funkcji. Zdefiniować można 63 operacje, w poniższej tabeli podanych jest kilka ważniejszych:

MATOP_MULT

MATOP_EQUAL

MATOP_GET_SIZE

MATOP_DESTROY

MATOP_VIEW

MATOP_NORM

MATOP_GET_DIAGONAL

MATOP_DIAGONAL_SCALE

MATOP_SET_VALUES

Przykłady

Przykład trywialny

#include "sles.h"
#define N=100

/****************************************************
                Procedura mnozenia przez macierz
****************************************************/

int mult(Mat A,Vec x,Vec y)
{
 VecCopy(x,y);
 return(0);
}

/****************************************************
                Program glowny
****************************************************/
main(int argc,char **args)
{
 Mat A;
 Vec v,b;
 PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);

 /* ...pominiete linie kodu: tworzenie wektora b i przypisanie mu
    wartosci... */

 /* tworzenie macierzy powłokowej */
 MatCreateShell(MPI_COMM_WORLD,N,N,N,N,PETSC_NULL,&A);
 MatShellSetOperation(A,MATOP_MULT,mult);

 /* zastosowanie */
 MatMult(A,v,b); 
 PetscFinalize();
}

Dyskretny laplasjan

#include "sles.h"
#include <stdio.h>

/****************************************************
      Procedura mnozenia przez macierz
****************************************************/
int mult(Mat A,Vec x,Vec y)
{
 int i,N,*idx;
 Scalar *wx,*wy;

 VecGetSize(x,&N);
 VecGetArray(x,&wx);
 VecGetArray(y,&wy);
 for(i=0;i<N;i++)
 { 
  wy[i] = 2*wx[i];
  if(i>0) wy[i] -= wx[i-1];
  if(i<N) wy[i] -= wx[i+1];
 }
 VecRestoreArray(x,&wx);
 VecRestoreArray(y,&wy); 
}
/****************************************************
                    Program glowny
****************************************************/
int main(int argc,char **args)
{
 Mat A;
 int N=10;
 Scalar val;
 Vec x,y;

 PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);

 /* tworzenie macierzy powłokowej */
 MatCreateShell(MPI_COMM_WORLD,N,N,N,N,0,&A);
 MatShellSetOperation(A,MATOP_MULT,mult);

 /* Wektory przydatne przy prezentacji */
 VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,N,&x);
 VecDuplicate(x,&y);
 val=1; 
 VecSet(&val,x);
 
 /* uzycie macierzy powlokowej */
 MatMult(A,x,y);
 VecView(y,VIEWER_STDOUT_SELF);

 /* Zakończenie */
 VecDestroy(x); VecDestroy(y);
 MatDestroy(A);
 PetscFinalize();
 return(0);
}

hydra:/usr/home/> przykład
1
0
0
0
0
0
0
0
0
1
hydra:/usr/home/> _

Ćwiczenia

1. Uruchom program z przykładu trywialnego dla różnych N i dodaj przeglądarki wektorów, dzięki której będzie można obejrzeć wektor przed i po mnożeniu przez macierz. Proponujemy także zaimplementowanie innej macierzy diagonalnej w postaci macierzy powłokowej, np.:

2. Wygeneruj macierz NxN:


obejrzyj ją w okienku, a następnie na konsoli. Naucz się odczytywać te widoki.

3. Przeskaluj macierz A o czynnik , gdzie . Napisz procedurę GenLap1D(Mat A, int h) która generuje macierz NxN dyskretyzacji 1-wymiarowego Laplasjanu.
4. Zrób to samo co w poprzednim ćwiczeniu, ale macierz zaimplementuj jako macierz powłokową.
5. Wygeneruj macierz dyskretyzacji operatora


metodą różnic skończonych. Wskazówka: użyj funkcji GenLap1D() a następnie dodaj wartości pochodzące z dyskretyzacji pierwszej pochodnej.

6. Obejrzyj macierz z poprzedniego ćwiczenia, następnie transponuj ją i znów obejrzyj.
7. Napisz funkcję tworzącą macierz dyskretyzacji pierwszej pochodnej różnicą centralną


Użyj obu macierzy do dyskretyzacji operatora L z ćwiczenia E) w wersji



(Najwygodniej zdefiniować go jako macierz powłokową!)

8. Napisz procedurę rozwiązywania układu z daną macierzą, metodą Richardsona z preconditionerem Jacobiego. Użyj funkcji MatGetDiag() do wyciągnięcia diagonali, a całą procedurę zaimplementuj używając funkcji wektorowych VecAXPY(), VecScale(), VecPointwiseDivide().
9. Zaprogramuj własną przeglądarkę do macierzy/wektorów, używając procedur graficznych PETSc, tak, aby każde wywołanie otwierało nowe okno, utrzymujące się do końca pracy programu.
10. Procedura z ćwiczenia G wymaga pewnych obiektów pomocniczych wektorów roboczych do przechowywania wyników pośrednich. Zastanów się jak to zgrabnie zaimplementować. Zajrzyj do rozdziału "Koncept kontekstu".
11. Trudne zadanie, ale warto. Zaimplementuj procedurę generującą macierz sztywności dla dyskretyzacji równania eliptycznego drugiego rzędu na zadanej triangulacji MES i rozmaitych elementów skończonych. Oczywiście, zgodnie z prawidłami sztuki, najpierw skonstruuj lokalne macierze sztywności, które potem złożysz w globalną macierz sztywności w procesie assemblacji. Znakomicie pasuje tu funkcja MatSetValues().