Wektory

Podstawy

...
Vec X;
...

Najprostszą funkcją służącą do wygenerowania wektora jest funkcja:

Pierwszy argument to informacja dla MPI (Message Passing Interface) o sposobie komunikacji między procesami. Nie warto zagłębiać się w subtelności z tym związane. Ważne jest, że z argumentem MPI_COMM_WORLD funkcja działa poprawnie zarówno na maszynach wieloprocesorowych, jak i jednoprocesorowych i we wszystkich podobnych przypadkach rozważanych w tym podręczniku będzie używana ta właśnie wartość. Drugi argument, n, jest liczbą elementów wektora przechowywanych w procesorze, który wywołał funkcję; możemy podać wartość PETSC_DECIDE, gdy chcemy aby pakiet ustawił tę wielkość automatycznie. PETSC_DECIDE jest idealną wartością dla n, gdy program nasz ma działać przede wszystkim w wersji sekwencyjnej. Trzeci argument, N, jest liczbą współrzędnych całego wektora (można wstawić PETSC_DETERMINE, wtedy zostanie ona określona na podstawie liczby procesorów i wartości n). Ostatnim argumentem jest wskaźnik do zmiennej typu wektor. Wektor będzie składał się ze współrzędnych typu Scalar.

...
VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,100,&X);
...

...
double x[100]; 
...

...
double *x;
...
x=(double*)calloc(100,sizeof(double));
...

Funkcja VecCreate() jest zaprojektowana tak, aby automatycznie rozproszyć wektor pomiędzy dostępne procesory (procesy). Pakiet sam zadecyduje czy wygenerować wektor sekwencyjny, czy równoległy, a jeśli równoległy, to ile elementów wektora będzie przechowywanych w poszczególnych procesorach. Dodajmy na marginesie, że dostępne są też dwie funkcje, dzięki którym użytkownik może zadecydować jakiego typu wektory zostaną wygenerowane: VecCreateSeq() lub VecCreateMPI(). Dzięki drugiej z nich można ręcznie ustawić jak wektor będzie rozproszony pomiędzy procesory.

Po wygenerowaniu wektora można przypisać wartości poszczególnym współrzędnym. Służą do tego dwie funkcje: VecSet(Scalar *val, Vec X) która nada wszystkim współrzędnym jednakową wartość val lub funkcja VecSetValue(Vec X, int idx, Scalar val, InsertMode mode). Pierwszy argument tej funkcji to wektor w który wstawiamy wartość, drugi to numer współrzędnej, trzeci wartość na tej współrzędnej. Czwarty argument to tryb wstawiania wartości, Możliwe są dwie wartości: ADD_VALUES - dodaj wartość do istniejącej, INSERT_VALUES - wstaw nową wartość.

Przykład

a na końcu wektor zostanie wypisany na konsolę.

#include "sles.h"
static char help[]="przyklad2\n";

int main(int argc,char **args)
{
  Vec     x; 
  Scalar val;
  int i,n=10;

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

  /* wygenerowanie wektora */
  VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,n,&x);
  
  /* przypisanie wartosci */
  for(i=0;i<n;i++)
  {
   if(i>0)val=(Scalar)1/i; else val=1;
   VecSetValue(x,i,val,INSERT_VALUES);
  }

  /* wypisuje wektor na konsolę */
  VecView(x,VIEWER_STDOUT_SELF);


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

hydra:/usr/home/bin/> przykład
1
1
0.5
0.333333
0.25
0.2
0.166667
0.142857
0.125
0.111111
hydra:/usr/home/bin/> _

Po zadeklarowaniu, wygenerowaniu i przypisaniu wartości wektorowi x, na zakończenie programu wypisujemy go na ekran przy użyciu przeglądarki do wektorów VecView().

Co jeszcze warto wiedzieć

Aby wygenerować kilka takich samych wektorów (na przykład wektor rozwiązania równania i wektor prawej strony - oba są tego samego rozmiaru) można użyć funkcji VecDuplicate(Vec stary,Vec *nowy) lub VecDuplicateVecs(Vec stary,int n, Vec **nowe) pierwsza z tych funkcji utworzy wektor nowy o takiej samej strukturze jak wektor stary, to znaczy o tym samym rozmiarze i tak samo rozproszony pomiędzy procesorami. Druga utworzy n nowych wektorów o takiej samej strukturze jak stary i zapisuje w tablicy wektorów nowe. Powielenia struktury wektora nie należy mylić z kopiowaniem jego wartości. Do tego służy funkcja VecCopy(Vec x,Vec y), kopiująca wartości z wektora x do y. Wektory x i y musza w tym wypadku mieć taką samą strukturę, co znaczy, że uprzednio muszą być utworzone, na przykład przy pomocy funkcji VecCreate() lub VecDuplicate().

Nie istnieje, niestety, funkcja VecGetValues(), która pozwalałaby na pobranie wartości z wektora. W przypadku równoległym trzeba stosować specjalne funkcje zbierające i rozsyłające dane po procesorach. W przypadku sekwencyjnym można używać funkcji VecGetArray(Vec x,Scalar **val), która zwraca wskaźnik do tablicy elementów wektora przechowywanych w procesorze który wywołał tę procedurę. Po zakończeniu wszelkich operacji na tablicy elementów trzeba ją z powrotem wstawić do wektora funkcją VecRestoreArray(Vec x, Scalar **val).

Przydatna jest też funkcja VecView(Vec v, Viewer viewer), która pokazuje wektor v przy pomocy przeglądarki viewer. Istnieją trzy predefiniowane przeglądarki:

