Dotychczas zdobyta przez nas wiedza pozwala na uruchomienie "prawdziwego" programu używającego PETSc. Jak to zostało zapowiedziane, program będzie rozwiązywał układ równań powstałych z dyskretyzacji metodą różnic skończonych zadania Poissona na kwadracie jednostkowym:
przy czym rozważymy
przypadek gdy 
.
Jak wiemy, dyskretyzacja taka prowadzi do następującej postaci macierzowej (patrz [6]):
gdzie macierz A ma postać:
,
static char help[] = "Rozwiązuje pięciodiagonalny układ z użyciem SLES.\n\n";
#include "sles.h"
#include <stdio.h>
int main(int argc,char **args)
{
Vec x, b, u;
Mat A; /* macierz układu równań */
SLES sles; /* kontekst SLESu */
int ierr, i, n = 10, m = 10,col[5], its,flg,k,l;
Scalar minus_one = -1.0, one = 1.0, val[5],h;
double norm;
PetscInitialize(&argc,&args,PETSC_NULL,help);
OptionsGetInt(PETSC_NULL,"n",&n,&flg);
OptionsGetInt(PETSC_NULL,"m",&m,&flg);
/* Tworzenie wektorów */
VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,n*m,&x);
VecDuplicate(x,&b);
VecDuplicate(x,&u);
/* Wartości wektora prawej strony */
val[0]=1;
VecSet(val,b);
/* Tworzenie macierzy laplasjanu na prostokacie mxn */
MatCreate(MPI_COMM_WORLD,n*m,n*m,&A);
h=1/(m*n+1);
val[0] = -1.0/h; val[1] = -1.0/h; val[2] = 4.0/h;
val[3] = -1.0/h; val[4] = -1.0/h;
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);
}
MatSetValues(A,1,&i,k,&col[l],&val[l],INSERT_VALUES);
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
/* Tworzenie kontekstu SLES */
SLESCreate(MPI_COMM_WORLD,&sles);
/*Ustawienie operatorów */
SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);
SLESSetFromOptions(sles);
/* rozwiązanie układu liniowego */
SLESSolve(sles,b,x,&its);
/* Sprawdzenie błedu */
VecAXPY(&minus_one,u,x);
VecNorm(x,NORM_2,&norm);
PetscPrintf(MPI_COMM_WORLD,"Norma błędu %g, %d iteracji\n",norm,its);
/* Zwolnienie przestrzeni roboczej */
VecDestroy(x);VecDestroy(u);
VecDestroy(b);MatDestroy(A);
SLESDestroy(sles);
PetscFinalize();
return(0);
}
Po uruchomieniu mamy:
hydra:/usr/home/students/bojko/> przyklad Norma błędu 0.00893566, 65 iteracji hydra:/usr/home/students/bojko/>
 |
 |
 |