Gdy wektor nie jest już potrzebny, należy zwolnić pamięć którą zajmuje funkcją VecDestroy(Vec x) lub - w przypadku tablicy wektorów - VecDestroyVecs(Vec *wektory,int n). 

Podstawowe operacje wektorowe

Nazwa funkcji	Operacja
VecAXPY(Scalar *a,Vec x,Vec y);
VecAYPX(Scalar *a,Vec x,Vec y);
VecWAXPY(Scalar *a,Vec x,Vec y,Vec w);
VecScale(Scalar *a,Vec x);
VecDot(Vec x,Vec y,Scalar *r);
VecTDot(Vec x,Vec y,Scalar *r);
VecNorm(Vec x,NormType norm,Double *r);  norm{NORM_1,NORM_2,NORM_INFINITY}
VecSum(Vec x, Scalar *r);
VecCopy(Vec x, Vec y);
VecSwap(Vec x, Vec y);
VecPointwiseMult(Vec x, Vec y, Vec w);
VecPointwiseDivide(Vec x, Vec y, Vec w);
VecMax(Vec x, int *idx, double *r);
VecMin(Vec x, int *idx, double *r);
VecAbs(Vec x);
VecReciprocal(Vec x);
VecShift(Scalar *s, Vec x);

Podsumowanie

• Zadeklarować zmienną typu Vec
• Utworzyć wektor funkcją VecCreate()
• Wstawić wartości funkcją VecSetValues(), po czym wywołać funkcje VecAssemblyBegin() i VecAssemblyEnd()
• Po użyciu wektora, zwolnić zajmowaną przezeń pamięć funkcją VecDestroy()

Przykład podsumowujący

static char help[]="przyklad 4.6";
#include "sles.h"
#include <stdio.h>

int main(int argc,char **args)
{
  Vec     x,y; 
  Scalar  val,*wartosci;
  int     i,n=5,*indeksy,flg;

  PetscInitialize(&argc,&args,PETSC_NULL,help);
  
  VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,n,&x);
  VecDuplicate(x,&y);
  indeksy=(int*)calloc(n/2,sizeof(int));
  wartosci=(Scalar*)calloc(n/2,sizeof(Scalar));
  val=0.5;
  VecSet(&val,x);
  VecView(x,VIEWER_STDOUT_SELF);
  PetscPrintf(MPI_COMM_WORLD,"\n");
  VecCopy(x,y);

/* przygotowanie danych */
  for(i=0;i<n/2;i++)
  {
   if(i>0)wartosci[i]=(Scalar)1/i; 
   else wartosci[i]=1;
   indeksy[i]=2*i;
  }

/* Wstawienie danych do wektora */
  VecSetValues(y,n/2,indeksy,wartosci,ADD_VALUES);
  VecAssemblyBegin(y);
  VecAssemblyEnd(y);
  VecView(y,VIEWER_STDOUT_SELF);
  PetscPrintf(MPI_COMM_WORLD,"\n");  

/* policzenie różnicy dwóch wektorów */
  val=-1;
  VecAXPY(&val,x,y);

/*wypisanie wektorów */
  VecView(y, VIEWER_STDOUT_SELF);
  PetscPrintf(MPI_COMM_WORLD,"\n");

/*zakończenie*/
  VecDestroy(x);
  VecDestroy(y);
  free(indeksy);
  free(wartosci);
  PetscFinalize();
  return(0);
 }

hydra:/usr/home/> przyklad
0.5
0.5
0.5
0.5
0.5

1.5
0.5
1.5
0.5
0.5

1
0
1
0
0

hydra:/usr/home/> _

Ćwiczenia

1. Utwórz na dwa sposoby (VecSetValue()+ pętla lub VecSet()) wektor o 10 współrzędnych wszystkich równych a) 1, b), Obejrzyj te wektory na konsoli i w okienku graficznym. (Uwaga: w pliku <math.h> jest zdefiniowana stała M_PI zawierająca przybliżenie )
2. Utwórz wektor x o N współrzędnych równych . (Zauważ, że taki wektor można potraktować jako dyskretyzację funkcji na odcinku [0,1]). Obejrzyj ten wektor na konsoli i w okienku graficznym. Eksperymentuj dla różnych wartości N.
3. Utwórz wektor y o NxN współrzędnych . Przez analogię z poprzednim można ten wektor traktować jako dyskretyzację funkcji na kwadracie [0,1]² z siatką o kroku h = 1/(N-1). Obejrzyj y w okienku graficznym. Zobacz jak trudno zgadnąć, jaką funkcję dyskretyzuje y, stąd wniosek, że do takiej wizualizacji trzeba używać innych narzędzi. Dowiedz się, jakie programy wizualizacyjne są dostępne na Twoim komputerze i jak z nich korzystać.
4. Oblicz normy a), b), c)wektorów x i y z poprzednich ćwiczeń. Jakich relacji między tymi normami można się spodziewać?
5. Zwróć uwagę, że dyskretna norma l₂ funkcji siatkowej reprezentowanej przez wektor x to odpowiednio przeskalowana norma euklidesowa tego wektora. Napisz procedurę, obliczającą normę l₂ wektora siatkowego x określonego na siatce o kroku h i N węzłach. Sprawdź na poprzednich przykładach, jak zachowuje się ta norma, gdy N rośnie. Porównaj to z normami: euklidesową i "max".
6. Wygeneruj wektor o długości 5 i losowych współrzędnych z przedziału [-1,1]. Będzie do tego potrzebna funkcja PetscRandom() której opis znajduje się w man pages [3].
7. Napisz funkcję VecGetValues(Vec v,int n,int *ind, Scalar *val), która z wektora wyciąga n kolejnych wartości.